Soil erosion is one of the most significant environmental problems in China. From 2010 to 2012, the fourth national census for soil erosion sampled 32 364 PSUs (Primary Sampling Units, small watersheds) with the areas of 0.2–3 km2. Land use and soil erosion controlling factors including rainfall erosivity, soil erodibility, slope length, slope steepness, biological practice, engineering practice, and tillage practice for the PSUs were surveyed, and the soil loss rate for each land use in the PSUs was estimated using an empirical model, the Chinese Soil Loss Equation (CSLE). Though the information collected from the sample units can be aggregated to estimate soil erosion conditions on a large scale; the problem of estimating soil erosion condition on a regional scale has not been addressed well. The aim of this study is to introduce a new model-based regional soil erosion assessment method combining a sample survey and geostatistics. We compared seven spatial interpolation models based on the bivariate penalized spline over triangulation (BPST) method to generate a regional soil erosion assessment from the PSUs. Shaanxi Province (3116 PSUs) in China was selected for the comparison and assessment as it is one of the areas with the most serious erosion problem. Ten-fold cross-validation based on the PSU data showed the model assisted by the land use, rainfall erosivity factor (R), soil erodibility factor (K), slope steepness factor (S), and slope length factor (L) derived from a 1 : 10 000 topography map is the best one, with the model efficiency coefficient (ME) being 0.75 and the MSE being 55.8 % of that for the model assisted by the land use alone. Among four erosion factors as the covariates, the S factor contributed the most information, followed by K and L factors, and R factor made almost no contribution to the spatial estimation of soil loss. The LS factor derived from 30 or 90 m Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) data worsened the estimation when used as the covariates for the interpolation of soil loss. Due to the unavailability of a 1 : 10 000 topography map for the entire area in this study, the model assisted by the land use, R, and K factors, with a resolution of 250 m, was used to generate the regional assessment of the soil erosion for Shaanxi Province. It demonstrated that 54.3 % of total land in Shaanxi Province had annual soil loss equal to or greater than 5 t ha−1 yr−1. High (20–40 t ha−1 yr−1), severe (40–80 t ha−1 yr−1), and extreme ( >  80 t ha−1 yr−1) erosion occupied 14.0 % of the total land. The dry land and irrigated land, forest, shrubland, and grassland in Shaanxi Province had mean soil loss rates of 21.77, 3.51, 10.00, and 7.27 t ha−1 yr−1, respectively. Annual soil loss was about 207.3 Mt in Shaanxi Province, with 68.9 % of soil loss originating from the farmlands and grasslands in Yan'an and Yulin districts in the northern Loess Plateau region and Ankang and Hanzhong districts in the southern Qingba mountainous region. This methodology provides a more accurate regional soil erosion assessment and can help policymakers to take effective measures to mediate soil erosion risks.