Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions changelog/887.fix.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,5 @@
Fabricated CMIP7 test data now reaches years without real CMIP6 source data by repeating the final year,
rather than relabelling the whole series onto the requested end date.
The real years keep their real dates, so a diagnostic can ask for a period that the source data genuinely covers.
This also fixes datasets split across several files, where each file was dragged onto the same end date
and the converted files overlapped.
92 changes: 55 additions & 37 deletions packages/climate-ref-core/src/climate_ref_core/cmip6_to_cmip7.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,17 +22,15 @@
from dataclasses import dataclass
from datetime import UTC, datetime, timedelta
from importlib import resources
from typing import TYPE_CHECKING, Any
from typing import Any

import attrs
import cftime
import numpy as np
import pandas as pd
import xarray as xr
from loguru import logger

if TYPE_CHECKING:
import xarray as xr


def suppress_bounds_coordinates(ds: xr.Dataset) -> xr.Dataset:
"""
Expand Down Expand Up @@ -496,15 +494,17 @@ def _month_index(t: Any) -> int:
_MONTHS_PER_YEAR = 12


def shift_time_axis_end(ds: xr.Dataset, end_year: int, end_month: int = 12) -> xr.Dataset:
def repeat_final_year_to(ds: xr.Dataset, end_year: int, end_month: int = 12) -> xr.Dataset:
"""
Relabel a monthly time axis so the final timestep lands on ``end_year``-``end_month``.
Extend a monthly series to ``end_year``-``end_month`` by repeating its final year.

The real timesteps keep their real dates and values. Only the tail is fabricated, by
tiling the last twelve months forward until the series reaches the requested end. This
gives CMIP7 ``historical`` coverage for years the CMIP6 source never ran (e.g. 2015-2021)
without moving the years that did.

The whole ``time`` coordinate (and any ``time_bnds``) is shifted by a constant
whole-month offset, so monthly spacing, calendar, and cftime types are preserved
and only the labels move. This is used to fabricate CMIP7 ``historical`` coverage
that reaches years for which no real data exists (e.g. 2002-2021 for the fire
diagnostic), without altering the underlying data values.
A series shorter than a year repeats whatever it has. A series that already reaches the
requested end is returned unchanged.

Only monthly data with a ``cftime`` time axis is supported; datasets without a
``time`` coordinate (fixed-frequency, e.g. ``sftlf``) are returned unchanged.
Expand All @@ -521,7 +521,7 @@ def shift_time_axis_end(ds: xr.Dataset, end_year: int, end_month: int = 12) -> x
Returns
-------
xr.Dataset
A shallow copy with the relabelled time axis.
A dataset whose time axis runs to the requested end.
"""
if not 1 <= end_month <= _MONTHS_PER_YEAR:
raise ValueError(f"end_month must be in 1..12, got {end_month}")
Expand All @@ -533,14 +533,13 @@ def shift_time_axis_end(ds: xr.Dataset, end_year: int, end_month: int = 12) -> x
last = time_values[-1]
if not isinstance(last, cftime.datetime):
raise TypeError(
"shift_time_axis_end requires a cftime time axis; "
"repeat_final_year_to requires a cftime time axis; "
f"got {type(last).__name__}. Decode with use_cftime=True."
)

# Whole-month offset that moves the final label onto the requested end.
target_index = end_year * 12 + (end_month - 1)
offset_months = target_index - _month_index(last)
if offset_months == 0:
target_index = end_year * _MONTHS_PER_YEAR + (end_month - 1)
missing = target_index - _month_index(last)
if missing <= 0:
return ds

calendar = last.calendar
Expand All @@ -554,31 +553,50 @@ def _days_in_month(year: int, month: int) -> int:
last_of_month = first_of_next - timedelta(days=1) # type: ignore[operator]
return int(last_of_month.day) # type: ignore[attr-defined]

def _shift(t: cftime.datetime) -> cftime.datetime:
total = _month_index(t) + offset_months
year, month = divmod(total, 12)
def _add_months(t: cftime.datetime, months: int) -> cftime.datetime:
year, month = divmod(_month_index(t) + months, _MONTHS_PER_YEAR)
month += 1
# Clamp the day to the target month/calendar: a non-multiple-of-12 offset
# can land a day-31 (or leap Feb-29) label on a shorter month, which cftime
# would reject. Our monthly data is mid-month so this is normally a no-op.
# Clamp the day to the target month/calendar: a day-31 (or leap Feb-29) label
# can land on a shorter month, which cftime would reject.
day = min(t.day, _days_in_month(year, month))
return cftime.datetime(year, month, day, t.hour, t.minute, t.second, t.microsecond, calendar=calendar)

ds = ds.copy(deep=False)
shifted = np.array([_shift(t) for t in time_values])
new_time = ds["time"].copy(data=shifted)
# Tile the final year, so each fabricated month reuses the same month one or more years back.
period = min(_MONTHS_PER_YEAR, len(time_values))
positions = [len(time_values) - period + (step % period) for step in range(missing)]
offsets = [
_month_index(last) + 1 + step - _month_index(time_values[position])
for step, position in enumerate(positions)
]

padding = ds.isel(time=positions)
padded_times = np.array(
[_add_months(time_values[position], offset) for position, offset in zip(positions, offsets)]
)
new_time = padding["time"].copy(data=padded_times)
new_time.encoding = dict(ds["time"].encoding)
ds = ds.assign_coords(time=new_time)
padding = padding.assign_coords(time=new_time)

# Shift the matching time bounds, if present, so the axis stays self-consistent.
# Move the matching time bounds with their timesteps, so the axis stays self-consistent.
bounds_name = ds["time"].attrs.get("bounds")
if bounds_name and bounds_name in ds:
bnds_values = ds[bounds_name].values
shifted_bnds = np.array([[_shift(v) for v in row] for row in bnds_values])
ds[bounds_name] = ds[bounds_name].copy(data=shifted_bnds)
padded_bnds = np.array(
[
[_add_months(bound, offset) for bound in bnds_values[position]]
for position, offset in zip(positions, offsets)
]
)
padding[bounds_name] = padding[bounds_name].copy(data=padded_bnds)

logger.debug(f"Shifted time axis by {offset_months} months so it ends {end_year:04d}-{end_month:02d}")
return ds
extended = xr.concat([ds, padding], dim="time", data_vars="minimal", coords="minimal")
extended["time"].attrs = dict(ds["time"].attrs)
extended["time"].encoding = dict(ds["time"].encoding)

logger.debug(
f"Repeated the final year for {missing} months so the series ends {end_year:04d}-{end_month:02d}"
)
return extended


def convert_cmip6_dataset(
Expand All @@ -604,10 +622,10 @@ def convert_cmip6_dataset(
inplace
If True, modify the dataset in place; otherwise return a copy
extend_historical_to
Opt-in ``(end_year, end_month)``. When set, the ``time`` axis is relabelled
via :func:`shift_time_axis_end` so the series ends on that month, fabricating
CMIP7 coverage for years without real data. Defaults to ``None`` (time axis
untouched), so existing conversions are byte-identical.
Opt-in ``(end_year, end_month)``. When set, the series is padded out to that month
via :func:`repeat_final_year_to`, fabricating CMIP7 coverage for years without real
data. Defaults to ``None`` (time axis untouched), so existing conversions are
byte-identical.

Returns
-------
Expand All @@ -619,7 +637,7 @@ def convert_cmip6_dataset(

if extend_historical_to is not None:
end_year, end_month = extend_historical_to
ds = shift_time_axis_end(ds, end_year=end_year, end_month=end_month)
ds = repeat_final_year_to(ds, end_year=end_year, end_month=end_month)

# Determine the primary variable (skip coordinates/bounds)
data_vars = [str(v) for v in ds.data_vars if not str(v).endswith("_bnds") and v not in ds.coords]
Expand Down
22 changes: 11 additions & 11 deletions packages/climate-ref-core/src/climate_ref_core/esgf/cmip7.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
format_cmip7_time_range,
get_dreq_entry,
get_frequency_from_table,
shift_time_axis_end,
repeat_final_year_to,
suppress_bounds_coordinates,
)
from climate_ref_core.data import resolve_cache_dir
Expand Down Expand Up @@ -51,8 +51,8 @@ def _convert_file_to_cmip7(
cmip7_facets
CMIP7 facets for the output path
extend_historical_to
Opt-in ``(end_year, end_month)`` passed to :func:`convert_cmip6_dataset`
to relabel the time axis so historical coverage reaches that month.
Opt-in ``(end_year, end_month)``. The series is padded out to that month by
repeating its final year, so historical coverage reaches it.
Defaults to ``None`` (time axis untouched).

Returns
Expand Down Expand Up @@ -85,12 +85,12 @@ def _convert_file_to_cmip7(

time_coder = xr.coders.CFDatetimeCoder(use_cftime=True)
with xr.open_dataset(cmip6_path, decode_times=time_coder) as ds:
# When fabricating extended historical coverage, relabel the time axis first
# so both the filename time range and the written data reflect the new dates.
# When fabricating extended historical coverage, pad the series first so both the
# filename time range and the written data reflect the added months.
source_ds = ds
if extend_historical_to is not None:
end_year, end_month = extend_historical_to
source_ds = shift_time_axis_end(ds, end_year=end_year, end_month=end_month)
source_ds = repeat_final_year_to(ds, end_year=end_year, end_month=end_month)

frequency = str(cmip7_facets.get("frequency", "mon"))
time_range = format_cmip7_time_range(source_ds, frequency)
Expand Down Expand Up @@ -199,11 +199,11 @@ def __init__(
time_span
Optional time range filter (start, end) in YYYY-MM format
extend_historical_to
Opt-in ``(end_year, end_month)``. When set, each converted CMIP7 file has
its time axis relabelled so historical coverage ends on that month, letting
us fabricate CMIP7 data for years without real CMIP6 source data (e.g. the
fire diagnostic's 2002-2021 window). Defaults to ``None`` (time axis
untouched), so other CMIP7 conversions are unchanged.
Opt-in ``(end_year, end_month)``. When set, a converted CMIP7 file that stops
short has its final year repeated until it ends on that month, letting us
fabricate CMIP7 data for years without real CMIP6 source data (e.g. the fire
diagnostic's 2002-2021 window). Files that already reach it are untouched, as
are all conversions when this defaults to ``None``.
"""
self.slug = slug
self.facets = facets
Expand Down
Loading
Loading