Abstract: In the late 1990s, the "Dryland to Paddy" project was implemented in the Sanjiang Plain. After planting rice in the Albic soil, the barrier soil layer turns into a favorable soil layer, the low-yield soil becomes high-yield soil, and the Albic soil phosphorus pool increases. After flooding, the availability of phosphorus (closed storage phosphorusO-P and iron-bound phosphorusFe-P) increases with the decrease in soil redox potential (Eh) and the increase in pH, which substantially affects soil phosphorus heterogeneity. Therefore, we urgently need an optimal interpolation method to improve the prediction accuracy of total phosphorus in the Albic soil of typical irrigation areas of the Sanjiang Plain. This will help evaluate the impact of climate change and land use on the soil phosphorus pools and provide a reference for estimating future soil phosphorus pools. This study used the inverse distance weighting (IDW) method, radial basis function (RBF), ordinary Kriging (OK), global polynomial method (GPI), local polynomial method (LPI), geographic weighted regression (GWR), and geographic weighting regression to Kriging (GWRK) to predict the distribution of soil phosphorus in the Bawusan, Qiliqin, and Daxing irrigation areas of the Sanjiang Plain. The cross-validation method was used to obtain the mean error (ME), root mean square error (RMSE), and relative improvement (RI) to compare the accuracies of the various methods to determine the best interpolation method for assessing the spatial heterogeneity of phosphorus in the same soil type with different sampling densities. Based on the assumptions of regression analysis, this study incorporated 24 environmental variables for exploratory regression analysis, including elevation, pH, organic matter, exchangeable sodium, total nitrogen, available phosphorus, available copper, cultivated layer bulk density, and rice yield. According to the regression results, the auxiliary variables that were significantly correlated with phosphorus were selected for least square analysis. Finally, exchangeable sodium, cation exchange capacity, and available phosphorus were selected as auxiliary variables for the Bawusan irrigation area; organic matter, available zinc, and available boron were selected as auxiliary variables for the Qiliqin irrigation area; and cation exchange capacity, available zinc and copper were selected as auxiliary variables for the Daxing irrigation area. Compared to OK, RI indicated that the GWRK method with environmental auxiliary variables significantly improved the simulation accuracy of the spatial distribution of phosphorus. The IDW, GPI, and LPI methods reduced the accuracy of phosphorus spatial distribution simulation, whereas the RBF method was inconsistent. When comparing the mapping effect and interpolation speed of the seven interpolation methods, LPI, GPI, GWR, and GWRK had better mapping effects, whereas IDW, RBF, LPI, GPI, and OK were faster. The GWRK method had a better mapping effect, but it should be combined with environmental auxiliary variables, and the operation was complicated and slow. Sampling evenness also affected the prediction results. Nonetheless, GWRK had the lowest ME and RMSE, indicating that it is the best interpolation method. RBF is an optional method when the sampling evenness is lower. GWRK is the best interpolation method, but the results are affected by the number of auxiliary variables and collinearity between the variables.