Based on the relatively low density sampling data from China's second national soil survey in northeast China, we compared the erformance of four spatial interpolation methods, inverse distance weighting(IDW), radial basis function(RBF), ordinary Kriging(OK) and regression Kriging(RK), under seven sample sizes in generating digital map of soil total nitrogen content at large scale based on the software ArcGIS and GS⁺. Results showed that soil total nitrogen content was in the range of 0.08~21.48 g/kg with great variability. The Nugget effect showed a medium strong spatial autocorrelation of soil total nitrogen content in the region, and the range of spatial autocorrelation was greater than research on smaller scale in the same region; Spatial variability of soil total nitrogen content changed when sample size was less than 171, in this case, spatial structure and accuracy test of interpolation were unbelievable. All of the four compared methods predicted the spatial pattern of soil total nitrogen content well with decreasing from northeast to southwest. The accuracy of interpolation changed in the order of RK >OK >RBF >IDW. With incorporated auxiliary variables of soil cation exchange capacity, soil depth, soil pH value and annual mean air temperature, the RK improved accuracy by 19.40%, 18.50% and 16.15% than IDW, RBF and OK, respectively. It also exhibited more details on soil total nitrogen content variation at the areas with sparse sample points. It suggested that the RK is a potential spatial interpolation method to improve the soil mapping accuracy at large area with low density sample sites.