← Back to examples
Agent skill

Compute area-by-threshold

Compute the area within an AOI where a variable from a Cecil dataset crosses a threshold (e.g., "hectares of forest with biomass above 100 Mg/ha"). Reports per-timestep pixel counts, percentage, and an approximate hectare figure.

MITSource on GitHub

Use this skill when the user asks "how much area meets condition X?" — e.g. forest with biomass above a level, cropland with NDVI below a value, pixels of a specific land-cover class.

Prerequisites

Run the subscribe-and-load skill first to get:

  • a loaded ds: xarray.Dataset
  • the source aoi: AOI (for the aoi.hectares figure)

Steps

  1. Pick the variable and threshold. Choose variable_name from ds.data_vars and a threshold appropriate to that variable's units. Pick an operator matching the user's question.

    import numpy as np

variable_name = "..."             # e.g. "loss_year", "tree_cover"
threshold = ...                    # numeric
operator = np.greater_equal        # or np.greater, np.less, np.less_equal, np.equal, np.not_equal

  2. Mask fill values and apply the threshold per-timestep.

    da = ds[variable_name]
nodata = da.attrs.get("_FillValue")
if nodata is not None:
    da = da.where(da != nodata)

time_values = da["time"].values if "time" in da.dims else [None]

results = []
for t in time_values:
    slice_da = da.sel(time=t) if t is not None else da
    values = np.ravel(slice_da.values)
    valid = values[~np.isnan(values)]
    if valid.size == 0:
        continue
    matched = int(np.count_nonzero(operator(valid, threshold)))
    pct = matched / valid.size * 100.0
    results.append({
        "time": str(t) if t is not None else None,
        "valid_pixels": int(valid.size),
        "matched_pixels": matched,
        "percentage": pct,
        "approx_hectares": pct / 100.0 * aoi.hectares,
    })

  3. Inspect or plot results. Each entry summarises one timestep (or the whole dataset if there's no time dim).

Important constraints

  • Equal-area approximation. This treats every pixel as having identical area. That is reasonable for small AOIs but biases outwards from the equator. For high-latitude AOIs or AOIs spanning many degrees of latitude, weight by actual pixel area (via rioxarray and the dataset's CRS) instead of using aoi.hectares directly.
  • _FillValue masking is essential — counting no-data pixels as zero or as a numeric value silently distorts the percentage.
  • Categorical variables (e.g. land-cover class) are valid here — use np.equal with the integer class code as the threshold.

References