This skill covers retrieving, processing, and analysing soil data from global and
national databases. Topics include the SoilGrids 2.0 REST API, soil organic carbon
stock calculation, USDA texture classification, simple kriging interpolation, and
depth-profile visualisation.
SoilGrids 2.0 (ISRIC) provides global predictions at 250 m resolution for key soil
properties at six standard depth intervals (0-5, 5-15, 15-30, 30-60, 60-100,
100-200 cm).
lat: float,
lon: float,
properties: list = None,
timeout: int = 30,
"""
Query SoilGrids 2.0 for a single point and return a tidy DataFrame.
Parameters
----------
lat, lon : float — WGS84 coordinates
properties : list — soil properties to retrieve (default: all seven)
timeout : int — HTTP timeout in seconds
Returns
-------
pd.DataFrame with columns: property, depth, mean, uncertainty_5pct,
uncertainty_95pct, unit
"""
if
is
None
"lon"
"lat"
"property"
"depth"
"value"
"mean"
"uncertainty"
for
in
"properties"
"layers"
"name"
"unit_measure"
"mapped_units"
""
"unit_measure"
"conversion_factor"
1.0
for
in
"depths"
"label"
"values"
"property"
"depth"
"mean"
"mean"
1
if
else
1
"Q0.05"
"Q0.05"
1
if
else
1
"Q0.95"
"Q0.95"
1
if
else
1
"unit"
return
1.2 Bounding-Box Grid Query
defget_soilgrids_bbox(
bbox: tuple,
property: str = "soc",
depth: str = "0-5cm",
n_points: int = 25,
timeout: int = 30,
) -> pd.DataFrame:
"""
Sample a soil property on a regular grid within a bounding box.
Parameters
----------
bbox : (min_lon, min_lat, max_lon, max_lat)
property : str — SoilGrids property code
depth : str — depth label, e.g. '0-5cm'
n_points : int — approximate number of grid points (square root taken per axis)
Returns
-------
pd.DataFrame with columns: lat, lon, property, depth, mean, unit
"""
min_lon, min_lat, max_lon, max_lat = bbox
side = max(2, int(np.sqrt(n_points)))
lons = np.linspace(min_lon, max_lon, side)
lats = np.linspace(min_lat, max_lat, side)
records = []
total = side * side
for i, lat inenumerate(lats):
for j, lon inenumerate(lons):
try:
df_pt = get_soilgrids_point(lat, lon, properties=[property], timeout=timeout)
row = df_pt[df_pt["depth"] == depth]
ifnot row.empty:
records.append({
"lat": lat,
"lon": lon,
"property": property,
"depth": depth,
"mean": float(row["mean"].iloc[0]),
"unit": row["unit"].iloc[0],
})
print(f" [{i*side+j+1}/{total}] ({lat:.3f}, {lon:.3f}): "f"{property} = {records[-1]['mean']:.2f}")
except Exception as exc:
print(f" WARN ({lat:.3f}, {lon:.3f}): {exc}")
return pd.DataFrame(records)
2. Derived Calculations
2.1 Soil Organic Carbon Stock
The SOC stock (kg C m⁻²) for a single depth interval is:
# Install once:# install.packages(c("soilDB", "aqp", "sf", "ggplot2"))
library(soilDB)
library(aqp)
library(sf)## Fetch SSURGO data for a set of map unit keys
fetch_ssurgo_profiles <-function(mukeys){# Build SQL to get horizon data
q <- sprintf("SELECT cokey, hzname, hzdept_r, hzdepb_r, sandtotal_r,
silttotal_r, claytotal_r, om_r, ph1to1h2o_r, dbthirdbar_r
FROM chorizon
WHERE cokey IN (
SELECT cokey FROM component WHERE mukey IN (%s)
)
ORDER BY cokey, hzdept_r",
paste(mukeys, collapse =","))
hz <- SDA_query(q)return(hz)}## Compute SOC stock for SSURGO horizons
compute_ssurgo_soc <-function(hz_df){
hz_df$thickness_cm <- hz_df$hzdepb_r - hz_df$hzdept_r
hz_df$soc_pct <- hz_df$om_r /1.724# SOM to SOC conversion
hz_df$bd_gcm3 <- hz_df$dbthirdbar_r
hz_df$soc_stock <- with(hz_df,(soc_pct /100)* bd_gcm3 * thickness_cm *10)# kg C m⁻²
total_stock <- tapply(hz_df$soc_stock, hz_df$cokey,sum, na.rm =TRUE)return(total_stock)}# Example usage:# mukeys <- c("2494753", "2494754", "2494755")# hz_data <- fetch_ssurgo_profiles(mukeys)# soc_by_component <- compute_ssurgo_soc(hz_data)# print(soc_by_component)
6. Examples
Example A — Map SOC Stocks Across a Regional Grid
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# Define study area (e.g., part of the US Corn Belt)
BBOX = (-91.0, 41.0, -88.0, 43.0) # (min_lon, min_lat, max_lon, max_lat)# Step 1: Sample SOC at a coarse grid (25 points) via SoilGrids
df_grid = get_soilgrids_bbox(
bbox=BBOX,
property="soc",
depth="0-5cm",
n_points=25,
)
print(df_grid.head())
# Step 2: Krige to a finer grid
known_coords = df_grid[["lon", "lat"]].values
known_values = df_grid["mean"].values
# Build 10x10 target grid
lons = np.linspace(BBOX[0], BBOX[2], 10)
lats = np.linspace(BBOX[1], BBOX[3], 10)
lon_grid, lat_grid = np.meshgrid(lons, lats)
target_coords = np.column_stack([lon_grid.ravel(), lat_grid.ravel()])
soc_predicted = krige_ordinary(known_coords, known_values, target_coords)
soc_map = soc_predicted.reshape(10, 10)
# Step 3: Visualise
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
sc = axes[0].scatter(df_grid["lon"], df_grid["lat"], c=df_grid["mean"],
cmap="YlOrBr", s=80, edgecolors="k")
plt.colorbar(sc, ax=axes[0], label="SOC (g/kg)")
axes[0].set_title("Observed SOC (0-5 cm)")
axes[0].set_xlabel("Longitude"); axes[0].set_ylabel("Latitude")
im = axes[1].imshow(
soc_map, origin="lower", cmap="YlOrBr",
extent=[BBOX[0], BBOX[2], BBOX[1], BBOX[3]], aspect="auto",
)
plt.colorbar(im, ax=axes[1], label="SOC (g/kg)")
axes[1].set_title("Kriged SOC (0-5 cm)")
axes[1].set_xlabel("Longitude"); axes[1].set_ylabel("Latitude")
plt.tight_layout()
plt.savefig("/tmp/soc_map.png", dpi=150)
plt.show()
print("Saved map to /tmp/soc_map.png")
Example B — Compare Soil Texture Across Multiple Field Sites
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
# Field site coordinates
sites = [
{"name": "Site A — Iowa", "lat": 42.0308, "lon": -93.6319},
{"name": "Site B — Kansas", "lat": 38.9717, "lon": -95.2353},
{"name": "Site C — Nebraska", "lat": 41.2565, "lon": -95.9345},
{"name": "Site D — Illinois", "lat": 40.6331, "lon": -89.3985},
]
# Retrieve sand/silt/clay at 0-30 cm for each site
records = []
for site in sites:
df = get_soilgrids_point(
lat=site["lat"],
lon=site["lon"],
properties=["sand", "silt", "clay"],
)
for depth in ["0-5cm", "5-15cm", "15-30cm"]:
row = {"site": site["name"], "depth": depth}
for prop in ["sand", "silt", "clay"]:
val = df[(df["property"] == prop) & (df["depth"] == depth)]["mean"]
row[prop] = float(val.iloc[0]) / 10.0ifnot val.empty else np.nan
records.append(row)
df_tex = pd.DataFrame(records)
# Add texture classification
df_tex["texture_class"] = df_tex.apply(
lambda r: classify_texture(r["sand"], r["clay"])
ifnot (np.isnan(r["sand"]) or np.isnan(r["clay"])) else"unknown",
axis=1,
)
print(df_tex.to_string(index=False))
# Ternary-style bar chart: stacked sand/silt/clay per site × depth
fig, ax = plt.subplots(figsize=(12, 5))
bar_width = 0.2
n_depths = 3
depth_labels = ["0-5cm", "5-15cm", "15-30cm"]
colours = {"sand": "#f4c542", "silt": "#a07850", "clay": "#c0392b"}
x = np.arange(len(sites))
for d_idx, depth inenumerate(depth_labels):
sub = df_tex[df_tex["depth"] == depth].set_index("site")
offset = (d_idx - 1) * bar_width
bottom = np.zeros(len(sites))
for frac in ["sand", "silt", "clay"]:
vals = [sub.loc[s["name"], frac] if s["name"] in sub.index else0for s in sites]
ax.bar(x + offset, vals, bar_width * 0.9, bottom=bottom,
color=colours[frac], label=frac if d_idx == 0else"")
bottom += np.array(vals)
ax.set_xticks(x)
ax.set_xticklabels([s["name"] for s in sites], rotation=15, ha="right")
ax.set_ylabel("Percentage (%)")
ax.set_title("Soil Texture by Site and Depth (SoilGrids 2.0)")
handles = [mpatches.Patch(color=v, label=k) for k, v in colours.items()]
ax.legend(handles=handles, loc="upper right")
plt.tight_layout()
plt.savefig("/tmp/texture_comparison.png", dpi=150)
plt.show()
print("Saved to /tmp/texture_comparison.png")
df_tex.to_csv("/tmp/soil_texture_sites.csv", index=False)
7. Tips and Gotchas
SoilGrids units: Values returned by the API are in mapped units (e.g., bulk
density in cg/cm³, SOC in dg/kg). Always divide by the conversion_factor in the
response or check the unit_measure field before using values.
Rate limiting: ISRIC does not publish a strict rate limit, but space requests 1-2 s
apart for large grids to avoid HTTP 429 errors.
SSURGO coverage: SSURGO covers the conterminous US only. Use SoilGrids or
FAO/HWSD for global work.
Texture triangle edge cases: Many simplified triangle implementations misclassify
points near class boundaries. Validate against the official USDA NRCS chart.
Kriging assumptions: Ordinary kriging assumes second-order stationarity. Always
inspect the experimental variogram before trusting kriged maps.
SOC stock uncertainty: Bulk density is often the largest source of error. Where
available, use locally measured BD rather than SoilGrids estimates.
8. References
Poggio et al. (2021). SoilGrids 2.0: producing soil information for the globe with
quantified spatial uncertainty. SOIL, 7, 217-240.
