Land subsidence in southwest Taiwan is mainly caused by ground-water withdrawal and a monitoring network of benchmarks has been established to observe this subsidence. Using results from differential leveling, the magnitudes of land subsidence from 1992 to 1994 were calculated for 50 benchmarks. The calculated maximum magnitude of land subsidence is 15.4 centimeters between these two years. Spatial interpolation is frequently required to provide information for assessing the distribution and extent of subsidence. In this study, geostatistics were applied to describe and characterize the spatial structure of sample data. Thus, an implementation procedure was proposed for the quality assessment of interpolation. The procedure includes geostatistical analysis, outlier detection, and comparison of different estimation techniques. The results show that the approach proposed in this study is effective to ensure the quality of interpolation.