SoilGrids depth profile — how a request expands¶

A single request expands to the Cartesian product of the properties, depths and quantiles you ask for — one coverage (one GeoTIFF) per (property, depth, quantile) cell. This notebook fetches clay across all six standard depths and builds a depth profile: how clay content changes with depth for one area, using the scaled-integer → percentage conversion.

Setup and fetch every depth¶

Omitting depths= fetches every depth the property publishes; we keep the mean layer. For clay that is six coverages — one per standard depth interval.

In [ ]:

import tempfile
from pathlib import Path

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

from earthlens.core import EarthLens

OUT_DIR = Path(tempfile.mkdtemp(prefix="soilgrids_profile_", dir="."))
NO_DATA = -32768
DEPTHS = ["0-5cm", "5-15cm", "15-30cm", "30-60cm", "60-100cm", "100-200cm"]

paths = EarthLens(
    data_source="soilgrids",
    variables=["clay"],
    quantiles=["mean"],
    lat_lim=[51.0, 51.3],
    lon_lim=[5.0, 5.3],
    path=str(OUT_DIR),
).download()

[p.name for p in paths]

import tempfile from pathlib import Path import matplotlib.pyplot as plt import numpy as np from pyramids.dataset import Dataset from earthlens.core import EarthLens OUT_DIR = Path(tempfile.mkdtemp(prefix="soilgrids_profile_", dir=".")) NO_DATA = -32768 DEPTHS = ["0-5cm", "5-15cm", "15-30cm", "30-60cm", "60-100cm", "100-200cm"] paths = EarthLens( data_source="soilgrids", variables=["clay"], quantiles=["mean"], lat_lim=[51.0, 51.3], lon_lim=[5.0, 5.3], path=str(OUT_DIR), ).download() [p.name for p in paths]

The six file names carry the depth token, confirming the one-file-per-cell expansion: clay_0-5cm_mean.tif, clay_5-15cm_mean.tif, ….

Mean clay per depth¶

For each depth we read the raster, mask the no-data sentinel, convert to a percentage (clay scale_factor is 10), and take the area mean — one number per depth interval.

In [ ]:

def mean_clay_pct(depth: str) -> float:
    arr = (
        Dataset.read_file(str(OUT_DIR / f"clay_{depth}_mean.tif"))
        .read_array()
        .astype("float64")
    )
    arr[arr == NO_DATA] = np.nan
    return float(np.nanmean(arr) / 10.0)


profile = [mean_clay_pct(d) for d in DEPTHS]
profile

def mean_clay_pct(depth: str) -> float: arr = ( Dataset.read_file(str(OUT_DIR / f"clay_{depth}_mean.tif")) .read_array() .astype("float64") ) arr[arr == NO_DATA] = np.nan return float(np.nanmean(arr) / 10.0) profile = [mean_clay_pct(d) for d in DEPTHS] profile

Plot the profile¶

A depth profile is read top-to-bottom, so we invert the y-axis: shallow soil at the top, deep soil at the bottom. The midpoint depth of each interval is a convenient x/y coordinate.

In [ ]:

midpoints = [2.5, 10, 22.5, 45, 80, 150]  # cm, centre of each interval

fig, ax = plt.subplots(figsize=(5, 6))
ax.plot(profile, midpoints, marker="o")
ax.invert_yaxis()
ax.set_xlabel("mean clay content (%)")
ax.set_ylabel("depth (cm)")
ax.set_title("Clay content vs depth (SoilGrids mean)")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

midpoints = [2.5, 10, 22.5, 45, 80, 150] # cm, centre of each interval fig, ax = plt.subplots(figsize=(5, 6)) ax.plot(profile, midpoints, marker="o") ax.invert_yaxis() ax.set_xlabel("mean clay content (%)") ax.set_ylabel("depth (cm)") ax.set_title("Clay content vs depth (SoilGrids mean)") ax.grid(True, alpha=0.3) plt.tight_layout() plt.show()

Takeaway¶

• Leaving depths= unset fetches every standard depth; the request expands to one GeoTIFF per (property, depth, quantile).
• Reducing each depth's raster to an area mean (after the /scale_factor conversion) turns the stack into a depth profile.
• Swap variables=["clay"] for any other property, or add quantiles=["Q0.05", "Q0.95"] to bracket the profile with the prediction interval.