SoilGrids quickstart — soil properties over a bounding box¶
ISRIC SoilGrids 2.0 is a global
250 m machine-learning map of soil properties. The earthlens soilgrids
backend subsets any property server-side over an OGC Web Coverage Service (WCS)
and writes it to GeoTIFF through pyramids
— earthlens never touches a competing array stack.
By the end of this notebook you will be able to fetch a soil property for a bounding box, read the result back, convert its scaled-integer values to a physical unit, and map it. We use a small onshore window over the Netherlands so every fetch returns in seconds.
Setup¶
earthlens provides the unified EarthLens entry point; pyramids reads the
written GeoTIFF; downloads go to a temporary directory next to this notebook.
import tempfile
from pathlib import Path
import numpy as np
from pyramids.dataset import Dataset, GeoReference
from earthlens.core import EarthLens
OUT_DIR = Path(tempfile.mkdtemp(prefix="soilgrids_"))
NO_DATA = -32768 # SoilGrids int16 no-data sentinel
Fetch a soil property¶
A request names one or more properties (variables=) plus a bounding box.
With no depths= / quantiles=, the backend fetches every standard depth at
the mean layer — here we pin a single topsoil depth (0-5cm) and the mean
layer to keep it to one file per property. download() returns the written
GeoTIFF paths.
paths = EarthLens(
data_source="soilgrids",
variables=["phh2o", "clay"],
depths=["0-5cm"],
quantiles=["mean"],
path=OUT_DIR,
lat_lim=[51.0, 51.5], # [south, north]
lon_lim=[5.0, 5.5], # [west, east]
).download()
[p.name for p in paths]
One GeoTIFF is written per (property, depth, quantile) cell, named
<property>_<depth>_<quantile>.tif. Two properties × one depth × one layer →
two files.
Scaled-integer units¶
SoilGrids stores every property as an integer to save space; divide by a
per-property scale_factor to recover the conventional unit. The backend keeps
the raw integers in the GeoTIFF (it never rescales), so you apply the factor
yourself.
| property | stored as | ÷ factor | physical unit |
|---|---|---|---|
phh2o |
pH ×10 | 10 | pH |
clay |
g/kg | 10 | % |
Read phh2o back and convert a stored value to pH.
ph = Dataset.read_file(OUT_DIR / "phh2o_0-5cm_mean.tif")
ph_scaled = ph.read_array().astype("float64")
ph_scaled[ph_scaled == NO_DATA] = np.nan
ph_real = ph_scaled / 10.0 # phh2o scale_factor is 10
float(np.nanmin(ph_real)), float(np.nanmax(ph_real))
The values land in a plausible topsoil range (roughly pH 4–7 for the Netherlands). Now map it.
ph_grid = Dataset.from_array(
ph_real,
no_data_value=np.nan,
geo_ref=GeoReference(geo=ph.geotransform, epsg=ph.epsg),
)
glyph = ph_grid.plot(cmap="RdYlGn", title="SoilGrids topsoil pH (H2O), 0–5 cm")
glyph.cbar.set_label("pH")
A second property — clay content¶
The same read-and-rescale recipe works for any property; only the
scale_factor changes. Clay is stored in g/kg, so dividing by 10 gives a mass
percentage.
clay = Dataset.read_file(OUT_DIR / "clay_0-5cm_mean.tif")
clay_pct = clay.read_array().astype("float64")
clay_pct[clay_pct == NO_DATA] = np.nan
clay_pct = clay_pct / 10.0 # clay scale_factor is 10 -> %
clay_grid = Dataset.from_array(
clay_pct,
no_data_value=np.nan,
geo_ref=GeoReference(geo=clay.geotransform, epsg=clay.epsg),
)
glyph = clay_grid.plot(cmap="YlOrBr", title="SoilGrids topsoil clay fraction, 0–5 cm")
glyph.cbar.set_label("clay (%)")
Takeaway¶
EarthLens(data_source="soilgrids", variables=[...], depths=[...], quantiles=[...], lat_lim=, lon_lim=)fetches one GeoTIFF per(property, depth, quantile).- Values are scaled integers — divide by the property's
scale_factor(next notebook lists them all) and mask the-32768no-data sentinel. - The
isricalias is equivalent tosoilgrids.
Next: the catalog explorer lists every property, depth, quantile and unit.