Nutrient cycling, phosphorus retention, and protection of water quality are critical ecosystem services provided by watersheds. To quantify these services in complex, multi-use landscapes a spatially-explicit framework is required that allows to derive their spatial patterns. In this study we integrate geographic information systems, global positioning systems, advanced statistical and geostatistical methods, and field data to assess soil phosphorus heterogeneity at regional scale. Our goal was to compare various digital soil mapping methods to predict soil phosphorus across the Santa Fe River Watershed (SFRW) [3,585 sq km] in north-central Florida. A stratified random sampling based on land use and soil order combinations was used to collect soil samples in four layers (0 to 30, 30 to 60, 60 to 120, and 120 to 180cm) at 141 sites in 2003 and 2004. An univariate interpolation method (Ordinary Kriging) and multivariate methods (Regression Kriging and Co-kriging) were used to predict and map geospatial distributions of soil phosphorus (Mehlich-1 P; MP) and soil total phosphorus (TP) using ancillary spatial environmental datasets across the SFRW. The functional relationships between soil properties and different environmental landscape properties were identified. Results showed that multivariate methods produced better predictions of soil MP and TP across the watershed. This study showed that digital soil mapping can predict soil properties and their spatial variation across large regions. The soil P distribution and variation maps will be helpful for management of fields in the SFRW to optimize nutrient status and minimize adverse impacts on the environment.