We analyze the accuracy of wall shear stress measurements in lattice Boltzmann simulations that are based on a voxel representation of the geometry and staircase approximation of boundaries. Such approximations are commonly used in the context of lattice Boltzmann simulations, because they favor the use of simple and highly efficient data structures. We show on several two-and three-dimensional simulations that this low-order approximation of the boundary affects the accuracy of wall shear stress measurements in areas directly adjacent to the wall. A few lattice nodes apart from the wall, the accuracy is however largely improved, and can be considered to be compatible with the overall accuracy of a simulation at a given coarseness level of the grid. This result is interpreted as a justification for the use of walls with staircase shape, even in simulations with high expectations regarding the level of accuracy. Furthermore, we propose a novel method for establishing the direction of the wall normal, a quantity which is required for the computation of the wall shear stress. With this method, the wall normal is computed from local data that is extracted from the results of the fluid flow simulation. Owing to the nature of the flow dynamics, which tends to smooth out the asperities of the wall, the information on the wall orientation obtained in this way is observed to be of high quality.