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.

In [1]:

import tempfile
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
from pyramids.dataset import Dataset

from earthlens import EarthLens

OUT_DIR = Path(tempfile.mkdtemp(prefix="soilgrids_", dir="."))
BBOX = dict(lat_lim=[51.0, 51.5], lon_lim=[5.0, 5.5])   # [south, north], [west, east]
NO_DATA = -32768   # SoilGrids int16 no-data sentinel

import tempfile from pathlib import Path import matplotlib.pyplot as plt import numpy as np from pyramids.dataset import Dataset from earthlens import EarthLens OUT_DIR = Path(tempfile.mkdtemp(prefix="soilgrids_", dir=".")) BBOX = dict(lat_lim=[51.0, 51.5], lon_lim=[5.0, 5.5]) # [south, north], [west, east] NO_DATA = -32768 # SoilGrids int16 no-data sentinel

2026-07-01 13:00:22 | INFO | pyramids.base.config | Logging is configured.

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.

In [2]:

paths = EarthLens(
    data_source="soilgrids",
    variables=["phh2o", "clay"],
    depths=["0-5cm"],
    quantiles=["mean"],
    path=str(OUT_DIR),
    **BBOX,
).download()

[p.name for p in paths]

paths = EarthLens( data_source="soilgrids", variables=["phh2o", "clay"], depths=["0-5cm"], quantiles=["mean"], path=str(OUT_DIR), **BBOX, ).download() [p.name for p in paths]

soilgrids:   0%|          | 0/2 [00:00<?, ?coverage/s]

2026-07-01 13:00:24.047 | INFO     | earthlens.soilgrids.backend:_fetch_one:341 - soilgrids phh2o_0-5cm_mean: WCS subset https://maps.isric.org/mapserv?map=/map/phh2o.map (values are scaled integers in pH*10; divide by 10 for pH)

soilgrids:  50%|█████     | 1/2 [00:07<00:07,  7.12s/coverage]

2026-07-01 13:00:31.163 | INFO     | earthlens.soilgrids.backend:_fetch_one:341 - soilgrids clay_0-5cm_mean: WCS subset https://maps.isric.org/mapserv?map=/map/clay.map (values are scaled integers in g/kg; divide by 10 for %)

soilgrids: 100%|██████████| 2/2 [00:17<00:00,  8.98s/coverage]

soilgrids: 100%|██████████| 2/2 [00:17<00:00,  8.70s/coverage]

2026-07-01 13:00:41.449 | INFO     | earthlens.soilgrids.backend:download:399 - soilgrids attribution: SoilGrids 2.0 (c) ISRIC - World Soil Information, licensed CC-BY 4.0; cite Poggio et al. (2021), SOIL 7, 217-240. https://soilgrids.org

Out[2]:

['phh2o_0-5cm_mean.tif', 'clay_0-5cm_mean.tif']

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.

In [3]:

ph = Dataset.read_file(str(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))

ph = Dataset.read_file(str(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))

Out[3]:

(4.0, 6.8)

The values land in a plausible topsoil range (roughly pH 4–7 for the Netherlands). Now map it.

In [4]:

fig, ax = plt.subplots(figsize=(6, 5))
im = ax.imshow(ph_real, cmap="RdYlGn", origin="upper")
ax.set_title("SoilGrids topsoil pH (H2O), 0–5 cm")
ax.set_xlabel("column")
ax.set_ylabel("row")
fig.colorbar(im, ax=ax, label="pH")
plt.tight_layout()
plt.show()

fig, ax = plt.subplots(figsize=(6, 5)) im = ax.imshow(ph_real, cmap="RdYlGn", origin="upper") ax.set_title("SoilGrids topsoil pH (H2O), 0–5 cm") ax.set_xlabel("column") ax.set_ylabel("row") fig.colorbar(im, ax=ax, label="pH") plt.tight_layout() plt.show()

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.

In [5]:

clay = Dataset.read_file(str(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 -> %

fig, ax = plt.subplots(figsize=(6, 5))
im = ax.imshow(clay_pct, cmap="YlOrBr", origin="upper")
ax.set_title("SoilGrids topsoil clay fraction, 0–5 cm")
fig.colorbar(im, ax=ax, label="clay (%)")
plt.tight_layout()
plt.show()

clay = Dataset.read_file(str(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 -> % fig, ax = plt.subplots(figsize=(6, 5)) im = ax.imshow(clay_pct, cmap="YlOrBr", origin="upper") ax.set_title("SoilGrids topsoil clay fraction, 0–5 cm") fig.colorbar(im, ax=ax, label="clay (%)") plt.tight_layout() plt.show()

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 -32768 no-data sentinel.
• The isric alias is equivalent to soilgrids.

Next: the catalog explorer lists every property, depth, quantile and unit.