A profound grasp of the quantitative spatial heterogeneity and distribution of the soil physicochemical attributes is crucial in understanding agricultural landscapes for ensuring the provisioning of soil ecosystem services. However, the analysis of data from remote sensing, like NDVI, can be of help in analysing the capacity of the landscape to provide supporting ecosystem services such as primary productivity. The research investigated and addressed the dispersion of important soil physico-chemical attributes in agricultural lands of the temperate Himalayan region of India using a geostatistical method and combining normalized difference vegetation index (NDVI) time-series data and the regression Kriging method. A 206 soil samples were gathered and assessed for soil parameters like pH, EC, OC, and available N, P, K, Ca, and Mg from Kishtwar district of Jammu. The coefficient of variation (CV) for pH and electrical conductivity (EC) ranged notably from 8.75 % to 118.98 %, highlighting diverse soil characteristics critical for local management practices. Mean elevation averaged 2743.32 m (m), with a moderate NDVI of 0.15, indicating dynamics in vegetation cover. Soil pH ranged from intensely acidic to marginally alkaline, with varying EC levels. Seemingly high organic carbon (OC), nitrogen (N), and potassium (K) levels, accompanied by medium phosphorus (P), calcium (Ca), and magnesium (Mg) levels were found in the region. The study employed ordinary kriging (OK) to map the spatial distribution of soil parameters, utilizing mean square error (MSE), root mean square error (RMSE), and the Moran’s I index. Exponential models were the best fit models for OC, while spherical models were fit for pH, EC, N, P, and Ca. Mathematical models were best fit for K and Mg. Spatial analysis using spherical and exponential models revealed distinct distribution patterns for pH, N, P, Ca, and Mg. The results of the degree of spatial dependence from the semi-variogram analyses indicated a strong (0.06 %) to moderate (0.51 %) to weak (2.81 %) dependence. The interpolated maps showed a distinct gradient in elevation (1053–4413 m), OC (0.13–2.80 %), NDVI (−0.16–0.54), pH (4.80–8.00), EC (0.03–9.80 dS m−1), N (201.15–993.19 kg ha−1), P (3.00–96.00 kg ha−1), K (124.88–1110.71 kg ha−1), Ca (7.00–46.00 meq 100 g soil−1), and Mg (2.30–21.50 meq 100 g soil−1) at the regional scale, indicating a wide range of spatial soil heterogeneity. The heterogeneity maps of soil parameters generated by this research can be effectively used by land planners and farm managers at a regional scale for crop nutrient management to reduce soil contamination risk. These maps serve as baseline materials and effective tools for suitable land management strategies such as conservation-effective tillage, integrated nutrient management, and organic farming based on the spatial distribution of soil properties and they can significantly enhance the long-term ecological sustainability of agro-ecosystems’ management.