Agriculture — crop water demand and soil moisture¶

Agroclimatic studies typically pair inputs (precipitation, irrigation potential), demand (reference evapotranspiration $ET_0$), and soil status (root-zone moisture). ERA5-Land carries all three at 0.1° resolution — finer than the standard ERA5 product and tuned for land-surface processes.

Domain context. A simplified agronomic water balance for a rain-fed field over one growing season:

$$ \text{deficit}(t) = \sum (P - |ET|)_t $$

When the cumulative deficit is negative, the crop is drawing on stored soil moisture; when positive, soil moisture should be recharging. Pairing the deficit time series with soil_moisture_layer_1 (top 7 cm) gives an immediate visual check on the model's coupling.

Setup¶

All imports up front. earthlens provides the unified EarthLens entry point plus the ecmwf Catalog and AggregationConfig; pyramids provides Dataset for reading the downloaded GeoTIFFs. OUT is the local folder the ERA5-Land files land in.

In [ ]:

from pathlib import Path

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

from earthlens import AggregationConfig, EarthLens
from earthlens.ecmwf import Catalog

OUT = Path("data/era5-land-corn-belt")
OUT.mkdir(parents=True, exist_ok=True)

from pathlib import Path import matplotlib.pyplot as plt import numpy as np import pandas as pd from pyramids.dataset import Dataset from earthlens import AggregationConfig, EarthLens from earthlens.ecmwf import Catalog OUT = Path("data/era5-land-corn-belt") OUT.mkdir(parents=True, exist_ok=True)

Step 1 — confirm the variables and their flux flags¶

Soil moisture is a state field (volumetric water content, instantaneous). Precipitation and evaporation are flux fields (per-step accumulations). op="auto" routes each correctly without the user picking.

In [ ]:

cat = Catalog()
for code in (
    "total-precipitation",
    "total-evaporation",
    "2m-temperature",
    "volumetric-soil-water-layer-1",
):
    spec = cat.get_variable("reanalysis-era5-land", code)
    print(
        f"{code:35s}  nc={spec.nc_variable:6s}  units={spec.units:25s}  is_flux={spec.is_flux}"
    )

cat = Catalog() for code in ( "total-precipitation", "total-evaporation", "2m-temperature", "volumetric-soil-water-layer-1", ): spec = cat.get_variable("reanalysis-era5-land", code) print( f"{code:35s} nc={spec.nc_variable:6s} units={spec.units:25s} is_flux={spec.is_flux}" )

Step 2 — retrieve a growing season at monthly resolution¶

Box: a ~1° area in the central US corn belt (40°–41°N, 92°–91°W). Range: May–October 2022 — the heart of the corn-belt growing season. ERA5-Land's monthly-means product gives one value per month at 0.1° native grid. Mixed flux + state retrieve in one call.

Build the request¶

Describe the retrieval declaratively: the ecmwf backend, monthly cadence, the May–October 2022 window, the ERA5-Land monthly-means dataset, the four variables, and the corn-belt bounding box. Nothing is fetched yet.

In [ ]:

earthlens = EarthLens(
    data_source="ecmwf",
    cadence="monthly",
    start="2022-05-01",
    end="2022-10-01",
    dataset="reanalysis-era5-land-monthly-means",
    variables=[
        "total-precipitation",
        "total-evaporation",
        "2m-temperature",
        "volumetric-soil-water-layer-1",
    ],
    aoi=[-92.0, 40.0, -91.0, 41.0],
    path=str(OUT),
)

earthlens = EarthLens( data_source="ecmwf", cadence="monthly", start="2022-05-01", end="2022-10-01", dataset="reanalysis-era5-land-monthly-means", variables=[ "total-precipitation", "total-evaporation", "2m-temperature", "volumetric-soil-water-layer-1", ], aoi=[-92.0, 40.0, -91.0, 41.0], path=str(OUT), )

Download and aggregate¶

download() pulls each variable from the CDS and chains the aggregator in the same call: monthly (1MS) means at the ERA5-Land native 0.1° cell size, with op="auto" routing flux fields (precip, evaporation) and state fields (temperature, soil moisture) correctly.

In [ ]:

earthlens.download(aggregate=AggregationConfig(freq="1MS", op="auto", cell_size=0.1))

earthlens.download(aggregate=AggregationConfig(freq="1MS", op="auto", cell_size=0.1))

Step 3 — assemble the deficit time series¶

Per-month catchment-mean for each variable, then derive the cumulative $P - |ET|$ deficit and pair it with soil-moisture L1.

A helper to stack the monthly GeoTIFFs¶

The aggregator writes one GeoTIFF per variable per month under OUT/aggregated. stack() reads every month for a given variable into a single time-stacked array.

In [ ]:

agg = OUT / "aggregated"


def stack(cds_variable: str) -> np.ndarray:
    paths = sorted(agg.glob(f"{cds_variable}_1MS_*.tif"))
    return np.stack([Dataset.read_file(str(p)).read_array() for p in paths])

agg = OUT / "aggregated" def stack(cds_variable: str) -> np.ndarray: paths = sorted(agg.glob(f"{cds_variable}_1MS_*.tif")) return np.stack([Dataset.read_file(str(p)).read_array() for p in paths])

Catchment-mean each variable¶

Spatially average each stacked field to one value per month, converting to convenient units: precipitation and evaporation to mm (with the sign flip that makes |ET| positive), temperature to °C, soil moisture left as volumetric water content.

In [ ]:

precip = np.nanmean(stack("total_precipitation"), axis=(1, 2)) * 1000  # mm
et = -np.nanmean(stack("total_evaporation"), axis=(1, 2)) * 1000  # mm (sign flip)
t_air = np.nanmean(stack("2m_temperature"), axis=(1, 2)) - 273.15  # °C
soil_moisture = np.nanmean(
    stack("volumetric_soil_water_layer_1"), axis=(1, 2)
)  # m^3/m^3

precip = np.nanmean(stack("total_precipitation"), axis=(1, 2)) * 1000 # mm et = -np.nanmean(stack("total_evaporation"), axis=(1, 2)) * 1000 # mm (sign flip) t_air = np.nanmean(stack("2m_temperature"), axis=(1, 2)) - 273.15 # °C soil_moisture = np.nanmean( stack("volumetric_soil_water_layer_1"), axis=(1, 2) ) # m^3/m^3

Assemble the deficit table¶

Collect the monthly series into a DataFrame, deriving the per-month and cumulative $P - |ET|$ deficit alongside air temperature and topsoil moisture.

In [ ]:

months = pd.date_range("2022-05-01", periods=len(precip), freq="MS")
df = pd.DataFrame(
    {
        "P [mm]": precip.round(1),
        "|ET| [mm]": et.round(1),
        "P-|ET| [mm]": (precip - et).round(1),
        "deficit cum [mm]": (precip - et).cumsum().round(1),
        "T_2m [°C]": t_air.round(1),
        "SM_L1": soil_moisture.round(3),
    },
    index=months,
)
df

months = pd.date_range("2022-05-01", periods=len(precip), freq="MS") df = pd.DataFrame( { "P [mm]": precip.round(1), "|ET| [mm]": et.round(1), "P-|ET| [mm]": (precip - et).round(1), "deficit cum [mm]": (precip - et).cumsum().round(1), "T_2m [°C]": t_air.round(1), "SM_L1": soil_moisture.round(3), }, index=months, ) df

Step 4 — soil moisture vs cumulative water surplus¶

If the model's land surface scheme is internally consistent, soil moisture should track the cumulative $P - |ET|$ trend (with the soil's natural drainage timescale).

In [ ]:

fig, ax1 = plt.subplots(figsize=(9, 5))
ax1.bar(
    months,
    precip - et,
    width=20,
    color="tab:green",
    alpha=0.5,
    label="P - |ET| [mm/month]",
)
ax1.set_ylabel("P - |ET| per month [mm]")
ax1.axhline(0, color="gray", lw=0.5)

ax2 = ax1.twinx()
ax2.plot(
    months,
    soil_moisture,
    marker="o",
    color="tab:brown",
    label="Soil moisture L1 [m³/m³]",
)
ax2.set_ylabel("Soil moisture L1 [m³/m³]", color="tab:brown")

ax1.set_title("Corn-belt growing season 2022 — water surplus and topsoil moisture")
fig.legend(loc="upper right", bbox_to_anchor=(0.95, 0.95))
plt.tight_layout()
plt.show()

fig, ax1 = plt.subplots(figsize=(9, 5)) ax1.bar( months, precip - et, width=20, color="tab:green", alpha=0.5, label="P - |ET| [mm/month]", ) ax1.set_ylabel("P - |ET| per month [mm]") ax1.axhline(0, color="gray", lw=0.5) ax2 = ax1.twinx() ax2.plot( months, soil_moisture, marker="o", color="tab:brown", label="Soil moisture L1 [m³/m³]", ) ax2.set_ylabel("Soil moisture L1 [m³/m³]", color="tab:brown") ax1.set_title("Corn-belt growing season 2022 — water surplus and topsoil moisture") fig.legend(loc="upper right", bbox_to_anchor=(0.95, 0.95)) plt.tight_layout() plt.show()

Notes¶

• reanalysis-era5-land is the ERA5 land-surface reanalysis at 0.1° resolution — the right product for agronomic studies. ERA5 single-levels lives at 0.25° and runs the same land scheme but with coarser surface fields.
• Reference ET ($ET_0$). ERA5 reports actual evaporation (model-computed). For Penman–Monteith $ET_0$, retrieve surface-net-solar-radiation, 2m-dewpoint-temperature, 10m-wind-speed, surface-pressure and run FAO-56 yourself; CDS doesn't expose $ET_0$ directly.
• Volumetric soil moisture layers. Layer 1 is 0–7 cm; layers 2/3/4 are 7–28, 28–100, 100–289 cm. Pull all four for a root-zone profile.
• Cell size. This notebook passes cell_size=0.1 (ERA5-Land native) on both EarthLens and AggregationConfig so the GeoTIFFs carry the right pixel-size tag in their metadata.