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
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
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
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
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
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
572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 | |
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
259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 | |
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
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.