This study focused on the geographical layout optimization of comprehensive observation stations for terrestrial surface natural resources including land, water, climate, and biology. The study was carried out following three steps according to the spatial statistic trinity. (1) The spatial distribution characteristics of terrestrial surface natural resources were analyzed, and the Eco-geographical Region System for China and the Normalized Difference Vegetation Index (NDVI) of the recent 10 years were adopted to characterize the spatially stratified heterogeneity and their spatial autocorrelation; (2) The point mean of the surface with stratified nonhomogeneity (P-MSN) was chosen as the inference method and its average estimation error variance was set as the objective function for the layout optimization; (3) Spatial simulated annealing was used to minimize the objective function to generate the geographical layout of comprehensive observation stations. The following conclusions were drawn: (1) Average NDVI of multiple years can characterize the spatial distribution characteristics of terrestrial surface natural resources; (2) P-MSN can adapt to the spatial distribution character of terrestrial surface natural resources and place dense stations in areas with large variance and sparse stations in areas with small variance; (3) The sample size-estimation error variance curve can be used to determine the best sample size and 1,000 stations are suggested in this study.