The spatial variation of soil test P (STP) in grassland soils is becoming important because of the use of STP as a basis for policies such as the recently EU-introduced Nitrate Directive. This research investigates the spatial variation of soil P in grazed grassland plots with a long-term (38 y) experiment. A total of 326 soil samples (including 14 samples from an adjacent grasswood buffer zone) were collected based on a 10 × 10 m 2 grid system. The samples were measured for STP and other nutrients. The results were analyzed using conventional statistics, geostatistics, and a geographic information system (GIS). Soil test P concentrations followed a lognormal distribution, with a median of 5.30 mg L -1 and a geometric mean of 5.35 mg L -1 . Statistically significant (p < 0.01) positive correlation between STP and pH was found. Spatial clusters and spatial outliers were detected using the local Moran's I index (a local indicator of spatial association) and were mapped using GIS. An obvious low-value spatial-cluster area was observed on the plots that received zero-P fertilizer application from 1968 to 1998 and a large high-value spatial-cluster area was found on the relatively high-P fertilizer application plots (15 kg ha -1 y -1 ). The local Moran's I index was also effective in detecting spatial outliers, especially at locations close to spatial-cluster areas. To obtain a reliable and stable spatial structure, semivariogram of soil-P data was produced after elimination of spatial outliers. A spherical model with a nugget effect was chosen to fit the experimental semivariogram. The spatial-distribution map of soil P was produced using the kriging interpolation method. The interpolated distribution map was dominated by medium STP values, ranging from 3 mg to 8 mg L -1 . An evidently low-P-value area was present in the upper side of the study area, as zero or short-term P fertilizer was applied on the plots. Meanwhile, high-P-value area was located mainly on the plots receiving 15 kg P ha -1 y -1 (for 38 y) as these plots accumulated excess P after a long-term P-fertilizer spreading. The high-or low-value patterns were in line with the spatial clusters. Geostatistics, combined with GIS and the local spatial autocorrelation index, provides a useful tool for analyzing the spatial variation in soil nutrients.