generate

historical_daily

drop_noncore_coords(ds: xr.Dataset) -> xr.Dataset

Drop coordinates the older and newer CDS extract formats disagree about.

Extracts pulled through the newer CDS API carry number (ensemble member) and expver ('0001' final ERA5, '0005' preliminary ERA5T) alongside the grid coords; older extracts carry neither, and are packed int16 rather than float32. xr.concat refuses to join datasets whose coordinates differ, so a look-ahead that crosses from an old extract into a new one has to be normalised first.

Which era a file belongs to tracks its download date, not its data year: surveying the archive, valid_time covers 1950-1989 and 2024 while time covers 1990-2023. So there are two seams, not the single 2023/2024 one this used to describe -- 1989->1990 crosses new format into old, 2023->2024 crosses old into new -- and the 1950-2023 regeneration passed through both.

Source code in src/climate_data/generate/historical_daily.py

def drop_noncore_coords(ds: xr.Dataset) -> xr.Dataset:
    """Drop coordinates the older and newer CDS extract formats disagree about.

    Extracts pulled through the newer CDS API carry `number` (ensemble member) and
    `expver` ('0001' final ERA5, '0005' preliminary ERA5T) alongside the grid coords;
    older extracts carry neither, and are packed int16 rather than float32. `xr.concat`
    refuses to join datasets whose coordinates differ, so a look-ahead that crosses from
    an old extract into a new one has to be normalised first.

    Which era a file belongs to tracks its *download* date, not its data year: surveying
    the archive, `valid_time` covers 1950-1989 and 2024 while `time` covers 1990-2023. So
    there are two seams, not the single 2023/2024 one this used to describe -- 1989->1990
    crosses new format into old, 2023->2024 crosses old into new -- and the 1950-2023
    regeneration passed through both.
    """
    extra = []
    for coord in ds.coords:
        if coord not in CORE_COORDS:
            extra.append(coord)
    return ds.drop_vars(extra)

load_variable_with_lookahead(cdata: ClimateData, variable: str, year: str, month: str, dataset: str) -> xr.Dataset

Load a month's hourly data plus the one sample that closes its final day.

Accumulation windows are stamped by their end, so the last day of a month is closed by the next month's first sample -- for December, by the next year's January. That single extra sample completes the final bin without opening a new one.

Raises if the look-ahead file is absent rather than degrading to the incomplete 23-hour window, which should stop the run rather than quietly shorten one day. Checked before loading the target month so the run fails fast. Note this reaches one year past the history range -- regenerating the last history year needs the following January -- which is why the extract tasks span cdc.EXTRACT_YEARS rather than cdc.HISTORY_YEARS.

Source code in src/climate_data/generate/historical_daily.py

def load_variable_with_lookahead(
    cdata: ClimateData,
    variable: str,
    year: str,
    month: str,
    dataset: str,
) -> xr.Dataset:
    """Load a month's hourly data plus the one sample that closes its final day.

    Accumulation windows are stamped by their end, so the last day of a month is closed
    by the next month's first sample -- for December, by the next *year's* January. That
    single extra sample completes the final bin without opening a new one.

    Raises if the look-ahead file is absent rather than degrading to the incomplete
    23-hour window, which should stop the run rather than quietly shorten one day.
    Checked before loading the target month so the run fails fast. Note this reaches one
    year past the history range -- regenerating the last history year needs the following
    January -- which is why the extract tasks span `cdc.EXTRACT_YEARS` rather than
    `cdc.HISTORY_YEARS`.
    """
    next_year, next_month = _next_month(year, int(month))
    lookahead_path = cdata.extracted_era5_path(dataset, variable, next_year, next_month)
    if not lookahead_path.exists():
        msg = (
            f"Cannot close {year}-{month} for {variable}: the accumulation window of its"
            f" final day ends in the following month, whose extract is missing at"
            f" {lookahead_path}. Extract it before generating {year}:\n"
            f"  cdtask extract era5_download -d {dataset} -x {variable}"
            f" -y {next_year} -m {next_month}\n"
            f"  cdtask extract era5_compress -d {dataset} -x {variable}"
            f" -y {next_year} -m {next_month}\n"
            f"Use the single-job tasks, not `cdrun extract era5`: the runner's --year"
            f" stops at the last history year, and it decides what to fetch by file"
            f" existence."
        )
        raise FileNotFoundError(msg)

    ds = drop_noncore_coords(load_variable(cdata, variable, year, month, dataset))
    lookahead = drop_noncore_coords(
        load_variable(cdata, variable, next_year, next_month, dataset).isel(time=[0])
    )

    # Taking sample zero assumes the next month opens on midnight, which is what closes
    # the target month's final day. ERA5-Land 1950_01 opens at 01:00 instead, so the
    # assumption does fail somewhere in the archive -- silently, at one day per month.
    stamp = pd.Timestamp(lookahead.time.to_index()[0])
    if (stamp.hour, stamp.minute) != (0, 0):
        msg = (
            f"Cannot close {year}-{month} for {variable}: the look-ahead at"
            f" {lookahead_path} opens at {stamp}, not midnight, so its first sample does"
            f" not close the final day of {year}-{month}."
        )
        raise ValueError(msg)

    # `join="exact"` rather than the default `"outer"`: a look-ahead on a different grid
    # should fail here, not silently NaN-pad both months onto a union grid. Note this
    # guards the single-level datasets only -- `load_variable` overwrites ERA5-Land's
    # coordinates from `cdc.ERA5_LAND_*`, so a genuine land-grid change would be relabelled
    # before reaching this point. The real land files differ across the CDS format change
    # by ~1e-5 degrees, which is why that overwrite exists.
    return xr.concat([ds, lookahead], dim="time", join="exact")

trim_to_month(ds: xr.Dataset, year: str, month: int) -> xr.Dataset

Drop collapsed bins that fall outside the target month.

The interval-aware collapse labels its first bin with the previous month's last day, since that bin holds that day's closing sample. Left in place, every month contributes a duplicate date at its seam.

Source code in src/climate_data/generate/historical_daily.py

def trim_to_month(ds: xr.Dataset, year: str, month: int) -> xr.Dataset:
    """Drop collapsed bins that fall outside the target month.

    The interval-aware collapse labels its first bin with the *previous* month's last
    day, since that bin holds that day's closing sample. Left in place, every month
    contributes a duplicate date at its seam.
    """
    dates = pd.to_datetime(ds.date.to_index())
    keep = np.flatnonzero((dates.year == int(year)) & (dates.month == month))
    return ds.isel(date=keep)

scenario_annual

bounded_source_variables(target_variable: str, anomaly_scheme: str) -> list[str]

Source variables of target_variable that anomaly_scheme cannot bias-correct.

Value-bounded variables are floored and clipped to cdc.VALUE_BOUNDS instead of being stabilised, so they run under cdc.BOUNDED_VARIABLE_SCHEMES only; see check_bounded_variable_scheme for why the other schemes are refused. The bounds are keyed on the daily variable, so the annual target has to be resolved through its sources the same way the de-bias check is.

Source code in src/climate_data/generate/scenario_annual.py

def bounded_source_variables(target_variable: str, anomaly_scheme: str) -> list[str]:
    """Source variables of `target_variable` that `anomaly_scheme` cannot bias-correct.

    Value-bounded variables are floored and clipped to `cdc.VALUE_BOUNDS` instead of being
    stabilised, so they run under `cdc.BOUNDED_VARIABLE_SCHEMES` only; see
    `check_bounded_variable_scheme` for why the other schemes are refused. The bounds are
    keyed on the daily variable, so the annual target has to be resolved through its
    sources the same way the de-bias check is.
    """
    if anomaly_scheme in cdc.BOUNDED_VARIABLE_SCHEMES:
        return []
    bounded = []
    for source_variable in TRANSFORM_MAP[target_variable].source_variables:
        if source_variable in cdc.VALUE_BOUNDS:
            bounded.append(source_variable)
    return bounded

forecast_jobs_for_anomaly_scheme(to_run: list[tuple[str, str, str, str]], anomaly_scheme: str) -> list[tuple[str, str, str, str]]

Drop forecast jobs whose variable the anomaly scheme cannot be applied to.

Filtered per job rather than per variable because the historical scenario never reaches generate_scenario_daily_main -- it reads the daily results off disk -- so the scheme does not constrain it, and an additive variable is perfectly runnable there. Only the forecast scenarios pass through compute_anomaly.

Two classes are dropped, and they are not mirror images. Additive variables run under the stabilised monthly scheme and no other. Value-bounded variables are the reverse: they run under cdc.BOUNDED_VARIABLE_SCHEMES only, which means they have to be dropped under the default scheme too -- this stage calls generate_scenario_daily_main in memory, so a bounded forecast job reaches check_bounded_variable_scheme and raises in the worker after being scheduled. That is the failure the daily launcher's filter already prevents.

A filter that removed every job is a usage error rather than an empty run, because the caller asked for work that cannot be done and an empty fan-out looks like success.

Source code in src/climate_data/generate/scenario_annual.py

def forecast_jobs_for_anomaly_scheme(
    to_run: list[tuple[str, str, str, str]],
    anomaly_scheme: str,
) -> list[tuple[str, str, str, str]]:
    """Drop forecast jobs whose variable the anomaly scheme cannot be applied to.

    Filtered per job rather than per variable because the `historical` scenario never
    reaches `generate_scenario_daily_main` -- it reads the daily results off disk -- so
    the scheme does not constrain it, and an additive variable is perfectly runnable
    there. Only the forecast scenarios pass through `compute_anomaly`.

    Two classes are dropped, and they are not mirror images. Additive variables run under
    the stabilised `monthly` scheme and no other. Value-bounded variables are the reverse:
    they run under `cdc.BOUNDED_VARIABLE_SCHEMES` only, which means they have to be dropped
    under the default scheme too -- this stage calls `generate_scenario_daily_main` in
    memory, so a bounded forecast job reaches `check_bounded_variable_scheme` and raises in
    the worker after being scheduled. That is the failure the daily launcher's filter
    already prevents.

    A filter that removed every job is a usage error rather than an empty run, because the
    caller asked for work that cannot be done and an empty fan-out looks like success.
    """
    keep = []
    skipped_additive = set()
    skipped_bounded = set()
    dropped_additive = 0
    dropped_bounded = 0
    for job in to_run:
        variable, scenario = job[0], job[1]
        forecast = scenario != "historical"
        if forecast and bounded_source_variables(variable, anomaly_scheme):
            skipped_bounded.add(variable)
            dropped_bounded += 1
        elif (
            forecast
            and anomaly_scheme != cdc.ANOMALY_SCHEME_MONTHLY
            and ANOMALY_TYPES[variable] != "multiplicative"
        ):
            skipped_additive.add(variable)
            dropped_additive += 1
        else:
            keep.append(job)

    if skipped_additive:
        print(
            f"Anomaly scheme '{anomaly_scheme}' applies to multiplicative variables only;"
            f" skipping {dropped_additive} forecast tasks for:"
            f" {', '.join(sorted(skipped_additive))}."
        )
    if skipped_bounded:
        print(
            f"Anomaly scheme '{anomaly_scheme}' does not apply to value-bounded variables,"
            f" which run under {', '.join(cdc.BOUNDED_VARIABLE_SCHEMES)} only;"
            f" skipping {dropped_bounded} forecast tasks for:"
            f" {', '.join(sorted(skipped_bounded))}."
        )
    if to_run and not keep:
        msg = (
            f"Anomaly scheme '{anomaly_scheme}' can run none of the selected variables, so"
            f" there is nothing to submit. Skipped:"
            f" {', '.join(sorted(skipped_additive | skipped_bounded))}."
        )
        if skipped_bounded:
            msg += (
                f" Value-bounded variables need"
                f" --anomaly-scheme {cdc.BOUNDED_VARIABLE_SCHEMES[0]}."
            )
        raise click.UsageError(msg)
    return keep

scenario_daily

apply_dry_day_rule(anomaly: xr.Dataset, target: xr.Dataset, dry_day_rule: str) -> xr.Dataset

Stop the eps offset from manufacturing rain on days the model reports as dry.

The multiplicative anomaly (T + eps)/(R + eps) is strictly positive even when the model reports no rain at all, so a rainless model day still receives E_ref(month) * a(d) > 0 of the ERA5 climatology. Worse, eps dominates the numerator for such a day, so every rainless day in a cell-month gets the identical anomaly 1/(R_m + eps) -- a flat positive floor rather than a dry spell. Wherever that floor clears the 0.1 mm/day wet-day threshold the pipeline reports a wet day that neither ERA5 nor the GCM has, which is what makes precipitation_days step at the boundary.

preserve zeroes the anomaly on those days and rescales the surviving days of the same cell-month so the month's summed anomaly is unchanged. Because that sum is preserved per GCM cell and interpolate_to_target_latlon is linear, the monthly -- and therefore the annual -- total on the target grid is untouched to floating point. Only the distribution across days moves. That is deliberate: this is a shape fix, exactly orthogonal to the Jensen de-bias, which is a level fix that leaves the shape alone. The two commute, because the de-bias scales a whole month uniformly and this rescale is scale-invariant.

A cell-month the model reports dry on every day has nothing to renormalise onto. Those are left exactly as they are rather than zeroed. Zeroing them is the variant that was measured and rejected: it loses up to 1.2% of the population-weighted annual total across the 1-3% of cell-months that are all-dry. Keeping them is what makes this rule total-preserving, and it is the whole difference between the two.

Source code in src/climate_data/generate/scenario_daily.py

def apply_dry_day_rule(
    anomaly: xr.Dataset,
    target: xr.Dataset,
    dry_day_rule: str,
) -> xr.Dataset:
    """Stop the ``eps`` offset from manufacturing rain on days the model reports as dry.

    The multiplicative anomaly ``(T + eps)/(R + eps)`` is strictly positive even when the
    model reports no rain at all, so a rainless model day still receives
    ``E_ref(month) * a(d) > 0`` of the ERA5 climatology. Worse, ``eps`` dominates the
    numerator for such a day, so *every* rainless day in a cell-month gets the identical
    anomaly ``1/(R_m + eps)`` -- a flat positive floor rather than a dry spell. Wherever that
    floor clears the 0.1 mm/day wet-day threshold the pipeline reports a wet day that neither
    ERA5 nor the GCM has, which is what makes ``precipitation_days`` step at the boundary.

    ``preserve`` zeroes the anomaly on those days and rescales the surviving days of the same
    cell-month so the month's *summed* anomaly is unchanged. Because that sum is preserved per
    GCM cell and ``interpolate_to_target_latlon`` is linear, the monthly -- and therefore the
    annual -- total on the target grid is untouched to floating point. Only the distribution
    across days moves. That is deliberate: this is a shape fix, exactly orthogonal to the
    Jensen de-bias, which is a level fix that leaves the shape alone. The two commute, because
    the de-bias scales a whole month uniformly and this rescale is scale-invariant.

    A cell-month the model reports dry on *every* day has nothing to renormalise onto. Those
    are left exactly as they are rather than zeroed. Zeroing them is the variant that was
    measured and rejected: it loses up to 1.2% of the population-weighted annual total across
    the 1-3% of cell-months that are all-dry. Keeping them is what makes this rule
    total-preserving, and it is the whole difference between the two.
    """
    if dry_day_rule == "none":
        return anomaly
    if dry_day_rule != "preserve":
        msg = f"Unknown dry-day rule: {dry_day_rule}"
        raise ValueError(msg)

    wet = target > cdc.DRY_DAY_THRESHOLD_MM
    kept = anomaly.where(wet, 0.0)

    month_total = anomaly.groupby("date.month").sum("date")
    kept_total = kept.groupby("date.month").sum("date")
    has_wet_day = kept_total > 0
    # 1.0 on an all-dry cell-month, so the restore below survives the rescale untouched.
    rescale = (month_total / kept_total.where(has_wet_day)).fillna(1.0)

    has_wet_day_daily = has_wet_day.sel(month=anomaly["date"].dt.month).drop_vars(
        "month"
    )
    kept = kept.where(has_wet_day_daily, anomaly)

    rescaled = kept.groupby("date.month") * rescale
    return rescaled.drop_vars("month")

apply_value_bounds(ds: xr.Dataset, bounds: tuple[float, float]) -> xr.Dataset

Clip every data variable to bounds (inclusive); NaN stays NaN.

Source code in src/climate_data/generate/scenario_daily.py

def apply_value_bounds(ds: xr.Dataset, bounds: tuple[float, float]) -> xr.Dataset:
    """Clip every data variable to ``bounds`` (inclusive); NaN stays NaN."""
    lo, hi = bounds
    if not lo < hi:
        msg = f"value bounds must satisfy lo < hi, got {bounds!r}."
        raise ValueError(msg)
    return ds.clip(lo, hi)

check_bounded_variable_scheme(target_variable: str, anomaly_scheme: str) -> None

Refuse a stabilised or yearly anomaly scheme for a value-bounded variable.

A bounded variable (relative humidity) has its raw inputs clipped to cdc.VALUE_BOUNDS before the ratio is formed, so the floor already does the job of the +1 / tapered eps in the eps-bearing schemes -- applying both would stabilise twice and shift the ratio away from the construction that was evaluated. The yearly family drops the monthly anchor that evaluation was made on. Called from the worker and, via variables_for_anomaly_scheme, from the launcher, so --target-variable all under the default scheme skips the variable with a message instead of queueing doomed jobs.

Source code in src/climate_data/generate/scenario_daily.py

def check_bounded_variable_scheme(target_variable: str, anomaly_scheme: str) -> None:
    """Refuse a stabilised or yearly anomaly scheme for a value-bounded variable.

    A bounded variable (relative humidity) has its raw inputs clipped to
    ``cdc.VALUE_BOUNDS`` before the ratio is formed, so the floor already does the job of
    the +1 / tapered eps in the eps-bearing schemes -- applying both would stabilise twice
    and shift the ratio away from the construction that was evaluated. The yearly family
    drops the monthly anchor that evaluation was made on. Called from the worker and, via
    ``variables_for_anomaly_scheme``, from the launcher, so ``--target-variable all`` under
    the default scheme skips the variable with a message instead of queueing doomed jobs.
    """
    bounds = cdc.VALUE_BOUNDS.get(target_variable)
    if bounds is None or anomaly_scheme in cdc.BOUNDED_VARIABLE_SCHEMES:
        return
    msg = (
        f"{target_variable!r} is value-bounded (inputs clipped to {bounds['input']}, output "
        f"to {bounds['output']}) and runs under {list(cdc.BOUNDED_VARIABLE_SCHEMES)} only; "
        f"got anomaly_scheme={anomaly_scheme!r}. Pass --anomaly-scheme "
        f"{cdc.BOUNDED_VARIABLE_SCHEMES[0]} for this variable."
    )
    raise ValueError(msg)

check_debias_variable(target_variable: str, debias_method: str, dry_day_rule: str = 'none') -> None

Refuse a correction for a variable it has not been validated against.

Called from the launchers as well as from the worker, so that --target-variable all fails in a second rather than after submitting thousands of doomed jobs.

Source code in src/climate_data/generate/scenario_daily.py

def check_debias_variable(
    target_variable: str, debias_method: str, dry_day_rule: str = "none"
) -> None:
    """Refuse a correction for a variable it has not been validated against.

    Called from the launchers as well as from the worker, so that ``--target-variable all``
    fails in a second rather than after submitting thousands of doomed jobs.
    """
    if debias_method != "none" and target_variable not in cdc.DEBIAS_VARIABLES:
        msg = (
            f"debias_method={debias_method!r} is not validated for {target_variable!r}. "
            f"Allowed: {list(cdc.DEBIAS_VARIABLES)}. Name the variable explicitly rather "
            "than using 'all'."
        )
        raise ValueError(msg)
    if dry_day_rule != "none" and target_variable not in cdc.DRY_DAY_VARIABLES:
        msg = (
            f"dry_day_rule={dry_day_rule!r} is not validated for {target_variable!r}. "
            f"Allowed: {list(cdc.DRY_DAY_VARIABLES)}. Name the variable explicitly rather "
            "than using 'all'."
        )
        raise ValueError(msg)

check_scheme_compatibility(anomaly_scheme: str, anomaly_type: str, debias_method: str, dry_day_rule: str) -> None

Reject combinations of the two correction axes that cannot both apply.

debias_method and dry_day_rule correct the (T + eps) / (R + eps) construction. Both monthly (constant eps) and monthly-taper (tapered eps) use it, so both accept them. The yearly and monthly-ratio families have no eps for those corrections to act on, so asking for either against them is a mistake rather than a no-op -- and silently ignoring the request would produce a file whose attrs claim a correction that was never applied.

Source code in src/climate_data/generate/scenario_daily.py

def check_scheme_compatibility(
    anomaly_scheme: str,
    anomaly_type: str,
    debias_method: str,
    dry_day_rule: str,
) -> None:
    """Reject combinations of the two correction axes that cannot both apply.

    `debias_method` and `dry_day_rule` correct the `(T + eps) / (R + eps)` construction.
    Both `monthly` (constant eps) and `monthly-taper` (tapered eps) use it, so both accept
    them. The yearly and monthly-ratio families have no eps for those corrections to act
    on, so asking for either against them is a mistake rather than a no-op -- and silently
    ignoring the request would produce a file whose attrs claim a correction that was
    never applied.
    """
    has_eps = anomaly_scheme in cdc.EPS_BEARING_SCHEMES
    if not has_eps and (debias_method != "none" or dry_day_rule != "none"):
        msg = (
            f"debias_method={debias_method!r} and dry_day_rule={dry_day_rule!r} cannot "
            f"be combined with anomaly_scheme={anomaly_scheme!r}: they correct the eps "
            f"stabiliser, which this scheme does not use. Schemes that carry an eps: "
            f"{list(cdc.EPS_BEARING_SCHEMES)}."
        )
        raise ValueError(msg)
    if anomaly_type != "multiplicative":
        msg = f"Anomaly scheme '{anomaly_scheme}' only applies to multiplicative variables."
        raise ValueError(msg)

compute_anomaly(reference: xr.Dataset, target: xr.Dataset, anomaly_type: str, *, debias_method: str, dry_day_rule: str, anomaly_scheme: str = cdc.ANOMALY_SCHEME_MONTHLY, eps_floor: float = cdc.DEFAULT_EPS_FLOOR, anomaly_cap: float | None = cdc.DEFAULT_ANOMALY_CAP) -> xr.Dataset

The forecast anomaly, optionally capped.

A thin wrapper so the ceiling applies to every scheme without threading it through each one: the scheme dispatch below has several return points, and duplicating the clip at each is how one of them ends up missing it.

Applied on the GCM grid, BEFORE interpolate_to_target_latlon. Capping afterwards would leave a blown-up cell already smeared across its neighbours by the regrid.

The cap bounds each cell's multiplier at the granularity the scheme anchors at -- annual for the yearly family, per calendar month for the monthly family -- and rescales the daily series to meet it.

It cannot be an elementwise clip. anomaly carries the target's daily date dimension, because the denominator is a reference-window mean with no date dim and the division broadcasts. Clipping elementwise bounds every day at anomaly_cap times the cell's reference mean, and against an annual denominator that ceiling sits inside the ordinary distribution of daily rainfall: a cell averaging 2 mm/day is cut at 40 mm/day, an unremarkable tropical wet day. The daily-over-annual ratio also carries the seasonal cycle, so such a clip bites hardest where seasonality is strongest rather than where the pathology is.

Measured on yearly-delta ssp126 before this change, an elementwise clip at 20 altered 20.5% of land pixels in 2083 and 15.5% in 2050; 82% of those had an annual anomaly at or below 2 -- cells projecting essentially no change -- and the global total moved -19.1% and -1.1% respectively.

Rescaling instead leaves any cell at or below the ceiling bit-identical to uncapped, brings a cell above it exactly to the ceiling, and preserves within-period shape so a wet day stays proportionally a wet day.

The cap still only removes precipitation -- it cannot conserve, and whatever it removes is gone from the total. That is intended: the values it removes are ones no cell's observed climatology supports. But it means a cap moves the level, so it must be chosen on evidence rather than set defensively.

Source code in src/climate_data/generate/scenario_daily.py

def compute_anomaly(
    reference: xr.Dataset,
    target: xr.Dataset,
    anomaly_type: str,
    *,
    debias_method: str,
    dry_day_rule: str,
    anomaly_scheme: str = cdc.ANOMALY_SCHEME_MONTHLY,
    eps_floor: float = cdc.DEFAULT_EPS_FLOOR,
    anomaly_cap: float | None = cdc.DEFAULT_ANOMALY_CAP,
) -> xr.Dataset:
    """The forecast anomaly, optionally capped.

    A thin wrapper so the ceiling applies to every scheme without threading it through
    each one: the scheme dispatch below has several return points, and duplicating the
    clip at each is how one of them ends up missing it.

    Applied on the GCM grid, BEFORE `interpolate_to_target_latlon`. Capping afterwards
    would leave a blown-up cell already smeared across its neighbours by the regrid.

    The cap bounds each cell's multiplier at the granularity the scheme anchors at --
    annual for the yearly family, per calendar month for the monthly family -- and
    rescales the daily series to meet it.

    It cannot be an elementwise clip. `anomaly` carries the target's daily date
    dimension, because the denominator is a reference-window mean with no date dim and
    the division broadcasts. Clipping elementwise bounds every *day* at `anomaly_cap`
    times the cell's reference mean, and against an annual denominator that ceiling sits
    inside the ordinary distribution of daily rainfall: a cell averaging 2 mm/day is cut
    at 40 mm/day, an unremarkable tropical wet day. The daily-over-annual ratio also
    carries the seasonal cycle, so such a clip bites hardest where seasonality is
    strongest rather than where the pathology is.

    Measured on `yearly-delta` ssp126 before this change, an elementwise clip at 20
    altered 20.5% of land pixels in 2083 and 15.5% in 2050; 82% of those had an annual
    anomaly at or below 2 -- cells projecting essentially no change -- and the global
    total moved -19.1% and -1.1% respectively.

    Rescaling instead leaves any cell at or below the ceiling bit-identical to uncapped,
    brings a cell above it exactly to the ceiling, and preserves within-period shape so
    a wet day stays proportionally a wet day.

    The cap still only removes precipitation -- it cannot conserve, and whatever it
    removes is gone from the total. That is intended: the values it removes are ones no
    cell's observed climatology supports. But it means a cap moves the level, so it must
    be chosen on evidence rather than set defensively.
    """
    anomaly = _compute_anomaly_uncapped(
        reference,
        target,
        anomaly_type,
        debias_method=debias_method,
        dry_day_rule=dry_day_rule,
        anomaly_scheme=anomaly_scheme,
        eps_floor=eps_floor,
    )
    if anomaly_cap is None:
        return anomaly
    if anomaly_type != "multiplicative":
        msg = (
            f"anomaly_cap={anomaly_cap!r} was requested for an additive anomaly. The cap "
            "bounds a multiplier; an additive anomaly is a difference in the variable's "
            "own units, where a ceiling of this kind is meaningless."
        )
        raise ValueError(msg)
    if anomaly_cap <= 0:
        msg = f"anomaly_cap must be positive, got {anomaly_cap!r}."
        raise ValueError(msg)
    if anomaly_scheme in cdc.YEARLY_ANOMALY_SCHEMES:
        # anchored to the annual level, so bound the annual multiplier
        over = anomaly.mean("date")
        factor = (anomaly_cap / over).where(over > anomaly_cap, 1.0)
        return anomaly * factor

    # the monthly family anchors each calendar month to ERA5 separately, so the ceiling
    # belongs per month. Bounding the annual mean instead would let one blown-up month
    # drag the rescale factor down and crush the other eleven, which are fine.
    over = anomaly.groupby("date.month").mean("date")
    factor = (anomaly_cap / over).where(over > anomaly_cap, 1.0)
    return anomaly.groupby("date.month") * factor

jensen_debias_factor(reference: xr.Dataset, reference_monthly: xr.Dataset, debias_method: str, eps: xr.Dataset | float = 1.0) -> xr.Dataset

The factor to divide a multiplicative anomaly by, per month and per GCM cell.

The anomaly is (T + 1) / (R + 1) with R a monthly mean over only five reference years. 1/(R + 1) is convex, so by Jensen's inequality the anomaly averages above 1 even when the target year is drawn from the same distribution as the reference period -- a level bias on every forecast year. This returns an estimate of that inflation.

loo -- leave-one-out. For each held-out reference year, form the multiplier of that year against the mean of the other years and average over the folds. Between reference years there is no climate signal, so an unbiased estimator would return 1; the excess is the bias, measured from the data with no series expansion. The held-out denominator averages n-1 years while the pipeline averages n, and the bias goes as 1/n, so the excess is rescaled by (n-1)/n.

This is provably >= 1: with u_y = T_y + 1 and S = sum_y u_y the held-out denominator is (S - u_y)/(n-1), so each fold is (n-1)*u_y/(S - u_y), convex in u_y, and Jensen gives mean_y >= f(S/n) = 1 with equality iff every reference year is identical. So dividing by it can only shrink the anomaly, never inflate it -- which is what makes the effect on a threshold count such as precipitation_days sign-definite.

analytic -- the second-order expansion 1 + Var(Rbar)/(R + 1)^2. Cheaper to reason about but it is a truncated series, and the neglected terms matter exactly where the correction is largest (near-zero R, where eps dominates the denominator). Kept for comparison; loo is the estimator this was built for.

reference_monthly is the pipeline's own denominator, passed in rather than recomputed so the analytic form squares precisely the value the anomaly divides by.

Source code in src/climate_data/generate/scenario_daily.py

def jensen_debias_factor(
    reference: xr.Dataset,
    reference_monthly: xr.Dataset,
    debias_method: str,
    eps: xr.Dataset | float = 1.0,
) -> xr.Dataset:
    """The factor to divide a multiplicative anomaly by, per month and per GCM cell.

    The anomaly is ``(T + 1) / (R + 1)`` with ``R`` a monthly mean over only five reference
    years. ``1/(R + 1)`` is convex, so by Jensen's inequality the anomaly averages above 1 even
    when the target year is drawn from the same distribution as the reference period -- a level
    bias on every forecast year. This returns an estimate of that inflation.

    ``loo`` -- leave-one-out. For each held-out reference year, form the multiplier of that
    year against the mean of the *other* years and average over the folds. Between reference
    years there is no climate signal, so an unbiased estimator would return 1; the excess is
    the bias, measured from the data with no series expansion. The held-out denominator
    averages ``n-1`` years while the pipeline averages ``n``, and the bias goes as ``1/n``, so
    the excess is rescaled by ``(n-1)/n``.

    This is provably ``>= 1``: with ``u_y = T_y + 1`` and ``S = sum_y u_y`` the held-out
    denominator is ``(S - u_y)/(n-1)``, so each fold is ``(n-1)*u_y/(S - u_y)``, convex in
    ``u_y``, and Jensen gives ``mean_y >= f(S/n) = 1`` with equality iff every reference year
    is identical. So dividing by it can only shrink the anomaly, never inflate it -- which is
    what makes the effect on a threshold count such as ``precipitation_days`` sign-definite.

    ``analytic`` -- the second-order expansion ``1 + Var(Rbar)/(R + 1)^2``. Cheaper to reason
    about but it is a truncated series, and the neglected terms matter exactly where the
    correction is largest (near-zero ``R``, where ``eps`` dominates the denominator). Kept for
    comparison; ``loo`` is the estimator this was built for.

    ``reference_monthly`` is the pipeline's own denominator, passed in rather than recomputed so
    the analytic form squares precisely the value the anomaly divides by.
    """
    by_year = _monthly_means_by_reference_year(reference)
    n_years = by_year.sizes["reference_year"]

    if debias_method == "loo":
        mean_year = by_year.mean("reference_year")
        folds = []
        valid = []
        for i in range(n_years):
            held_out = by_year.isel(reference_year=i)
            others = (n_years * mean_year - held_out) / (n_years - 1)
            fold_denominator = others + eps
            folds.append(
                (held_out + eps) / fold_denominator.where(fold_denominator > 0)
            )
            valid.append((fold_denominator > 0).astype("int8"))
        raw = xr.concat(folds, dim="reference_year").mean("reference_year")
        factor = 1.0 + ((n_years - 1) / n_years) * (raw - 1.0)
        # Where any fold is undefined the estimate is undefined, so apply no correction.
        # Averaging the surviving folds is NOT a repair: the undefined one is precisely
        # the large fold that lifts the mean above 1, so dropping it collapses the factor
        # (0.2 for a one-wet-four-dry cell-month) and the de-bias INFLATES the anomaly --
        # the opposite of its purpose. Only reachable when eps can be zero, i.e. under the
        # taper; the constant eps = 1 floors every denominator at 1.
        factor = factor.where(sum(valid) == n_years, 1.0).fillna(1.0)
    elif debias_method == "analytic":
        variance = by_year.var("reference_year", ddof=1)
        factor = 1.0 + (variance / n_years) / (reference_monthly + eps) ** 2
    else:
        msg = f"Unknown debias method: {debias_method}"
        raise ValueError(msg)

    factor = factor.drop_vars("reference_year", errors="ignore")
    values = factor.to_dataarray()
    if not bool(np.isfinite(values).all()):
        msg = (
            f"Non-finite value in the {debias_method} de-bias factor. Interpolation would "
            "silently fill it from a neighbour rather than surface it, so refusing to proceed."
        )
        raise ValueError(msg)
    minimum = float(values.min())
    if minimum < 1.0 - _DEBIAS_FACTOR_TOLERANCE:
        # Both estimators are >= 1 by construction -- loo by Jensen on a convex fold,
        # analytic because it adds a non-negative term. A value below 1 means the
        # construction's precondition failed, and dividing by it would inflate rather
        # than shrink the anomaly. Refuse rather than ship an anti-correction.
        msg = (
            f"The {debias_method} de-bias factor reached {minimum!r}, below its "
            "guaranteed lower bound of 1. Dividing by it would inflate the anomaly."
        )
        raise ValueError(msg)
    return factor

load_and_shift_longitude_and_correct_time(member_path: str | Path, year: str) -> xr.Dataset

Put a member's year onto the real Gregorian calendar, day by day.

The conversion is by DATE, never by value. interp_calendar used to resample onto the target axis by linear interpolation, so whenever the source calendar's year length differed from the target's -- a noleap member in a leap year -- every target day became a blend of two source days. convert_calendar maps each source day onto its own date instead and leaves 29 February missing; the reindex holds the output axis to exactly this year's days whatever the source calendar, and interpolate_na fills the gap from the nearest real day.

Every variable comes through here, not just precipitation, but what the blending cost depends on the field. Temperature is smooth, so blending barely moves it and the threshold measures it feeds -- days_over_30C, the suitability maps -- were knocked either way and largely cancelled in the spatial mean. Precipitation is spiky and mostly zero against a 0.1 mm cut sitting on the floor, so smearing a wet day onto its dry neighbours could only push them UP over the line, never below it. That one-way ratchet is why precipitation_days took a coherent ~13.5 d per noleap member in all 19 leap years 2024-2096 while everything else came out as cancelling noise. (CLIMATE-35)

No align_on is passed: it only takes effect when a 360_day calendar is involved, and no member in the extracts uses one. If that changes, align_on needs a deliberate choice rather than xarray's default.

Source code in src/climate_data/generate/scenario_daily.py

def load_and_shift_longitude_and_correct_time(
    member_path: str | Path,
    year: str,
) -> xr.Dataset:
    """Put a member's year onto the real Gregorian calendar, day by day.

    The conversion is by DATE, never by value. `interp_calendar` used to resample onto the
    target axis by linear interpolation, so whenever the source calendar's year length
    differed from the target's -- a `noleap` member in a leap year -- every target day
    became a blend of two source days. `convert_calendar` maps each source day onto its own
    date instead and leaves 29 February missing; the `reindex` holds the output axis to
    exactly this year's days whatever the source calendar, and `interpolate_na` fills the
    gap from the nearest real day.

    Every variable comes through here, not just precipitation, but what the blending cost
    depends on the field. Temperature is smooth, so blending barely moves it and the
    threshold measures it feeds -- `days_over_30C`, the suitability maps -- were knocked
    either way and largely cancelled in the spatial mean. Precipitation is spiky and mostly
    zero against a 0.1 mm cut sitting on the floor, so smearing a wet day onto its dry
    neighbours could only push them UP over the line, never below it. That one-way ratchet
    is why `precipitation_days` took a coherent ~13.5 d per noleap member in all 19 leap
    years 2024-2096 while everything else came out as cancelling noise. (CLIMATE-35)

    No `align_on` is passed: it only takes effect when a `360_day` calendar is involved,
    and no member in the extracts uses one. If that changes, `align_on` needs a deliberate
    choice rather than xarray's default.
    """
    time_slice = slice(f"{year}-01-01", f"{year}-12-31")
    time_range = pd.date_range(f"{year}-01-01", f"{year}-12-31")
    ds = load_and_shift_longitude(member_path, time_slice)
    ds = (
        ds.assign_coords(time=ds.time.dt.floor("D"))
        .convert_calendar("standard")
        .reindex(time=time_range)
        .interpolate_na(dim="time", method="nearest", fill_value="extrapolate")
        .rename({"time": "date"})
    )
    return ds

variables_for_anomaly_scheme(target_variables: list[str], anomaly_scheme: str, anomaly_types: dict[str, str]) -> list[str]

Drop target variables the anomaly scheme cannot be applied to.

The yearly schemes are defined for multiplicative variables only -- compute_anomaly raises for anything else -- while --target-variable defaults to ALL. Without this filter the default invocation of a yearly scheme submits additive jobs that are certain to fail, and they fail only after being scheduled and retried.

Skipped variables are named rather than dropped quietly, and selecting nothing but additive variables is an error rather than an empty run.

Source code in src/climate_data/generate/scenario_daily.py

def variables_for_anomaly_scheme(
    target_variables: list[str],
    anomaly_scheme: str,
    anomaly_types: dict[str, str],
) -> list[str]:
    """Drop target variables the anomaly scheme cannot be applied to.

    The yearly schemes are defined for multiplicative variables only -- `compute_anomaly`
    raises for anything else -- while `--target-variable` defaults to ALL. Without this
    filter the default invocation of a yearly scheme submits additive jobs that are
    certain to fail, and they fail only after being scheduled and retried.

    Skipped variables are named rather than dropped quietly, and selecting nothing but
    additive variables is an error rather than an empty run.
    """
    keep = []
    skip = []
    skip_bounded = []
    for variable in target_variables:
        if (
            variable in cdc.VALUE_BOUNDS
            and anomaly_scheme not in cdc.BOUNDED_VARIABLE_SCHEMES
        ):
            # Value-bounded variables are floored explicitly; see
            # check_bounded_variable_scheme for why the other schemes are refused.
            skip_bounded.append(variable)
        elif (
            anomaly_scheme != cdc.ANOMALY_SCHEME_MONTHLY
            and anomaly_types[variable] != "multiplicative"
        ):
            skip.append(variable)
        else:
            keep.append(variable)

    if skip:
        print(
            f"Anomaly scheme '{anomaly_scheme}' applies to multiplicative variables only;"
            f" skipping {len(skip)}: {', '.join(skip)}."
        )
    if skip_bounded:
        print(
            f"Anomaly scheme '{anomaly_scheme}' does not apply to value-bounded variables,"
            f" which run under {', '.join(cdc.BOUNDED_VARIABLE_SCHEMES)} only;"
            f" skipping {len(skip_bounded)}: {', '.join(skip_bounded)}."
        )
    if not keep:
        msg = (
            f"No multiplicative variables selected, so anomaly scheme '{anomaly_scheme}'"
            f" has nothing to run. Selected: {', '.join(target_variables)}."
        )
        raise click.UsageError(msg)
    return keep

utils

annual_mean_from_monthly(ds: xr.Dataset) -> xr.Dataset

Day-weighted annual mean of a 12-month climatology.

Weights are bound by month label, not position, so a permuted month coordinate still weights correctly; anything other than a full 1..12 coordinate is rejected rather than misweighted.

Source code in src/climate_data/generate/utils.py

def annual_mean_from_monthly(ds: xr.Dataset) -> xr.Dataset:
    """Day-weighted annual mean of a 12-month climatology.

    Weights are bound by month label, not position, so a permuted month
    coordinate still weights correctly; anything other than a full 1..12
    coordinate is rejected rather than misweighted.
    """
    months = np.arange(1, 13)
    if not np.array_equal(np.sort(ds["month"].to_numpy()), months):
        msg = f"expected a full 1..12 month coordinate, got {ds['month'].to_numpy()}"
        raise ValueError(msg)
    weights = xr.DataArray(cdc.DAYS_IN_MONTH, dims=["month"], coords={"month": months})
    return ds.weighted(weights).mean("month")

buck_vapor_pressure(temperature_c: xr.Dataset) -> xr.Dataset

Approximate vapor pressure of water.

https://en.wikipedia.org/wiki/Arden_Buck_equation https://journals.ametsoc.org/view/journals/apme/20/12/1520-0450_1981_020_1527_nefcvp_2_0_co_2.xml

Parameters

temperature_c Temperature in Celsius

Returns

xr.Dataset Vapor pressure in hPa

Source code in src/climate_data/generate/utils.py

def buck_vapor_pressure(temperature_c: xr.Dataset) -> xr.Dataset:
    """Approximate vapor pressure of water.

    https://en.wikipedia.org/wiki/Arden_Buck_equation
    https://journals.ametsoc.org/view/journals/apme/20/12/1520-0450_1981_020_1527_nefcvp_2_0_co_2.xml

    Parameters
    ----------
    temperature_c
        Temperature in Celsius

    Returns
    -------
    xr.Dataset
        Vapor pressure in hPa
    """
    over_water = 6.1121 * np.exp(
        (18.678 - temperature_c / 234.5) * (temperature_c / (257.14 + temperature_c))
    )
    over_ice = 6.1115 * np.exp(
        (23.036 - temperature_c / 333.7) * (temperature_c / (279.82 + temperature_c))
    )
    vp = xr.where(temperature_c > 0, over_water, over_ice)  # type: ignore[no-untyped-call]
    return vp  # type: ignore[no-any-return]

daily_accumulation_last(ds: xr.Dataset) -> xr.Dataset

Collapse a within-day cumulative variable to its daily total.

For ERA5-Land total_precipitation, which accumulates since 00Z, the day's total is the window's closing sample, so last reads it directly. Prefer this to a maximum: the two agree only while the window rises monotonically, and int16 packing can make it tick down in its final step, in which case the maximum is an earlier, quantisation-inflated sample. Measured over 741,270 land day-pixels, last matched the true close for 100.0000% against 99.9864% for max.

Note resample emits a regular axis with NaN for empty bins, unlike groupby, which emits observed groups only. Callers must trim the partial bins at each end.

skipna=False because the default steps back past an absent closing sample to hour 23, reintroducing the incomplete window this function exists to eliminate. NaN is the honest answer for a window that never closed. It is not, however, a detectable one: generate_historical_daily_main fills ERA5-Land NaNs from the interpolated single-level field -- that is how ocean pixels are supplied -- so such a pixel carries a complete 0.25 degree value rather than a truncated 0.1 degree one, and never reaches validate_output. Preferring the coarse-but-whole value is the point; detecting the substitution would need a separate check against the sea mask.

Source code in src/climate_data/generate/utils.py

def daily_accumulation_last(ds: xr.Dataset) -> xr.Dataset:
    """Collapse a within-day cumulative variable to its daily total.

    For ERA5-Land `total_precipitation`, which accumulates since 00Z, the day's total is
    the window's closing sample, so `last` reads it directly. Prefer this to a maximum:
    the two agree only while the window rises monotonically, and int16 packing can make
    it tick down in its final step, in which case the maximum is an earlier,
    quantisation-inflated sample. Measured over 741,270 land day-pixels, `last` matched
    the true close for 100.0000% against 99.9864% for `max`.

    Note `resample` emits a regular axis with NaN for empty bins, unlike `groupby`, which
    emits observed groups only. Callers must trim the partial bins at each end.

    `skipna=False` because the default steps back past an absent closing sample to hour
    23, reintroducing the incomplete window this function exists to eliminate. NaN is the
    honest answer for a window that never closed. It is not, however, a detectable one:
    `generate_historical_daily_main` fills ERA5-Land NaNs from the interpolated
    single-level field -- that is how ocean pixels are supplied -- so such a pixel carries
    a complete 0.25 degree value rather than a truncated 0.1 degree one, and never reaches
    `validate_output`. Preferring the coarse-but-whole value is the point; detecting the
    substitution would need a separate check against the sea mask.
    """
    resampled = ds.resample(
        time=ACCUMULATION_FREQ,
        closed=ACCUMULATION_CLOSED,
        label=ACCUMULATION_LABEL,
    ).last(skipna=False)
    return resampled.rename({"time": "date"})

daily_accumulation_sum(ds: xr.Dataset) -> xr.Dataset

Collapse a variable of per-hour increments to its daily total.

For ERA5 single-levels, "accumulations are over the hour ending at the validity date/time", so the increments sum. The same interval convention applies as for daily_accumulation_last, so the two stay on a common date axis -- they are merged downstream in generate_historical_daily_main.

Source code in src/climate_data/generate/utils.py

def daily_accumulation_sum(ds: xr.Dataset) -> xr.Dataset:
    """Collapse a variable of per-hour increments to its daily total.

    For ERA5 single-levels, "accumulations are over the hour ending at the validity
    date/time", so the increments sum. The same interval convention applies as for
    `daily_accumulation_last`, so the two stay on a common `date` axis -- they are merged
    downstream in `generate_historical_daily_main`.
    """
    resampled = ds.resample(
        time=ACCUMULATION_FREQ,
        closed=ACCUMULATION_CLOSED,
        label=ACCUMULATION_LABEL,
    ).sum()
    return resampled.rename({"time": "date"})

identity(ds: xr.Dataset) -> xr.Dataset

Identity transformation

Source code in src/climate_data/generate/utils.py

def identity(ds: xr.Dataset) -> xr.Dataset:
    """Identity transformation"""
    return ds

interpolate_to_target_latlon(ds: xr.Dataset, method: str = 'nearest', target_lon: xr.DataArray = cdc.TARGET_LONGITUDE, target_lat: xr.DataArray = cdc.TARGET_LATITUDE) -> xr.Dataset

Interpolate a dataset to a target latitude and longitude grid.

Parameters

ds Dataset to interpolate method Interpolation method target_lon Target longitude grid target_lat Target latitude grid

Returns

xr.Dataset Interpolated dataset

Source code in src/climate_data/generate/utils.py

def interpolate_to_target_latlon(
    ds: xr.Dataset,
    method: str = "nearest",
    target_lon: xr.DataArray = cdc.TARGET_LONGITUDE,
    target_lat: xr.DataArray = cdc.TARGET_LATITUDE,
) -> xr.Dataset:
    """Interpolate a dataset to a target latitude and longitude grid.

    Parameters
    ----------
    ds
        Dataset to interpolate
    method
        Interpolation method
    target_lon
        Target longitude grid
    target_lat
        Target latitude grid

    Returns
    -------
    xr.Dataset
        Interpolated dataset
    """
    return (
        ds.interp(longitude=target_lon, latitude=target_lat, method=method)  # type: ignore[arg-type]
        .interpolate_na(dim="longitude", method="nearest", fill_value="extrapolate")
        .sortby("latitude")
        .interpolate_na(dim="latitude", method="nearest", fill_value="extrapolate")
        .sortby("latitude", ascending=False)
    )

kelvin_to_celsius(temperature_k: xr.Dataset) -> xr.Dataset

Convert temperature from Kelvin to Celsius

Parameters

temperature_k Temperature in Kelvin

Returns

xr.Dataset Temperature in Celsius

Source code in src/climate_data/generate/utils.py

def kelvin_to_celsius(temperature_k: xr.Dataset) -> xr.Dataset:
    """Convert temperature from Kelvin to Celsius

    Parameters
    ----------
    temperature_k
        Temperature in Kelvin

    Returns
    -------
    xr.Dataset
        Temperature in Celsius
    """
    return temperature_k - 273.15

meter_to_millimeter(rainfall_m: xr.Dataset) -> xr.Dataset

Convert rainfall from meters to millimeters

Parameters

rainfall_m Rainfall in meters

Returns

xr.Dataset Rainfall in millimeters

Source code in src/climate_data/generate/utils.py

def meter_to_millimeter(rainfall_m: xr.Dataset) -> xr.Dataset:
    """Convert rainfall from meters to millimeters

    Parameters
    ----------
    rainfall_m
        Rainfall in meters

    Returns
    -------
    xr.Dataset
        Rainfall in millimeters
    """
    return 1000 * rainfall_m

parse_reference_years(reference_years: str) -> slice

Turn a START-END year string (inclusive, four-digit years) into a date slice.

Source code in src/climate_data/generate/utils.py

def parse_reference_years(reference_years: str) -> slice:
    """Turn a START-END year string (inclusive, four-digit years) into a date slice."""
    start, sep, end = reference_years.partition("-")
    valid = (
        bool(sep)
        and len(start) == 4  # noqa: PLR2004
        and len(end) == 4  # noqa: PLR2004
        and start.isdigit()
        and end.isdigit()
        and int(start) <= int(end)
    )
    if not valid:
        msg = (
            "reference-years must be 'START-END' with four-digit years and "
            f"START <= END, got {reference_years!r}"
        )
        raise ValueError(msg)
    return slice(f"{start}-01-01", f"{end}-12-31")

precipitation_flux_to_rainfall(precipitation_flux: xr.Dataset) -> xr.Dataset

Convert precipitation flux to rainfall

Parameters

precipitation_flux Precipitation flux in kg m-2 s-1

Returns

xr.Dataset Rainfall in mm/day

Source code in src/climate_data/generate/utils.py

def precipitation_flux_to_rainfall(precipitation_flux: xr.Dataset) -> xr.Dataset:
    """Convert precipitation flux to rainfall

    Parameters
    ----------
    precipitation_flux
        Precipitation flux in kg m-2 s-1

    Returns
    -------
    xr.Dataset
        Rainfall in mm/day
    """
    seconds_per_day = 86400
    mm_per_kg_m2 = 1
    return seconds_per_day * mm_per_kg_m2 * precipitation_flux

rh_percent(temperature_c: xr.Dataset, dewpoint_temperature_c: xr.Dataset) -> xr.Dataset

Calculate relative humidity from temperature and dewpoint temperature.

Parameters

temperature_c Temperature in Celsius dewpoint_temperature_c Dewpoint temperature in Celsius

Returns

xr.Dataset Relative humidity as a percentage

Source code in src/climate_data/generate/utils.py

def rh_percent(
    temperature_c: xr.Dataset, dewpoint_temperature_c: xr.Dataset
) -> xr.Dataset:
    """Calculate relative humidity from temperature and dewpoint temperature.

    Parameters
    ----------
    temperature_c
        Temperature in Celsius
    dewpoint_temperature_c
        Dewpoint temperature in Celsius

    Returns
    -------
    xr.Dataset
        Relative humidity as a percentage
    """
    # saturation vapour pressure
    svp = buck_vapor_pressure(temperature_c)
    # actual vapour pressure
    vp = buck_vapor_pressure(dewpoint_temperature_c)
    return 100 * vp / svp

scale_wind_speed_height(wind_speed_10m: xr.Dataset) -> xr.Dataset

Scaling wind speed from a height of 10 meters to a height of 2 meters

Reference: Bröde et al. (2012) https://doi.org/10.1007/s00484-011-0454-1

Parameters

wind_speed_10m The 10m wind speed [m/s]. May be signed (ie a velocity component)

Returns

xr.DataSet The 2m wind speed [m/s]. May be signed (ie a velocity component)

Source code in src/climate_data/generate/utils.py

def scale_wind_speed_height(wind_speed_10m: xr.Dataset) -> xr.Dataset:
    """Scaling wind speed from a height of 10 meters to a height of 2 meters

    Reference: Bröde et al. (2012)
    https://doi.org/10.1007/s00484-011-0454-1

    Parameters
    ----------
    wind_speed_10m
        The 10m wind speed [m/s]. May be signed (ie a velocity component)

    Returns
    -------
    xr.DataSet
        The 2m wind speed [m/s]. May be signed (ie a velocity component)
    """
    scale_factor = np.log10(2 / 0.01) / np.log10(10 / 0.01)
    return scale_factor * wind_speed_10m  # type: ignore[no-any-return]

vector_magnitude(x: xr.Dataset, y: xr.Dataset) -> xr.Dataset

Calculate the magnitude of a vector.

Source code in src/climate_data/generate/utils.py

def vector_magnitude(x: xr.Dataset, y: xr.Dataset) -> xr.Dataset:
    """Calculate the magnitude of a vector."""
    return np.sqrt(x**2 + y**2)  # type: ignore[return-value]