Abstract: Soil sampling design is an important issue for agricultural management and environmental monitoring. In this study, a total of 537 soil samples were collected based on a 30 × 30–m grid from a permanent dairy farm in southeastern Ireland. Five different subsample experiments at lower densities based on the original data set were performed to study the optimal soil sampling design for soil P and Mg using geostatistics and a GIS (geographical information system). Soil P ranged from 1.3 to 35.7 mg L−1, with a CV value of 0.68. Soil Mg ranged from 134.7 to 685.2 mg L−1, with a small CV value of 0.28. Soil P followed neither a normal nor a lognormal distribution. Box-Cox transformation was applied to achieve normality. On the other hand, soil Mg followed a normal distribution, as did its subdata. For soil P, an omnidirectional spherical model was used to describe the spatial autocorrelation. For soil Mg, a nested model (an exponential model combined with a Gaussian model) was used to fit the variograms. Further soil P interpolated maps revealed that soil grid sampling interval could increase to 90 m without a significant loss of spatial information, whereas soil Mg sampling interval could increase to 120 m, confirming that soil Mg had much stronger spatial structure than soil P. According to this study, a grid of 90 × 90 m was recommended for soil sampling, which was confirmed in other practical grassland farms. The spatial structure information was very useful to optimize soil sampling design for practical grassland management.