The estimation of regional gross primary production (GPP) is a crucial issue in carbon cycle studies. One commonly used way to estimate the characteristics of GPP is to infer the total amount of GPP by collecting field samples. In this process, the spatial sampling design will affect the error variance of GPP estimation. This letter uses geostatistical model-based sampling to optimize the sampling locations in a spatial heterogeneous area. The approach is illustrated with a real-world application of designing a sampling strategy for estimating the regional GPP in the Babao river basin, China. By considering the heterogeneities in the spatial distribution of the GPP, the sampling locations were optimized by minimizing the spatially averaged interpolation error variance. To accelerate the optimization process, a spatial simulated annealing search algorithm was employed. Compared with a sampling design without considering stratification and anisotropies, the proposed sampling method reduced the error variance of regional GPP estimation.