Two geostatistical interpolation methods were studied for rock joint modeling; ordinary kriging and indicator kriging. Geostatistics was extended to improve the spatial interpretation of joint parameters, especially for pole vectors, which was kriged on the sphere. A matrix approach was introduced for probabilistic block failure analysis, and applied to study the stability of pit slopes and subway tunnels. The localized structural stability was achieved in probabilistic terms based on the geostatistical joint model.