Three-dimensional surface geometry of rock discontinuities and its evolution with shearing are of great importance in understanding the deformability and hydro-mechanical behavior of rock masses. In the present research, surfaces of three natural rock fractures were digitized and studied before and after the direct shear test. The variography analysis of the surfaces indicated a strong non-linear trend in the data. Therefore, the spatial variability of rock fracture surfaces was decomposed to one deterministic component characterized by a base polynomial function, and one stochastic component described by the variogram of residuals. By using an image-processing technique, 343 damaged zones with different sizes, shapes, initial roughness characteristics, local stress fields, and asperity strength values were spatially located and clustered. In order to characterize the overall spatial structure of the degraded zones, the concept of ‘pseudo-zonal variogram’ was introduced. The results showed that the spatial continuity at the damage locations increased due to asperity degradation. The increase in the variogram range was anisotropic and tended to be higher in the shear direction; thus, the direction of maximum continuity rotated towards the shear direction. Finally, the regression-kriging method was used to reconstruct the morphology of the intact surfaces and degraded areas. The cross-validation error of interpolation for the damaged zones was found smaller than that obtained for the intact surface