From 9e767428f6e0c165b451a742433617ae5109d3d0 Mon Sep 17 00:00:00 2001 From: maritsandstad Date: Fri, 13 Mar 2026 09:27:42 +0100 Subject: [PATCH 01/11] Fix to ensure path exists for netcdf dump (cherry picked from commit fc591cac8a2f6d926bc4ff9d5802342ee70b53b1) --- src/meteor/ensemble_output.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/meteor/ensemble_output.py b/src/meteor/ensemble_output.py index e6de46d..507e3e7 100644 --- a/src/meteor/ensemble_output.py +++ b/src/meteor/ensemble_output.py @@ -5,6 +5,8 @@ climate variables, spatial aggregations, and impact metrics. """ +from pathlib import Path + import xarray as xr @@ -177,6 +179,8 @@ def to_netcdf(self, path, include_impacts=True): combined.attrs[key] = value # Save to netCDF + path = Path(path) + path.parent.mkdir(exist_ok=True, parents=True) combined.to_netcdf(path) print(f"✅ Saved ensemble to {path}") From b90204c4541ea131c3e0700bf9701e595236281d Mon Sep 17 00:00:00 2001 From: maritsandstad Date: Fri, 13 Mar 2026 14:59:28 +0100 Subject: [PATCH 02/11] Various fixes for scaling timeseries and slightly more varried input data formats --- CHANGELOG.rst | 3 +++ src/meteor/geo_data_utils.py | 29 +++++++++++++++++++++++++++++ src/meteor/meteor_interface.py | 7 ++++++- src/meteor/noise_generator.py | 7 +++---- tests/unit/test_geo_data_utils.py | 20 ++++++++++++++++++++ 5 files changed, 61 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 14e57c7..8cae05c 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -17,6 +17,9 @@ The changes listed in this file are categorised as follows: [Unreleased] --------------------- +### Changed +- Now possibly to send variable length temperature scaling timeseries, fixed noise generator for wrong ordering of base data dimensions + [Version 1.6.0] ----------------------------- diff --git a/src/meteor/geo_data_utils.py b/src/meteor/geo_data_utils.py index b493741..a999b0a 100644 --- a/src/meteor/geo_data_utils.py +++ b/src/meteor/geo_data_utils.py @@ -529,3 +529,32 @@ def extend_temeperature_anomaly_timeseries_for_scaling( dims=annual_temp_prediction_anomaly_gm.dims, ) return temperature_input_anomaly_extended + + +def find_time_dim_and_cut(base_array, n_time, n_lat, n_lon): + """ + Base_array should have dimensions (time, lat, lon), but the ordering can be off + """ + base_shape = base_array.shape + lat_dim = None + lon_dim = None + time_dim = 0 + for i, dim_size in enumerate(base_shape): + if dim_size == n_lat and lat_dim is None: + lat_dim = i + elif dim_size == n_lon: + lon_dim = i + else: + time_dim = i + if lat_dim is None or lon_dim is None: + raise ValueError( + f"Couldn't find lat/lon dimensions in base array with shape {base_shape}" + ) + if base_shape[time_dim] != n_time: + if time_dim == 0: + base_array = base_array[:n_time, :, :] + elif time_dim == 1: + base_array = base_array[:, :n_time, :] + elif time_dim == 2: + base_array = base_array[:, :, :n_time] + return base_array.transpose(time_dim, lat_dim, lon_dim) diff --git a/src/meteor/meteor_interface.py b/src/meteor/meteor_interface.py index a898733..12f3774 100644 --- a/src/meteor/meteor_interface.py +++ b/src/meteor/meteor_interface.py @@ -786,6 +786,9 @@ def _get_or_compute_pattern_scaling( base_year = 1750 # Default assumption if temp_scaling_ts is not None: + print( + f" Ts scaleing before calculating annual pred {temp_scaling_ts['year']}" + ) annual_prediction = self._compute_timeseries_scaling( variable, annual_prediction, @@ -914,7 +917,9 @@ def _compute_timeseries_scaling( annual_temp_prediction_gm_anomaly = global_mean( annual_temp_prediction - annual_temp_prediction_base ) - temperature_input_anomaly = temp_scaling_ts - temperature_input_base + print(temp_scaling_ts.shape) + print(temperature_input_base.shape, annual_temp_prediction_gm_anomaly.shape) + temperature_input_anomaly = temp_scaling_ts - temperature_input_base.values if len(temperature_input_anomaly) != len(annual_temp_prediction_gm_anomaly): temperature_input_anomaly = ( extend_temeperature_anomaly_timeseries_for_scaling( diff --git a/src/meteor/noise_generator.py b/src/meteor/noise_generator.py index 4275819..8eee5a1 100644 --- a/src/meteor/noise_generator.py +++ b/src/meteor/noise_generator.py @@ -23,7 +23,7 @@ from sklearn.linear_model import LinearRegression from statsmodels.tsa.api import VAR -from .geo_data_utils import global_mean +from .geo_data_utils import find_time_dim_and_cut, global_mean class MeteorNoiseGenerator: @@ -553,9 +553,8 @@ def generate_realization( # Now reshape to (n_time, n_lat, n_lon) # If the time dimension doesn't match, select the first n_time steps - if base_values.shape[0] != n_time: - base_values = base_values[:n_time, :, :] - + print(n_lat, n_lon, n_time) + base_values = find_time_dim_and_cut(base_values, n_time, n_lat, n_lon) base_clim_np = base_values.reshape(n_time, n_lat, n_lon) # Generate realizations (only stochastic component varies) diff --git a/tests/unit/test_geo_data_utils.py b/tests/unit/test_geo_data_utils.py index 6ece028..000d95d 100644 --- a/tests/unit/test_geo_data_utils.py +++ b/tests/unit/test_geo_data_utils.py @@ -266,3 +266,23 @@ def test_extend_temperature_anomaly_timeseries_for_scaling(): extended[4] == 0.5 ) # Year after last should be same as annual_temp_prediction_gm_anomaly assert extended.coords["year"].values.tolist() == target_years.tolist() + + +def test_find_time_dim_and_cut(): + # Create a DataArray with time dimension and extra dimensions + data = np.random.rand(10, 5, 5) + n_time = 8 + n_lat = 5 + n_lon = 5 + + cut_data = geo_data_utils.find_time_dim_and_cut(data, n_time, n_lat, n_lon) + assert cut_data.shape == (n_time, n_lat, n_lon) + # Test with time as second dimension + data2 = np.random.rand(5, 10, 5) + cut_data2 = geo_data_utils.find_time_dim_and_cut(data2, n_time, n_lat, n_lon) + assert cut_data2.shape == (n_time, n_lat, n_lon) + + # Test with time as third dimension + data3 = np.random.rand(5, 5, 10) + cut_data3 = geo_data_utils.find_time_dim_and_cut(data3, n_time, n_lat, n_lon) + assert cut_data3.shape == (n_time, n_lat, n_lon) From 5fb128e0f8a2ea4589426debcd1ad9cb2e59604f Mon Sep 17 00:00:00 2001 From: benmsanderson Date: Thu, 11 Jun 2026 11:42:31 +0200 Subject: [PATCH 03/11] Preserve low-frequency global variability: default to pure VAR (use_exog='none') The EOF-weighting fix (PR #74) corrected the *monthly* global tas variance, but the *annual/decadal* global variability was still ~2.5x too small: the generated noise was temporally white/anti-persistent (lag-1 autocorr ~-0.1) whereas real global-mean tas is strongly red (lag-1 ~+0.57). Root cause is the VAR-X exogenous regressor, not the EOF basis or the seasonal t_glob removal. Stage-by-stage (CanESM5 tas, validated against piControl): target (piControl) annual std 0.089 K lag-1 +0.57 training anomaly global mean 0.081 K +0.56 (red) in-sample retained PCs -> global 0.074 K +0.51 (red) VAR-X generated, use_exog='all' 0.038 K +0.18 (white) VAR-X generated, use_exog='none' 0.076 K +0.53 (red) Using the smoothed global temperature (t_glob) as a VAR-X exogenous regressor absorbs the persistent low-frequency global variability into the deterministic forced term. At generation t_glob is the prescribed (smooth, internally invariant) trajectory, so that power is never regenerated -> white noise. With use_exog='none' the persistence stays in the AR dynamics and is reproduced (CanESM5 86% of target annual std, ACCESS-ESM1-5 79%; monthly variance unchanged). The temperature-dependent mean/seasonal response is already captured by the seasonal model, so the exog regressor is redundant. Changes: - MeteorNoiseGenerator and train_*_from_cmip6 default use_exog 'temp_only'->'none' - interface _get_default_config: tas 'all'->'none' (pr was already 'none') - 'temp_only'/'all' retained for backward compatibility (documented as suppressing low-frequency global variability; 'temp_only' is also unstable at high lag_order) Co-Authored-By: Claude Opus 4.8 (1M context) (cherry picked from commit 75759b3650a8b0a3a1e5fbcff32e6f284cd67062) --- src/meteor/meteor_interface.py | 9 +++++- src/meteor/noise_generator.py | 44 ++++++++++++++++++++--------- tests/unit/test_meteor_interface.py | 4 ++- tests/unit/test_noise_generator.py | 3 ++ 4 files changed, 45 insertions(+), 15 deletions(-) diff --git a/src/meteor/meteor_interface.py b/src/meteor/meteor_interface.py index 12f3774..0464249 100644 --- a/src/meteor/meteor_interface.py +++ b/src/meteor/meteor_interface.py @@ -45,7 +45,14 @@ def _get_default_config(variable): # Variable-specific defaults if variable == "tas": - config["use_exog"] = "all" + # 'none' (pure VAR): using t_glob as a VAR-X exogenous regressor absorbs + # the persistent low-frequency global variability into the deterministic + # forced term, so it is lost at generation (the prescribed trajectory has + # no internal variability) -- the noise becomes white and annual/decadal + # global variance collapses ~2.5x. The temperature-dependent mean/seasonal + # response is already captured by the seasonal model, so the exog is + # redundant here. See MeteorNoiseGenerator.use_exog. + config["use_exog"] = "none" config["transform"] = False elif variable == "pr": config["use_exog"] = "none" diff --git a/src/meteor/noise_generator.py b/src/meteor/noise_generator.py index 8eee5a1..cc0d8f0 100644 --- a/src/meteor/noise_generator.py +++ b/src/meteor/noise_generator.py @@ -55,7 +55,7 @@ class MeteorNoiseGenerator: Whether the model has been fitted """ - def __init__(self, n_modes=40, lag_order=2, use_exog="temp_only", weight_eofs=True): + def __init__(self, n_modes=40, lag_order=2, use_exog="none", weight_eofs=True): """ Initialize the noise generator. @@ -69,11 +69,23 @@ def __init__(self, n_modes=40, lag_order=2, use_exog="temp_only", weight_eofs=Tr less sensitive to truncation). lag_order : int, default 2 Lag order for VARX model - use_exog : str, default 'temp_only' + use_exog : str, default 'none' Exogenous variables to use in VARX model: - - 'all': Use temperature, annual_cos, annual_sin (original behavior) - - 'temp_only': Use only temperature (recommended to avoid spurious seasonality) - - 'none': Pure VAR with no exogenous variables + - 'none': Pure VAR with no exogenous variables (recommended). + - 'temp_only': Use only temperature. + - 'all': Use temperature, annual_cos, annual_sin (original behavior). + + 'none' is the default because using the smoothed global temperature + (t_glob) as an exogenous regressor absorbs the persistent + low-frequency global variability into the deterministic forced term. + At generation t_glob is the prescribed (smooth, internally + invariant) trajectory, so that power is not regenerated: the noise + becomes temporally white and annual/decadal global-mean variance + collapses (~2.5x too small for tas). The temperature-dependent + mean and seasonal response is already captured by the seasonal + model, so the exog regressor is redundant as well as harmful to + internal variability. 'temp_only'/'all' are retained for backward + compatibility but suppress low-frequency global variability. weight_eofs : bool, default True If True, area-weight the anomaly field by sqrt(cos(latitude)) before fitting the EOF/PCA basis, so PCA optimizes area-weighted variance @@ -1180,7 +1192,7 @@ def train_noise_model_from_cmip6( custom_global_temp=None, use_picontrol_baseline=True, save_diagnostics=False, - use_exog="temp_only", + use_exog="none", weight_eofs=True, verbose=False, ): @@ -1221,11 +1233,14 @@ def train_noise_model_from_cmip6( model.diagnostic_X_features, model.diagnostic_t_glob, model.diagnostic_time, model.diagnostic_seasonal_coef, model.diagnostic_seasonal_intercept, and model.diagnostic_Y_data. - use_exog : str, default 'temp_only' + use_exog : str, default 'none' Exogenous variables to use in VARX model: + - 'none': Pure VAR with no exogenous variables (recommended; preserves + low-frequency global variability) + - 'temp_only': Use only temperature - 'all': Use temperature, annual_cos, annual_sin (may cause spurious seasonality) - - 'temp_only': Use only temperature (recommended) - - 'none': Pure VAR with no exogenous variables + See MeteorNoiseGenerator for why t_glob as an exog regressor suppresses + global-mean internal variability. weight_eofs : bool, default True If True, area-weight the anomaly field by sqrt(cos(latitude)) before fitting the EOF basis so global-mean variability is preserved. See @@ -1305,7 +1320,7 @@ def train_multiple_noise_models_from_cmip6( cache_dir=None, custom_global_temp=None, use_picontrol_baseline=True, - use_exog="temp_only", + use_exog="none", weight_eofs=True, ): """ @@ -1337,11 +1352,14 @@ def train_multiple_noise_models_from_cmip6( Whether to use piControl data as baseline for temperature anomalies. This ensures consistency with pattern scaling. If False, falls back to using first 42 years of training data. - use_exog : str, default 'temp_only' + use_exog : str, default 'none' Exogenous variables to use in VARX model: + - 'none': Pure VAR with no exogenous variables (recommended; preserves + low-frequency global variability) + - 'temp_only': Use only temperature - 'all': Use temperature, annual_cos, annual_sin (may cause spurious seasonality) - - 'temp_only': Use only temperature (recommended) - - 'none': Pure VAR with no exogenous variables + See MeteorNoiseGenerator for why t_glob as an exog regressor suppresses + global-mean internal variability. weight_eofs : bool, default True If True, area-weight the anomaly field by sqrt(cos(latitude)) before fitting the EOF basis so global-mean variability is preserved. See diff --git a/tests/unit/test_meteor_interface.py b/tests/unit/test_meteor_interface.py index de3625d..a022f0d 100644 --- a/tests/unit/test_meteor_interface.py +++ b/tests/unit/test_meteor_interface.py @@ -58,7 +58,9 @@ def test_get_default_config(): assert tas_config["n_modes_noise"] == 40 assert tas_config["lag_order"] == 2 assert tas_config["training_scenario"] == "ssp245" - assert tas_config["use_exog"] == "all" + # 'none': t_glob as a VAR-X exog regressor whitens the noise and collapses + # low-frequency global variability (see MeteorNoiseGenerator.use_exog) + assert tas_config["use_exog"] == "none" assert not tas_config["transform"] pr_config = _get_default_config("pr") diff --git a/tests/unit/test_noise_generator.py b/tests/unit/test_noise_generator.py index dd12c0f..e0dde0b 100644 --- a/tests/unit/test_noise_generator.py +++ b/tests/unit/test_noise_generator.py @@ -26,6 +26,9 @@ def test_meteor_noise_generator_initialization(): assert generator.n_modes == 40 # default value assert generator.lag_order == 2 # default value assert generator.weight_eofs is True # area-weighting on by default + # pure VAR by default: t_glob exog whitens the noise and collapses + # low-frequency global variability + assert MeteorNoiseGenerator().use_exog == "none" # Test error handling with uninitiated state with pytest.raises(ValueError, match="Invalid use_exog value:"): From 585614aaba2aaa90821fcedd4e5ba5620ff98bdd Mon Sep 17 00:00:00 2001 From: maritsandstad Date: Thu, 11 Jun 2026 23:20:04 +0200 Subject: [PATCH 04/11] Fixing scaling (cherry picked from commit 4f5df801a0ca3efeb083ee797e9745fc9273092f) --- src/meteor/geo_data_utils.py | 22 ++++++++++++++++------ src/meteor/meteor_interface.py | 4 ---- 2 files changed, 16 insertions(+), 10 deletions(-) diff --git a/src/meteor/geo_data_utils.py b/src/meteor/geo_data_utils.py index a999b0a..fefb7aa 100644 --- a/src/meteor/geo_data_utils.py +++ b/src/meteor/geo_data_utils.py @@ -34,6 +34,18 @@ def get_time_name(ds): return time_name raise RuntimeError("Couldn't find a time coordinate") +def get_year_series(ds: xr.Dataset) -> np.ndarray: + """Extract a DataFrame of just the year columns from a larger DataFrame.""" + time_name = get_time_name(ds) + if time_name == "year": + return ds[time_name].values + elif time_name == "time": + return ds[time_name].dt.year.values + elif time_name == "month": + return ds[time_name].dt.year.values + (ds[time_name].dt.month.values - 1) / 12.0 + else: + raise ValueError(f"Unexpected time dimension name '{time_name}'") + def get_lat_name(ds): """ @@ -509,12 +521,8 @@ def extend_temeperature_anomaly_timeseries_for_scaling( temperature_input_anomaly : xarray.DataArray Temperature anomaly time series used for scaling (relative to base year) """ - years_prediction = annual_temp_prediction_anomaly_gm[ - get_time_name(annual_temp_prediction_anomaly_gm) - ].values - years_input = temperature_input_anomaly[ - get_time_name(temperature_input_anomaly) - ].values + years_prediction = get_year_series(annual_temp_prediction_anomaly_gm) + years_input = get_year_series(temperature_input_anomaly) temperature_input_anomaly_extended = np.copy( annual_temp_prediction_anomaly_gm.values ) @@ -523,11 +531,13 @@ def extend_temeperature_anomaly_timeseries_for_scaling( temperature_input_anomaly_extended[i] = temperature_input_anomaly.sel( {get_time_name(temperature_input_anomaly): year} ).values + temperature_input_anomaly_extended = xr.DataArray( data=temperature_input_anomaly_extended, coords={get_time_name(annual_temp_prediction_anomaly_gm): years_prediction}, dims=annual_temp_prediction_anomaly_gm.dims, ) + return temperature_input_anomaly_extended diff --git a/src/meteor/meteor_interface.py b/src/meteor/meteor_interface.py index 0464249..75068b1 100644 --- a/src/meteor/meteor_interface.py +++ b/src/meteor/meteor_interface.py @@ -924,8 +924,6 @@ def _compute_timeseries_scaling( annual_temp_prediction_gm_anomaly = global_mean( annual_temp_prediction - annual_temp_prediction_base ) - print(temp_scaling_ts.shape) - print(temperature_input_base.shape, annual_temp_prediction_gm_anomaly.shape) temperature_input_anomaly = temp_scaling_ts - temperature_input_base.values if len(temperature_input_anomaly) != len(annual_temp_prediction_gm_anomaly): temperature_input_anomaly = ( @@ -934,8 +932,6 @@ def _compute_timeseries_scaling( temperature_input_anomaly, ) ) - print(temp_scaling_ts.shape) - print(temperature_input_anomaly.shape, annual_temp_prediction_gm_anomaly.shape) temp_scaling = np.where( annual_temp_prediction_gm_anomaly.values != 0, temperature_input_anomaly.values / annual_temp_prediction_gm_anomaly.values, From 16ad71736c7ca6f3c95d06fdbc196419f2fd953c Mon Sep 17 00:00:00 2001 From: maritsandstad Date: Fri, 12 Jun 2026 10:08:03 +0200 Subject: [PATCH 05/11] Tests, fixes for more diverse data (cherry picked from commit 25f0686c887976692247459b0503eca57d4f9e46) --- CHANGELOG.rst | 7 +++- src/meteor/geo_data_utils.py | 14 +++++-- tests/unit/test_geo_data_utils.py | 67 +++++++++++++++++++++++++++++++ 3 files changed, 83 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 8cae05c..8cac9b6 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -18,7 +18,12 @@ The changes listed in this file are categorised as follows: --------------------- ### Changed -- Now possibly to send variable length temperature scaling timeseries, fixed noise generator for wrong ordering of base data dimensions + +- Now possibly to send variable length temperature scaling timeseries, fixed noise generator for wrong ordering of base data dimensions + +### Fixed + +- Variable length timeseries now works also when don't have "year" as time dimension. [Version 1.6.0] diff --git a/src/meteor/geo_data_utils.py b/src/meteor/geo_data_utils.py index fefb7aa..743c849 100644 --- a/src/meteor/geo_data_utils.py +++ b/src/meteor/geo_data_utils.py @@ -34,17 +34,17 @@ def get_time_name(ds): return time_name raise RuntimeError("Couldn't find a time coordinate") + def get_year_series(ds: xr.Dataset) -> np.ndarray: """Extract a DataFrame of just the year columns from a larger DataFrame.""" time_name = get_time_name(ds) if time_name == "year": return ds[time_name].values - elif time_name == "time": + if time_name == "time": return ds[time_name].dt.year.values - elif time_name == "month": + if time_name == "month": return ds[time_name].dt.year.values + (ds[time_name].dt.month.values - 1) / 12.0 - else: - raise ValueError(f"Unexpected time dimension name '{time_name}'") + raise ValueError(f"Unexpected time dimension name '{time_name}'") def get_lat_name(ds): @@ -526,6 +526,12 @@ def extend_temeperature_anomaly_timeseries_for_scaling( temperature_input_anomaly_extended = np.copy( annual_temp_prediction_anomaly_gm.values ) + input_time_name = get_time_name(temperature_input_anomaly) + if input_time_name != "year": + temperature_input_anomaly = temperature_input_anomaly.rename( + {input_time_name: "year"} + ) + temperature_input_anomaly.coords["year"] = years_input for i, year in enumerate(years_prediction): if year in years_input: temperature_input_anomaly_extended[i] = temperature_input_anomaly.sel( diff --git a/tests/unit/test_geo_data_utils.py b/tests/unit/test_geo_data_utils.py index 000d95d..ecc3a52 100644 --- a/tests/unit/test_geo_data_utils.py +++ b/tests/unit/test_geo_data_utils.py @@ -1,4 +1,5 @@ import numpy as np +import pandas as pd import pytest import xarray as xr @@ -266,6 +267,72 @@ def test_extend_temperature_anomaly_timeseries_for_scaling(): extended[4] == 0.5 ) # Year after last should be same as annual_temp_prediction_gm_anomaly assert extended.coords["year"].values.tolist() == target_years.tolist() + temp_anomaly_time = temp_anomaly.rename({"year": "time"}) + temp_anomaly_time.coords["time"] = pd.to_datetime( + { + "year": temp_anomaly_time.coords["time"].values, + "month": [6] * 3, + "day": [30] * 3, + } + ).values + + print( + pd.to_datetime({"year": target_years, "month": [6] * 5, "day": [30] * 5}).values + ) + print( + annual_temp_prediction_gm_anomaly.rename({"year": "time"}).assign_coords( + time=pd.to_datetime( + {"year": target_years, "month": [6] * 5, "day": [30] * 5} + ).values + ) + ) + extended_time = geo_data_utils.extend_temeperature_anomaly_timeseries_for_scaling( + annual_temp_prediction_gm_anomaly.rename({"year": "time"}).assign_coords( + time=pd.to_datetime( + {"year": target_years, "month": [6] * 5, "day": [30] * 5} + ).values + ), + temp_anomaly_time, + ) + assert extended_time.shape == (5,) + assert extended_time[1] == 0.1 + assert extended_time[0] == 0.1 + assert extended_time[2] == 0.2 + assert extended_time[3] == 0.3 + assert extended_time[4] == 0.5 + + +def test_get_year_series(): + # Create a simple dataset with different time coordinate names + data = np.zeros(5) + ds_year = xr.Dataset( + {"var": (("year"), data)}, + coords={"year": np.arange(2000, 2005)}, + ) + ds_time = xr.Dataset( + {"var": (("time"), data)}, + coords={ + "time": pd.to_datetime( + {"year": range(2000, 2005), "month": [6] * 5, "day": [30] * 5} + ).values + }, + ) + print(ds_time["time"].dt) + ds_month = ds_year.rename({"year": "month"}) + ds_month.coords["month"] = pd.to_datetime( + {"year": range(2000, 2005), "month": [6] * 5, "day": [30] * 5} + ).values # Monthly data + + years_from_year = geo_data_utils.get_year_series(ds_year) + years_from_time = geo_data_utils.get_year_series(ds_time) + years_from_month = geo_data_utils.get_year_series(ds_month) + + assert years_from_year.shape == (5,) + assert years_from_time.shape == (5,) + assert years_from_month.shape == (5,) + assert np.array_equal(years_from_year, np.arange(2000, 2005)) + assert np.array_equal(years_from_time, np.arange(2000, 2005)) + assert not np.array_equal(years_from_month, np.arange(2000, 2005)) def test_find_time_dim_and_cut(): From b3e1b75125a0742975f3ff65d42d0d80af9379db Mon Sep 17 00:00:00 2001 From: maritsandstad Date: Sat, 13 Jun 2026 22:35:28 +0200 Subject: [PATCH 06/11] interface field generation fix on base (cherry picked from commit e16947bbd09d434b99198b417171ab148b4579e7) --- src/meteor/meteor_interface.py | 728 ++++++++++++++++++++-------- src/meteor/noise_generator.py | 27 +- tests/unit/test_meteor_interface.py | 379 ++++++++++++++- 3 files changed, 918 insertions(+), 216 deletions(-) diff --git a/src/meteor/meteor_interface.py b/src/meteor/meteor_interface.py index 75068b1..560f8a6 100644 --- a/src/meteor/meteor_interface.py +++ b/src/meteor/meteor_interface.py @@ -6,6 +6,8 @@ """ import os +from dataclasses import dataclass +from typing import Any import numpy as np import xarray as xr @@ -66,6 +68,90 @@ def _get_default_config(variable): return config +@dataclass +class PatternScalingResult: + """ + Result of a pattern scaling computation. + + Holds both the time-sliced pattern scaling output used directly by the + generators and the full-trajectory warming plus slice bookkeeping needed to + generate spun-up stochastic PCs over the full trajectory and slice them to + the requested output window. + + Attributes + ---------- + monthly_prediction : xr.DataArray + Monthly pattern prediction sliced to the requested output window. + monthly_warming : np.ndarray + Global-mean monthly warming sliced to the requested output window. + em_data : Any + Emissions data used to drive the pattern model. + conc_data : Any + Concentration data used to drive the pattern model. + full_monthly_warming : np.ndarray + Global-mean monthly warming over the full (un-sliced) trajectory. Used to + generate stochastic PCs so the autoregressive spin-up transient is parked + at the trajectory start rather than inside the output window. + base_year : int + First year of the full monthly trajectory (origin for all month indexing). + start_month_idx : int + Month index (relative to ``base_year``) of the first output month. + end_month_idx : int + Month index (relative to ``base_year``) one past the last output month. + """ + + monthly_prediction: xr.DataArray + monthly_warming: np.ndarray + em_data: Any + conc_data: Any + full_monthly_warming: np.ndarray + base_year: int + start_month_idx: int + end_month_idx: int + + +@dataclass +class GenerationInputs: + """ + Shared inputs prepared once per variable and consumed by both generators. + + Produced by :meth:`MeteorInterface._prepare_generation` and passed to both + ``_generate_timeseries`` and ``_generate_gridded`` so the two paths share the + same pattern scaling and the same spun-up stochastic PC realisations. + + Attributes + ---------- + pattern : PatternScalingResult + Pattern scaling result (sliced prediction/warming + slice bookkeeping). + stochastic_pcs : np.ndarray or None + Stochastic PCs generated once over the full trajectory and sliced to the + output window. Shape ``(n_realizations, n_months, n_modes)``. ``None`` when + ``include_noise`` is False. + """ + + pattern: PatternScalingResult + stochastic_pcs: Any + + +def _stack_realizations(realizations): + """ + Stack a list of per-realization DataArrays along a ``realization`` dimension. + + Parameters + ---------- + realizations : list of xr.DataArray + One DataArray per ensemble member. + + Returns + ------- + xr.DataArray + Concatenated array with a leading ``realization`` dimension. + """ + if len(realizations) > 1: + return xr.concat(realizations, dim="realization") + return realizations[0].expand_dims(realization=[0]) + + class MeteorInterface: """ High-level interface for training and generating METEOR emulators. @@ -596,6 +682,24 @@ def generate_ensemble_outputs( var_output = VariableOutput(variable) + # Prepare shared generation inputs ONCE when both output types are + # requested, so the time series and gridded paths are driven by the + # same pattern scaling and the same spun-up stochastic PCs (mutually + # consistent noise). When only one type is requested there is nothing + # to be consistent with, so that generator builds its own inputs. + gen_inputs = None + if timeseries and gridded: + gen_inputs = self._prepare_generation( + variable, + scenario, + start_year, + end_year, + n_realizations, + include_noise=include_noise, + temp_scaling_ts=temp_scaling_ts, + verbose=verbose, + ) + # Generate timeseries if requested if timeseries: if verbose: # pragma: no cover @@ -607,6 +711,7 @@ def generate_ensemble_outputs( end_year, n_realizations, timeseries, + gen_inputs=gen_inputs, custom_regions=custom_regions, include_noise=include_noise, temp_scaling_ts=temp_scaling_ts, @@ -624,6 +729,7 @@ def generate_ensemble_outputs( end_year, n_realizations, gridded, + gen_inputs=gen_inputs, include_noise=include_noise, temp_scaling_ts=temp_scaling_ts, verbose=verbose, @@ -700,8 +806,10 @@ def _get_or_compute_pattern_scaling( Returns ------- - tuple - (monthly_prediction_sliced, monthly_warming_sliced, em_data, conc_data) + PatternScalingResult + Sliced monthly prediction/warming plus the full-trajectory warming and + slice bookkeeping (``base_year``, ``start_month_idx``, ``end_month_idx``) + needed to generate spun-up stochastic PCs aligned to the output window. """ # Parse scenario input scenario_info = parse_scenario_input(scenario) @@ -834,7 +942,16 @@ def _get_or_compute_pattern_scaling( ) monthly_warming_sliced = full_monthly_warming[start_month_idx:end_month_idx] - return monthly_prediction_sliced, monthly_warming_sliced, em_data, conc_data + return PatternScalingResult( + monthly_prediction=monthly_prediction_sliced, + monthly_warming=monthly_warming_sliced, + em_data=em_data, + conc_data=conc_data, + full_monthly_warming=full_monthly_warming, + base_year=base_year, + start_month_idx=start_month_idx, + end_month_idx=end_month_idx, + ) # TODO - possibly add verbosity? def _compute_timeseries_scaling( @@ -945,61 +1062,50 @@ def _compute_timeseries_scaling( ) return annual_prediction_base + annual_prediction_anomaly * temp_scaling - def _generate_timeseries( + def _prepare_generation( self, variable, scenario, start_year, end_year, n_realizations, - aggregations, - custom_regions=None, include_noise=True, temp_scaling_ts=None, verbose=True, ): """ - Generate time series outputs with spatial aggregations. + Prepare the shared inputs consumed by both output generators. + + Computes pattern scaling once and generates the stochastic PCs once over + the FULL trajectory (so the autoregressive spin-up transient is parked at + the trajectory start, not inside the output window), then slices the PCs to + the requested output window. The resulting :class:`GenerationInputs` is + passed to both ``_generate_timeseries`` and ``_generate_gridded`` so the two + paths are driven by identical pattern scaling and identical noise draws. Parameters ---------- variable : str - Climate variable to generate - scenario : str - Emission scenario ('ssp245', 'ssp585', etc.) - start_year : int - Start year - end_year : int - End year (inclusive) + Climate variable ('tas', 'pr'). + scenario : str or dict + Scenario specification (see :meth:`generate_ensemble_outputs`). + start_year, end_year : int + Output window (inclusive). n_realizations : int - Number of ensemble members - aggregations : list of str - Spatial aggregations to compute: - - 'global': Global mean - - 'regional:CODE': AR6 region mean (e.g., 'regional:NEU') - - 'regional:custom:NAME': Custom region (requires custom_regions dict) - - 'point:LAT,LON': Single grid point (e.g., 'point:59.9,10.8') - custom_regions : dict, optional - Custom region definitions: {'name': {'lat': (min, max), 'lon': (min, max)}} + Number of ensemble members. include_noise : bool - If True, add stochastic noise; if False, return forced response only + If False, no PCs are generated (climatology only). + temp_scaling_ts : xr.DataArray, optional + Optional global-mean temperature trajectory to scale the pattern to. verbose : bool - Print progress messages + Print progress messages. Returns ------- - dict - Dictionary mapping aggregation names to xarray DataArrays - with shape (n_realizations, n_months) + GenerationInputs + Shared pattern scaling result and (window-sliced) stochastic PCs. """ - # Always use the TRAINING scenario for transform fitting, not the - # prediction scenario. default is ssp245 - transform_training_scenario = self._training_config.get(variable, {}).get( - "training_scenario", "ssp245" - ) - - # Get pattern scaling results - pattern_result = self._get_or_compute_pattern_scaling( + pattern = self._get_or_compute_pattern_scaling( variable, scenario, start_year, @@ -1007,10 +1113,75 @@ def _generate_timeseries( temp_scaling_ts=temp_scaling_ts, verbose=verbose, ) - monthly_prediction = pattern_result[0] - monthly_warming = pattern_result[1] - # Get CMIP6 data for transform fitting + stochastic_pcs = None + if include_noise: + noise_model = self.noise_models[variable] + if verbose: # pragma: no cover + print( + f" → Generating {n_realizations} stochastic PC realizations " + "(full trajectory, spun-up)..." + ) + # Generate over the FULL trajectory so spin-up is resolved before the + # output window, then slice to the window on a January boundary + # (start_month_idx is always a multiple of 12). + full_pcs = noise_model.generate_stochastic_pcs( + pattern.full_monthly_warming, + n_realizations=n_realizations, + random_seed=None, + ) + if full_pcs.ndim == 2: + # Single realization -> add leading realization axis + full_pcs = full_pcs[np.newaxis, ...] + stochastic_pcs = full_pcs[ + :, pattern.start_month_idx : pattern.end_month_idx, : + ] + + return GenerationInputs(pattern=pattern, stochastic_pcs=stochastic_pcs) + + def _get_transform_config(self, variable): + """ + Return the resolved transform config for a variable, or None. + + Handles both the fitted case (stored as a dict with a ``'config'`` entry) + and the not-yet-fitted case (stored as a ``VariableTransformConfig``). + """ + transform_info = self.transforms.get(variable, None) + if isinstance(transform_info, dict): + return transform_info.get("config") + return transform_info + + def _load_transform_reference(self, variable, start_year, end_year, verbose=True): + """ + Load and prepare CMIP6 reference data for distribution-transform fitting. + + Only needed for variables that have a distribution transform (e.g. ``pr``). + Returns the gridded CMIP6 reference field sliced to the output window plus, + for precipitation, the gridded first-year baseline field. The timeseries + path aggregates these per requested region; the gridded path uses them + directly with the per-gridpoint (3D) transform. + + Parameters + ---------- + variable : str + Climate variable. + start_year, end_year : int + Output window (inclusive). + verbose : bool + Print progress messages. + + Returns + ------- + tuple + ``(ssp_data, pr_first_year_mean)`` where ``ssp_data`` is the gridded + reference field (anomalies for ``tas``, absolute for ``pr``) and + ``pr_first_year_mean`` is the gridded first-year baseline field for + ``pr`` (``None`` otherwise). + """ + transform_training_scenario = self._training_config.get(variable, {}).get( + "training_scenario", "ssp245" + ) + if verbose: # pragma: no cover print( f" → Loading CMIP6 training data for {transform_training_scenario}..." @@ -1019,47 +1190,37 @@ def _generate_timeseries( ["historical", transform_training_scenario], self.model, monthly=True )[variable] - # Load piControl data for baseline (used for temperature anomalies) if verbose: # pragma: no cover print(f" → Loading piControl baseline for {variable}...") picontrol_data = self.data_getter.make_meteor_training_data_composite( ["piControl"], self.model, monthly=True )[variable] - # Compute piControl climatology (mean across all time) picontrol_mean = picontrol_data.mean(dim="month") - # For precipitation: use first-year (2015) baseline instead of piControl - # This is because CMIP6 scenarios in 2015 already include ~1°C of historical - # warming effects on precipitation, so using piControl would create a ~2-4% bias. - # We use the first 12 months of the prediction period (start_year) as the baseline. - # Note: ssp_data is a composite starting from historical (~1850), so we need to - # find the correct index for start_year. + pr_first_year_mean = None + + # For precipitation: use first-year baseline instead of piControl. CMIP6 + # scenarios already include ~1°C of historical warming effects on + # precipitation, so piControl would create a ~2-4% bias. Use the first 12 + # months of the prediction period (start_year) as the baseline. ssp_data is + # a composite starting from historical (~1850), so find the index for + # start_year. if variable == "pr": - # Infer composite start year from the data length and structure - # Historical experiments in CMIP6 typically start at 1850 - # We can infer this from the data by checking if it includes historical n_months = len(ssp_data.month) - # Check if ssp_data has a 'start_year' attribute (set by data getter) - # Otherwise infer from experiment structure if hasattr(ssp_data, "start_year"): composite_start_year = int(ssp_data.start_year) else: # Default assumption: historical+scenario composite starts at 1850 - # If the data is shorter than expected, calculate backwards from end_year expected_months_from_1850 = (end_year - 1850 + 1) * 12 if n_months < expected_months_from_1850: - # Data is shorter - calculate start year from data length composite_start_year = end_year - (n_months // 12) + 1 else: composite_start_year = 1850 start_year_idx = (start_year - composite_start_year) * 12 - end_year_idx = ( - end_year - composite_start_year + 1 - ) * 12 # +1 for inclusive + end_year_idx = (end_year - composite_start_year + 1) * 12 # inclusive - # Validate indices are within bounds - fail loudly if not composite_end_year = composite_start_year + n_months // 12 - 1 if start_year_idx < 0: @@ -1083,9 +1244,9 @@ def _generate_timeseries( month=slice(start_year_idx, start_year_idx + 12) ).mean(dim="month") - # CRITICAL: Slice ssp_data to only the prediction period (start_year to end_year) - # for Gamma transform fitting. Using the full historical+scenario composite - # would result in a lower mean distribution, causing negative bias. + # CRITICAL: Slice ssp_data to only the prediction period for Gamma + # transform fitting. Using the full historical+scenario composite would + # result in a lower mean distribution, causing negative bias. ssp_data = ssp_data.isel(month=slice(start_year_idx, end_year_idx)) if verbose: # pragma: no cover @@ -1093,9 +1254,8 @@ def _generate_timeseries( f" → Using {start_year} baseline for PR instead of piControl" ) - # For temperature: convert CMIP6 to anomalies (pattern scaling outputs anomalies) - # For precipitation: keep CMIP6 as absolute values (for Gamma transform fitting) - # but we'll add first-year baseline to pattern output below + # For temperature: convert CMIP6 to anomalies (pattern scaling outputs + # anomalies). For precipitation keep absolute values for Gamma fitting. if variable == "tas": if verbose: # pragma: no cover print( @@ -1103,36 +1263,92 @@ def _generate_timeseries( ) ssp_data = ssp_data - picontrol_mean - # ✅ Generate stochastic PCs (or skip if climatology only) - noise_model = self.noise_models[variable] - stochastic_pcs = None + return ssp_data, pr_first_year_mean - if include_noise: - # CRITICAL: Generate stochastic PCs ONCE for all aggregations - # This ensures all spatial scales share the same underlying variability - if verbose: # pragma: no cover - print( - f" → Generating {n_realizations} stochastic PC realizations..." - ) + def _generate_timeseries( + self, + variable, + scenario, + start_year, + end_year, + n_realizations, + aggregations, + gen_inputs=None, + custom_regions=None, + include_noise=True, + temp_scaling_ts=None, + verbose=True, + ): + """ + Generate time series outputs with spatial aggregations. - stochastic_pcs = noise_model.generate_stochastic_pcs( - monthly_warming, - n_realizations=n_realizations, - random_seed=None, # Can expose this as parameter if needed + Parameters + ---------- + variable : str + Climate variable to generate + scenario : str + Emission scenario ('ssp245', 'ssp585', etc.) + start_year : int + Start year + end_year : int + End year (inclusive) + n_realizations : int + Number of ensemble members + aggregations : list of str + Spatial aggregations to compute: + - 'global': Global mean + - 'regional:CODE': AR6 region mean (e.g., 'regional:NEU') + - 'regional:custom:NAME': Custom region (requires custom_regions dict) + - 'point:LAT,LON': Single grid point (e.g., 'point:59.9,10.8') + custom_regions : dict, optional + Custom region definitions: {'name': {'lat': (min, max), 'lon': (min, max)}} + include_noise : bool + If True, add stochastic noise; if False, return forced response only + verbose : bool + Print progress messages + + Returns + ------- + dict + Dictionary mapping aggregation names to xarray DataArrays + with shape (n_realizations, n_months) + """ + # Build shared generation inputs if not provided by the caller. + if gen_inputs is None: + gen_inputs = self._prepare_generation( + variable, + scenario, + start_year, + end_year, + n_realizations, + include_noise=include_noise, + temp_scaling_ts=temp_scaling_ts, + verbose=verbose, ) - else: - if verbose: # pragma: no cover + + monthly_prediction = gen_inputs.pattern.monthly_prediction + monthly_warming = gen_inputs.pattern.monthly_warming + stochastic_pcs = gen_inputs.stochastic_pcs + noise_model = self.noise_models[variable] + + # Resolve transform and (only if needed) load CMIP6 reference data. + # tas has no transform, so its reference data is never loaded. + transform_config = self._get_transform_config(variable) + ssp_data = None + pr_first_year_mean = None + if transform_config and transform_config.transform_type: + ssp_data, pr_first_year_mean = self._load_transform_reference( + variable, start_year, end_year, verbose=verbose + ) + + if verbose: # pragma: no cover + if include_noise: + print(" → Using shared stochastic PC realizations") + else: print(" → Climatology only (no stochastic variability)") # Generate outputs for each aggregation results = {} - transform_info = self.transforms.get(variable, None) - - # Handle both dict (fitted) and VariableTransformConfig (not fitted) cases - if isinstance(transform_info, dict): - transform_config = transform_info.get("config") - else: - transform_config = transform_info # It's a VariableTransformConfig object for agg in aggregations: if verbose: # pragma: no cover @@ -1143,6 +1359,8 @@ def _generate_timeseries( # Pattern scaling outputs anomalies, but Gamma transform needs absolute values # We use first-year (2015) baseline instead of piControl to match CMIP6 starting point pr_baseline_agg = None + # CMIP6 reference aggregation, only needed/available when a transform exists + cmip6_agg = None if agg == "global": # Global mean @@ -1164,7 +1382,8 @@ def _generate_timeseries( else: # Climatology only: return pattern scaling with shape (1, time) raw_ensemble = pattern_agg[np.newaxis, :] - cmip6_agg = global_mean(ssp_data) + if ssp_data is not None: + cmip6_agg = global_mean(ssp_data) elif agg.startswith("regional:"): # Check if it's a custom region @@ -1188,7 +1407,8 @@ def _generate_timeseries( pattern_agg = regional_mean( monthly_prediction, region_mask=region_mask ).values - cmip6_agg = regional_mean(ssp_data, region_mask=region_mask) + if ssp_data is not None: + cmip6_agg = regional_mean(ssp_data, region_mask=region_mask) if variable == "pr": pr_baseline_agg = float( regional_mean( @@ -1234,7 +1454,8 @@ def _generate_timeseries( else: # Climatology only: return pattern scaling with shape (1, time) raw_ensemble = pattern_agg[np.newaxis, :] - cmip6_agg = regional_mean(ssp_data, region_code=region_code) + if ssp_data is not None: + cmip6_agg = regional_mean(ssp_data, region_code=region_code) elif agg.startswith("point:"): # Point extraction @@ -1262,7 +1483,8 @@ def _generate_timeseries( else: # Climatology only: return pattern scaling with shape (1, time) raw_ensemble = pattern_agg[np.newaxis, :] - cmip6_agg = extract_point(ssp_data, lat, lon) + if ssp_data is not None: + cmip6_agg = extract_point(ssp_data, lat, lon) else: raise ValueError(f"Unknown aggregation type: {agg}") @@ -1343,6 +1565,7 @@ def _generate_gridded( end_year, n_realizations, gridded_spec, + gen_inputs=None, include_noise=True, temp_scaling_ts=None, verbose=True, @@ -1364,41 +1587,62 @@ def _generate_gridded( Dictionary with keys 'annual', 'monthly', 'climatology' containing xarray DataArrays with gridded fields """ - # Get pattern scaling results (from cache or compute) - pattern_result = ( - self._get_or_compute_pattern_scaling( # pylint: disable=unused-variable + # Build shared generation inputs if not provided by the caller. The PCs + # are generated once over the full trajectory (spun-up) and sliced to the + # output window, then sliced again per output field below. + if gen_inputs is None: + gen_inputs = self._prepare_generation( variable, scenario, start_year, end_year, + n_realizations, + include_noise=include_noise, temp_scaling_ts=temp_scaling_ts, verbose=verbose, ) - ) - monthly_prediction = pattern_result[0] - monthly_warming = pattern_result[1] - # Get noise model + monthly_prediction = gen_inputs.pattern.monthly_prediction + monthly_warming = gen_inputs.pattern.monthly_warming + stochastic_pcs = gen_inputs.stochastic_pcs noise_model = self.noise_models[variable] - # Generate stochastic PCs (or skip if climatology only) - if include_noise: - if verbose: # pragma: no cover - print(f" → Generating {n_realizations} gridded realizations") - else: + if not include_noise: if verbose: # pragma: no cover print(" → Generating gridded climatology (no noise)") - n_realizations = 1 # Force to 1 for climatology + elif verbose: # pragma: no cover + print(" → Using shared stochastic PC realizations (gridded)") + + # Resolve transform and (only for transform variables, e.g. pr) load the + # CMIP6 reference field and fit the per-gridpoint target distribution once. + transform_config = self._get_transform_config(variable) + target_params = None + pr_baseline_field = None + if transform_config and transform_config.transform_type: + ssp_data, pr_baseline_field = self._load_transform_reference( + variable, start_year, end_year, verbose=verbose + ) + if verbose: # pragma: no cover + print( + f" → Fitting per-gridpoint {transform_config.transform_type} " + "target distribution..." + ) + target_params = transform_config.fit_3d_func( + ssp_data, transform_config.transform_type + ) # Extract requested time slices results = {} n_months = len(monthly_warming) - # Helper to convert year to month index + # Window-relative month index. monthly_prediction, monthly_warming and the + # sliced stochastic PCs all share the same origin (start_year), which is + # derived from base_year inside _get_or_compute_pattern_scaling, so all + # three stay aligned. def year_to_month_idx(year): return (year - start_year) * 12 - # Annual means + # Annual means (12-month average of each year) if "annual" in gridded_spec: if verbose: # pragma: no cover print( @@ -1409,38 +1653,19 @@ def year_to_month_idx(year): start_idx = year_to_month_idx(year) end_idx = start_idx + 12 if start_idx >= 0 and end_idx <= n_months: - # Generate realizations for this year - year_realizations = [] - for i in range(n_realizations): # pylint: disable=unused-variable - if include_noise: - # Generate full field with noise - realization = noise_model.generate_realization( - monthly_warming[start_idx:end_idx], - n_realizations=1, - noise_only=True, - add_base=monthly_prediction.isel( - month=slice(start_idx, end_idx) - ), - ) - else: - # Just use pattern scaling - realization = monthly_prediction.isel( - month=slice(start_idx, end_idx) - ) - - # Average over 12 months - annual_mean = realization.mean(dim="month") - year_realizations.append(annual_mean) - - # Stack realizations - if len(year_realizations) > 1: - annual_fields[year] = xr.concat( - year_realizations, dim="realization" - ) - else: - annual_fields[year] = year_realizations[0].expand_dims( - realization=[0] - ) + annual_fields[year] = self._generate_gridded_slice( + monthly_prediction, + monthly_warming, + noise_model, + stochastic_pcs, + start_idx, + end_idx, + include_noise, + reduce_time=True, + transform_config=transform_config, + target_params=target_params, + pr_baseline_field=pr_baseline_field, + ) else: if verbose: # pragma: no cover print( @@ -1448,7 +1673,7 @@ def year_to_month_idx(year): ) results["annual"] = annual_fields - # Monthly fields + # Monthly fields (all 12 months retained) if "monthly" in gridded_spec: if verbose: # pragma: no cover print( @@ -1459,34 +1684,19 @@ def year_to_month_idx(year): start_idx = year_to_month_idx(year) end_idx = start_idx + 12 if start_idx >= 0 and end_idx <= n_months: - # Generate realizations for this year - year_realizations = [] - for i in range(n_realizations): - if include_noise: - # Generate full field with noise - realization = noise_model.generate_realization( - monthly_warming[start_idx:end_idx], - n_realizations=1, - noise_only=True, - add_base=monthly_prediction.isel( - month=slice(start_idx, end_idx) - ), - ) - else: - # Just use pattern scaling - realization = monthly_prediction.isel( - month=slice(start_idx, end_idx) - ) - - year_realizations.append(realization) - - # Stack realizations (shape: realizations, month, lat, lon) - if len(year_realizations) > 1: - year_months = xr.concat(year_realizations, dim="realization") - else: - year_months = year_realizations[0].expand_dims(realization=[0]) - - monthly_fields[year] = year_months + monthly_fields[year] = self._generate_gridded_slice( + monthly_prediction, + monthly_warming, + noise_model, + stochastic_pcs, + start_idx, + end_idx, + include_noise, + reduce_time=False, + transform_config=transform_config, + target_params=target_params, + pr_baseline_field=pr_baseline_field, + ) else: if verbose: # pragma: no cover print( @@ -1507,38 +1717,21 @@ def year_to_month_idx(year): start_idx = year_to_month_idx(clim_start) end_idx = year_to_month_idx(clim_end + 1) # +1 to include end year if start_idx >= 0 and end_idx <= n_months: - # Generate realizations for this period - clim_realizations = [] - for i in range(n_realizations): - if include_noise: - # Generate full field with noise - realization = noise_model.generate_realization( - monthly_warming[start_idx:end_idx], - n_realizations=1, - noise_only=True, - add_base=monthly_prediction.isel( - month=slice(start_idx, end_idx) - ), - ) - else: - # Just use pattern scaling - realization = monthly_prediction.isel( - month=slice(start_idx, end_idx) - ) - - # Average over all months in period - clim_mean = realization.mean(dim="month") - clim_realizations.append(clim_mean) - - # Stack realizations - if len(clim_realizations) > 1: - climatology_fields[f"{clim_start}-{clim_end}"] = xr.concat( - clim_realizations, dim="realization" - ) - else: - climatology_fields[f"{clim_start}-{clim_end}"] = ( - clim_realizations[0].expand_dims(realization=[0]) + climatology_fields[f"{clim_start}-{clim_end}"] = ( + self._generate_gridded_slice( + monthly_prediction, + monthly_warming, + noise_model, + stochastic_pcs, + start_idx, + end_idx, + include_noise, + reduce_time=True, + transform_config=transform_config, + target_params=target_params, + pr_baseline_field=pr_baseline_field, ) + ) else: if verbose: # pragma: no cover print( @@ -1551,6 +1744,145 @@ def year_to_month_idx(year): return results + def _generate_gridded_slice( + self, + monthly_prediction, + monthly_warming, + noise_model, + stochastic_pcs, + start_idx, + end_idx, + include_noise, + reduce_time, + transform_config=None, + target_params=None, + pr_baseline_field=None, + ): + """ + Generate one gridded output slice (annual / monthly / climatology). + + Unifies the three previously-duplicated gridded loops. The only behavioural + differences between output types are the month range (``start_idx`` / + ``end_idx``, window-relative) and whether the time axis is averaged away + (``reduce_time``). + + Noise is taken from the shared spun-up PCs (sliced to this window) so that + gridded fields are consistent with the time series outputs and carry full + stationary variability. Any distribution transform (e.g. precipitation + Gamma) is applied to the MONTHLY field *before* time-averaging, since the + transform is a per-gridpoint quantile map over the sample axis. + + Parameters + ---------- + monthly_prediction : xr.DataArray + Window-sliced monthly pattern prediction (month, lat, lon). + monthly_warming : np.ndarray + Window-sliced global-mean monthly warming. + noise_model : object + Fitted noise model for this variable. + stochastic_pcs : np.ndarray or None + Window-sliced shared PCs (n_realizations, n_months, n_modes), or None + when ``include_noise`` is False. + start_idx, end_idx : int + Window-relative month indices bounding this slice. + include_noise : bool + Whether to add stochastic noise. + reduce_time : bool + If True, average over the month axis (annual / climatology); if False, + keep all months (monthly fields). + transform_config : VariableTransformConfig, optional + Distribution transform configuration (e.g. for ``pr``). + target_params : dict, optional + Pre-fitted per-gridpoint target distribution parameters. + pr_baseline_field : xr.DataArray, optional + Gridded first-year baseline (lat, lon) added to anomalies before the + transform to obtain absolute precipitation. + + Returns + ------- + xr.DataArray + Stacked realizations with a leading ``realization`` dimension. + """ + base_slice = monthly_prediction.isel(month=slice(start_idx, end_idx)) + + if include_noise: + pcs_slice = stochastic_pcs[:, start_idx:end_idx, :] + realizations = noise_model.generate_realization( + monthly_warming[start_idx:end_idx], + noise_only=True, + add_base=base_slice, + stochastic_pcs=pcs_slice, + ) + if not isinstance(realizations, list): + realizations = [realizations] + else: + realizations = [base_slice] + + ensemble = _stack_realizations(realizations) + + # Apply the distribution transform on the monthly field, before averaging. + if transform_config and transform_config.transform_type: + ensemble = self._apply_gridded_transform( + ensemble, transform_config, target_params, pr_baseline_field + ) + + if reduce_time: + ensemble = ensemble.mean(dim="month") + + return ensemble + + def _apply_gridded_transform( + self, ensemble, transform_config, target_params, pr_baseline_field + ): + """ + Apply a per-gridpoint distribution transform to a gridded ensemble. + + Mirrors the time series transform but uses the 3D (per-gridpoint) fitting + and application path. The input must still retain its month axis so that + each gridpoint has a sample distribution to map. + + Parameters + ---------- + ensemble : xr.DataArray + Generated ensemble (realization, month, lat, lon). + transform_config : VariableTransformConfig + Transform configuration providing ``fit_3d_func`` / ``apply_func``. + target_params : dict + Pre-fitted per-gridpoint target distribution parameters. + pr_baseline_field : xr.DataArray or None + Gridded first-year baseline added to convert anomalies to absolute + values before the transform (precipitation). + + Returns + ------- + xr.DataArray + Transformed ensemble with the same coords/dims as the input. + """ + data = ensemble + if pr_baseline_field is not None: + # The CMIP6 reference field carries a singleton "ens" dimension (added + # by the data getter via expand_dims). Reduce the baseline to its + # spatial (lat, lon) grid so it broadcasts cleanly over the ensemble's + # (realization, month, lat, lon) dims instead of appending a spurious + # trailing axis. + extra_dims = [d for d in pr_baseline_field.dims if d not in ("lat", "lon")] + if extra_dims: + pr_baseline_field = pr_baseline_field.isel( + {d: 0 for d in extra_dims}, drop=True + ) + # Broadcast (lat, lon) baseline over realization and month + data = data + pr_baseline_field + + # Fit Gaussian per gridpoint to the generated ensemble, then map to target. + gaussian_params = transform_config.fit_3d_func(data.values, "gaussian") + transformed = transform_config.apply_func( + data.values, + gaussian_params, + target_params, + target_dist=transform_config.transform_type, + ) + return xr.DataArray(transformed, coords=data.coords, dims=data.dims) + def _apply_impacts( self, var_output, variable, impact_configs, custom_regions=None, verbose=True ): diff --git a/src/meteor/noise_generator.py b/src/meteor/noise_generator.py index cc0d8f0..3462cc7 100644 --- a/src/meteor/noise_generator.py +++ b/src/meteor/noise_generator.py @@ -487,6 +487,7 @@ def generate_realization( random_seed=None, noise_only=False, add_base=None, + stochastic_pcs=None, ): """ Generate stochastic climate realizations. @@ -508,6 +509,12 @@ def generate_realization( Base climatology to add to each realization. If provided, the addition is done efficiently in NumPy before XArray conversion, avoiding expensive XArray operations. Must have compatible shape with the output. + stochastic_pcs : np.ndarray, optional + Pre-generated stochastic PCs from :meth:`generate_stochastic_pcs`. If + provided, these PCs are used instead of generating fresh ones, enabling + self-consistent ensemble generation shared with regional/global means. + Shape ``(n_time, n_modes)`` or ``(n_realizations, n_time, n_modes)``. + When provided, ``n_realizations`` is inferred from the array. Returns ------- @@ -518,7 +525,7 @@ def generate_realization( if not self.fitted: raise ValueError("Model must be fitted before generating realizations") - if random_seed is not None: + if random_seed is not None and stochastic_pcs is None: np.random.seed(random_seed) # Create time coordinate (shared across all realizations) @@ -569,11 +576,25 @@ def generate_realization( base_values = find_time_dim_and_cut(base_values, n_time, n_lat, n_lon) base_clim_np = base_values.reshape(n_time, n_lat, n_lon) + # Use pre-generated PCs if provided (self-consistent ensembles), else + # generate a fresh stochastic component per realization. + if stochastic_pcs is not None: + if stochastic_pcs.ndim == 2: + pcs_to_use = [stochastic_pcs] + else: + pcs_to_use = list(stochastic_pcs) + n_realizations = len(pcs_to_use) + else: + pcs_to_use = None + # Generate realizations (only stochastic component varies) realizations = [] - for _ in range(n_realizations): + for i in range(n_realizations): # Generate stochastic component (this is the only unique part per realization) - synthetic_pcs = self._generate_stochastic_pcs(X_exog, n_time) + if pcs_to_use is not None: + synthetic_pcs = pcs_to_use[i] + else: + synthetic_pcs = self._generate_stochastic_pcs(X_exog, n_time) # Reconstruct anomalies (NumPy) reconstructed_anomalies = synthetic_pcs @ self._physical_components() diff --git a/tests/unit/test_meteor_interface.py b/tests/unit/test_meteor_interface.py index a022f0d..5e117dd 100644 --- a/tests/unit/test_meteor_interface.py +++ b/tests/unit/test_meteor_interface.py @@ -13,7 +13,18 @@ import xarray as xr from meteor.ensemble_output import EnsembleOutput -from meteor.meteor_interface import MeteorInterface, _get_default_config +from meteor.meteor_interface import ( + GenerationInputs, + MeteorInterface, + PatternScalingResult, + _get_default_config, + _stack_realizations, +) +from meteor.precipitation_transform import ( + apply_distribution_transform, + fit_distribution_parameters_3d, +) +from meteor.variable_transforms import VariableTransformConfig def _set_trained(interface, variables): @@ -79,8 +90,10 @@ def test_get_default_config(): assert not generic_config["transform"] -def test_tas_converted_to_anomalies(mock_interface): - """Test that tas data is converted to anomalies from piControl baseline.""" +def test_tas_skips_reference_data_loading(mock_interface): + """tas has no distribution transform, so no CMIP6/piControl reference data + should be loaded during time series generation (the former anomaly-conversion + path only fed transform fitting, which tas does not perform).""" # Create mock piControl data with known mean picontrol_mean = 288.0 # K picontrol_data = xr.Dataset( @@ -153,7 +166,7 @@ def test_tas_converted_to_anomalies(mock_interface): np.zeros(100), dims=["month"], coords={"month": range(100)} ) - # This should trigger the anomaly conversion for tas + # This should NOT trigger any reference-data loading for tas _ = mock_interface._generate_timeseries( variable="tas", scenario="ssp245", @@ -165,16 +178,8 @@ def test_tas_converted_to_anomalies(mock_interface): verbose=False, ) - # Verify that make_meteor_training_data_composite was called for both scenario and piControl - calls = [ - call[0] - for call in mock_interface.data_getter.make_meteor_training_data_composite.call_args_list - ] - - # Should have been called with scenario data - assert any("ssp245" in str(call) or "historical" in str(call) for call in calls) - # Should have been called with piControl data (for tas only) - assert any("piControl" in str(call) for call in calls) + # tas has no transform, so no CMIP6/piControl reference data is loaded. + mock_interface.data_getter.make_meteor_training_data_composite.assert_not_called() def test_pr_not_converted_to_anomalies(mock_interface): @@ -577,7 +582,16 @@ def test_generate_gridded_climatology_no_noise_single_realization(interface_fact coords={"month": np.arange(24), "lat": [0, 1], "lon": [0, 1]}, ) interface._get_or_compute_pattern_scaling = MagicMock( - return_value=(monthly_prediction, monthly_warming) + return_value=PatternScalingResult( + monthly_prediction=monthly_prediction, + monthly_warming=monthly_warming, + em_data=None, + conc_data=None, + full_monthly_warming=monthly_warming, + base_year=2000, + start_month_idx=0, + end_month_idx=24, + ) ) gridded = interface._generate_gridded( @@ -679,3 +693,338 @@ def test_compute_timeseries_scaling(): ) assert scaling_factor.shape == (1, 1, 1) assert np.isclose(scaling_factor.values[0, 0, 0], 1.0) + +def test_interface_with_tabids_in_data_getter_kwargs(): + """MeteorInterface passes tabids from data_getter_kwargs to Cmip6MeteorDataGetter.""" + with patch("meteor.meteor_interface.Cmip6MeteorDataGetter") as mock_getter_class: + mock_getter_class.return_value = MagicMock() + + MeteorInterface( + model="TestModel", + variables=["tas"], + cache_dir="/tmp", + data_getter_kwargs={"tabids": "Amon"}, + ) + + _, call_kwargs = mock_getter_class.call_args + assert call_kwargs["tabids"] is not None + + +def test_train_with_custom_training_scenario(interface_factory): + """train() with a non-default training_scenario stores it in _training_config.""" + interface, _ = interface_factory(model="TestModel", variables=("tas",)) + + with ( + patch.object(interface, "_train_pattern_scaling"), + patch.object(interface, "_train_noise_model"), + ): + interface.train(training_scenario="ssp370", verbose=False) + + assert interface._training_config["tas"]["training_scenario"] == "ssp370" + + +def test_train_with_variable_configs(interface_factory): + """train() with variable_configs applies per-variable overrides.""" + interface, _ = interface_factory(model="TestModel", variables=("tas",)) + + with ( + patch.object(interface, "_train_pattern_scaling"), + patch.object(interface, "_train_noise_model"), + ): + interface.train(variable_configs={"tas": {"n_modes_noise": 20}}, verbose=False) + + assert interface._training_config["tas"]["n_modes_noise"] == 20 + + +def test_generate_saves_to_file(interface_factory): + """generate_ensemble_outputs() calls ensemble.to_netcdf when save_to is given.""" + interface, _ = interface_factory(model="TestModel", variables=("tas",)) + _set_trained(interface, ["tas"]) + + interface._generate_timeseries = MagicMock( + return_value={"global": xr.DataArray([290.0], dims=["time"])} + ) + + with patch.object(EnsembleOutput, "to_netcdf") as mock_save: + interface.generate_ensemble_outputs( + scenario="ssp245", + start_year=2020, + end_year=2020, + n_realizations=1, + timeseries=["global"], + save_to="/tmp/test_output.nc", + verbose=False, + ) + + mock_save.assert_called_once_with("/tmp/test_output.nc") + + +# ============================================================================= +# Shared generation helpers (introduced by the unified-generation refactor) +# ============================================================================= + + +def test_stack_realizations_single_adds_realization_dim(): + """A single realization is promoted to a length-1 ``realization`` dimension.""" + field = xr.DataArray( + np.ones((3, 2, 2)), + dims=["month", "lat", "lon"], + coords={"month": range(3), "lat": [0, 1], "lon": [0, 1]}, + ) + + stacked = _stack_realizations([field]) + + assert "realization" in stacked.dims + assert stacked.sizes["realization"] == 1 + # The underlying field is unchanged aside from the new leading axis. + assert np.array_equal(stacked.isel(realization=0).values, field.values) + + +def test_stack_realizations_multiple_concatenates(): + """Multiple realizations are concatenated along the ``realization`` dim.""" + fields = [ + xr.DataArray( + np.full((3, 2, 2), float(i)), + dims=["month", "lat", "lon"], + coords={"month": range(3), "lat": [0, 1], "lon": [0, 1]}, + ) + for i in range(4) + ] + + stacked = _stack_realizations(fields) + + assert stacked.sizes["realization"] == 4 + # Each member retains its distinct values. + for i in range(4): + assert np.all(stacked.isel(realization=i).values == float(i)) + + +def test_get_transform_config_handles_dict_and_direct(interface_factory): + """_get_transform_config resolves both fitted (dict) and unfitted configs.""" + interface, _ = interface_factory(model="TestModel", variables=("pr",)) + config = VariableTransformConfig("pr", "gamma", "positivity") + + # Fitted case: stored as a dict with a 'config' entry. + interface.transforms["pr"] = {"config": config, "target_params": {}} + assert interface._get_transform_config("pr") is config + + # Unfitted case: stored as the config object directly. + interface.transforms["pr"] = config + assert interface._get_transform_config("pr") is config + + # Missing case: variable with no transform registered. + assert interface._get_transform_config("tas") is None + + +def test_prepare_generation_slices_full_trajectory_pcs(interface_factory): + """PCs are generated once over the full trajectory then sliced to the window. + + The autoregressive spin-up transient must be parked at the trajectory start, + so PCs are generated over ``full_monthly_warming`` and only afterwards sliced + to ``[start_month_idx:end_month_idx]``. + """ + interface, _ = interface_factory(model="TestModel", variables=("pr",)) + _set_trained(interface, ["pr"]) + + n_months_full = 36 # 3 years from base_year + start_idx, end_idx = 12, 24 # output window = second year + full_warming = np.linspace(0.0, 3.0, n_months_full) + + pattern = PatternScalingResult( + monthly_prediction=xr.DataArray(np.zeros(12), dims=["month"]), + monthly_warming=full_warming[start_idx:end_idx], + em_data=None, + conc_data=None, + full_monthly_warming=full_warming, + base_year=2000, + start_month_idx=start_idx, + end_month_idx=end_idx, + ) + interface._get_or_compute_pattern_scaling = MagicMock(return_value=pattern) + + full_pcs = np.arange(3 * n_months_full * 4).reshape(3, n_months_full, 4) + interface.noise_models["pr"].generate_stochastic_pcs.return_value = full_pcs + + gen_inputs = interface._prepare_generation( + "pr", "ssp245", 2001, 2001, n_realizations=3, verbose=False + ) + + assert isinstance(gen_inputs, GenerationInputs) + # PCs generated over the FULL trajectory (length 36), not the window. + call_args, _ = interface.noise_models["pr"].generate_stochastic_pcs.call_args + assert len(call_args[0]) == n_months_full + # Returned PCs are sliced to the output window. + assert gen_inputs.stochastic_pcs.shape == (3, end_idx - start_idx, 4) + assert np.array_equal(gen_inputs.stochastic_pcs, full_pcs[:, start_idx:end_idx, :]) + + +def test_prepare_generation_normalizes_2d_pcs(interface_factory): + """A 2D (single-realization) PC array is promoted to a leading realization axis.""" + interface, _ = interface_factory(model="TestModel", variables=("pr",)) + _set_trained(interface, ["pr"]) + + pattern = PatternScalingResult( + monthly_prediction=xr.DataArray(np.zeros(12), dims=["month"]), + monthly_warming=np.zeros(12), + em_data=None, + conc_data=None, + full_monthly_warming=np.zeros(24), + base_year=2000, + start_month_idx=0, + end_month_idx=12, + ) + interface._get_or_compute_pattern_scaling = MagicMock(return_value=pattern) + interface.noise_models["pr"].generate_stochastic_pcs.return_value = np.zeros( + (24, 4) + ) + + gen_inputs = interface._prepare_generation( + "pr", "ssp245", 2000, 2000, n_realizations=1, verbose=False + ) + + assert gen_inputs.stochastic_pcs.shape == (1, 12, 4) + + +def test_prepare_generation_no_noise_skips_pcs(interface_factory): + """With include_noise=False no PCs are generated and the field is None.""" + interface, _ = interface_factory(model="TestModel", variables=("pr",)) + _set_trained(interface, ["pr"]) + + pattern = PatternScalingResult( + monthly_prediction=xr.DataArray(np.zeros(12), dims=["month"]), + monthly_warming=np.zeros(12), + em_data=None, + conc_data=None, + full_monthly_warming=np.zeros(24), + base_year=2000, + start_month_idx=0, + end_month_idx=12, + ) + interface._get_or_compute_pattern_scaling = MagicMock(return_value=pattern) + + gen_inputs = interface._prepare_generation( + "pr", "ssp245", 2000, 2000, n_realizations=5, include_noise=False, verbose=False + ) + + assert gen_inputs.stochastic_pcs is None + interface.noise_models["pr"].generate_stochastic_pcs.assert_not_called() + + +def test_generate_gridded_slice_no_noise_reduces_time(interface_factory): + """Without noise the slice is the time-mean of the base pattern, no noise calls.""" + interface, _ = interface_factory(model="TestModel", variables=("tas",)) + noise_model = MagicMock() + + monthly_prediction = xr.DataArray( + np.arange(24 * 2 * 2, dtype=float).reshape(24, 2, 2), + dims=["month", "lat", "lon"], + coords={"month": np.arange(24), "lat": [0, 1], "lon": [0, 1]}, + ) + + result = interface._generate_gridded_slice( + monthly_prediction, + np.zeros(24), + noise_model, + stochastic_pcs=None, + start_idx=0, + end_idx=12, + include_noise=False, + reduce_time=True, + ) + + noise_model.generate_realization.assert_not_called() + assert result.dims == ("realization", "lat", "lon") + assert result.sizes["realization"] == 1 + expected = monthly_prediction.isel(month=slice(0, 12)).mean(dim="month") + assert np.allclose(result.isel(realization=0).values, expected.values) + + +def test_generate_gridded_slice_passes_sliced_pcs(interface_factory): + """With noise the window-sliced PCs and base climatology are forwarded.""" + interface, _ = interface_factory(model="TestModel", variables=("tas",)) + noise_model = MagicMock() + + monthly_prediction = xr.DataArray( + np.zeros((24, 2, 2)), + dims=["month", "lat", "lon"], + coords={"month": np.arange(24), "lat": [0, 1], "lon": [0, 1]}, + ) + noise_model.generate_realization.return_value = [ + xr.DataArray( + np.zeros((12, 2, 2)), + dims=["month", "lat", "lon"], + coords={"month": np.arange(12), "lat": [0, 1], "lon": [0, 1]}, + ) + for _ in range(2) + ] + + stochastic_pcs = np.arange(2 * 24 * 4).reshape(2, 24, 4) + + result = interface._generate_gridded_slice( + monthly_prediction, + np.zeros(24), + noise_model, + stochastic_pcs=stochastic_pcs, + start_idx=12, + end_idx=24, + include_noise=True, + reduce_time=False, + ) + + _, call_kwargs = noise_model.generate_realization.call_args + # PCs are sliced to the requested window before being passed on. + assert np.array_equal(call_kwargs["stochastic_pcs"], stochastic_pcs[:, 12:24, :]) + assert call_kwargs["noise_only"] is True + # Monthly output retains the month axis and both realizations. + assert result.sizes["realization"] == 2 + assert "month" in result.dims + + +def test_apply_gridded_transform_drops_singleton_ens_dim(interface_factory): + """The (lat, lon) baseline must not append a spurious axis to the ensemble. + + The CMIP6 data getter adds a singleton ``ens`` dimension to its fields. When + that baseline is added to the generated ensemble it must be reduced to its + spatial grid first, otherwise xarray broadcasting produces a 5D array that the + per-gridpoint transform rejects. This is a regression test for that bug. + """ + interface, _ = interface_factory(model="TestModel", variables=("pr",)) + + rng = np.random.default_rng(0) + ensemble = xr.DataArray( + rng.normal(0.0, 1e-6, size=(2, 12, 2, 2)), + dims=["realization", "month", "lat", "lon"], + coords={ + "realization": [0, 1], + "month": np.arange(12), + "lat": [0, 1], + "lon": [0, 1], + }, + ) + # Baseline carries a singleton ``ens`` dim as produced by the data getter. + pr_baseline_field = xr.DataArray( + rng.uniform(1e-5, 2e-5, size=(1, 2, 2)), + dims=["ens", "lat", "lon"], + coords={"ens": [1], "lat": [0, 1], "lon": [0, 1]}, + ) + + target_ref = rng.uniform(1e-5, 3e-5, size=(24, 2, 2)) + target_params = fit_distribution_parameters_3d(target_ref, "gamma") + + config = VariableTransformConfig( + "pr", + "gamma", + "positivity", + fit_3d_func=fit_distribution_parameters_3d, + apply_func=apply_distribution_transform, + ) + + result = interface._apply_gridded_transform( + ensemble, config, target_params, pr_baseline_field + ) + + # Dimensions are preserved (no spurious ``ens`` axis) ... + assert result.dims == ("realization", "month", "lat", "lon") + assert result.shape == (2, 12, 2, 2) + # ... and the gamma transform guarantees non-negative precipitation. + assert np.all(result.values >= 0) From a5c56564fcfd8435fdb40e61a6da539d49e91c68 Mon Sep 17 00:00:00 2001 From: maritsandstad Date: Sat, 13 Jun 2026 22:49:12 +0200 Subject: [PATCH 07/11] Linting, Changelog (cherry picked from commit 402e20bc5ea89f97e1fae42e85af78df4c83c2cb) --- CHANGELOG.rst | 1 + tests/unit/test_meteor_interface.py | 1 + 2 files changed, 2 insertions(+) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 8cac9b6..46a5197 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -24,6 +24,7 @@ The changes listed in this file are categorised as follows: ### Fixed - Variable length timeseries now works also when don't have "year" as time dimension. +- Fixes to generate annual and monthly gridded ensembles with unified noise and preserving more of the variance. [Version 1.6.0] diff --git a/tests/unit/test_meteor_interface.py b/tests/unit/test_meteor_interface.py index 5e117dd..c03df01 100644 --- a/tests/unit/test_meteor_interface.py +++ b/tests/unit/test_meteor_interface.py @@ -694,6 +694,7 @@ def test_compute_timeseries_scaling(): assert scaling_factor.shape == (1, 1, 1) assert np.isclose(scaling_factor.values[0, 0, 0], 1.0) + def test_interface_with_tabids_in_data_getter_kwargs(): """MeteorInterface passes tabids from data_getter_kwargs to Cmip6MeteorDataGetter.""" with patch("meteor.meteor_interface.Cmip6MeteorDataGetter") as mock_getter_class: From 9250e37738a5095e7c1fa678cd05e734ac33f655 Mon Sep 17 00:00:00 2001 From: Ben Sanderson Date: Tue, 16 Jun 2026 12:42:35 +0200 Subject: [PATCH 08/11] Fit precipitation quantile map per month-of-year, full window MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The Gamma transform was conflating the seasonal cycle and the climate-change trend with internal variability: - Time-series path: the Gaussian was fit on the whole flattened (n_real, n_time) ensemble, so sigma was dominated by the wet/dry-season swing and the inter-realization noise band got squashed in the CDF→PPF roundtrip. - Gridded path: the transform was applied per year-slice, so the per- gridpoint Gaussian mean was the local 1-year climatology, which normalized every year to its own local distribution and wiped the climate-change trend at each gridpoint. Fit and apply per month-of-year (12 fits) so the seasonal cycle lives in the per-month means instead of sigma, and on the gridded side build a single full-window transformed ensemble up front so the annual / monthly / climatology slicers just read from it. Old pooled fit/apply functions are kept for backwards compatibility; the seasonal variants are wired through VariableTransformConfig and only used when present. Co-Authored-By: Claude Opus 4.7 (1M context) (cherry picked from commit 706a94a24b194aa2efd15ec9c88a2498d14fe56d) --- src/meteor/meteor_interface.py | 231 +++++++++++++++++++------- src/meteor/precipitation_transform.py | 193 +++++++++++++++++++++ src/meteor/variable_transforms.py | 16 ++ 3 files changed, 376 insertions(+), 64 deletions(-) diff --git a/src/meteor/meteor_interface.py b/src/meteor/meteor_interface.py index 560f8a6..f637bac 100644 --- a/src/meteor/meteor_interface.py +++ b/src/meteor/meteor_interface.py @@ -1513,23 +1513,40 @@ def _generate_timeseries( if variable == "pr" and pr_baseline_agg is not None: ensemble_for_transform = ensemble_for_transform + pr_baseline_agg - # Fit Gaussian to generated data - gaussian_params = transform_config.fit_1d_func( - ensemble_for_transform, "gaussian" + # Fit per-month-of-year (seasonal) when the transform exposes it + # — otherwise σ_gaussian is dominated by the seasonal cycle and + # the quantile map squashes the inter-realization noise band. + use_seasonal = ( + transform_config.fit_1d_seasonal_func is not None + and transform_config.apply_seasonal_func is not None ) - # Fit target distribution to CMIP6 data - target_params = transform_config.fit_1d_func( - cmip6_agg, transform_config.transform_type - ) - - # Apply transform - transformed_ensemble = transform_config.apply_func( - ensemble_for_transform, - gaussian_params, - target_params, - target_dist=transform_config.transform_type, - ) + if use_seasonal: + gaussian_params = transform_config.fit_1d_seasonal_func( + ensemble_for_transform, "gaussian" + ) + target_params = transform_config.fit_1d_seasonal_func( + cmip6_agg, transform_config.transform_type + ) + transformed_ensemble = transform_config.apply_seasonal_func( + ensemble_for_transform, + gaussian_params, + target_params, + target_dist=transform_config.transform_type, + ) + else: + gaussian_params = transform_config.fit_1d_func( + ensemble_for_transform, "gaussian" + ) + target_params = transform_config.fit_1d_func( + cmip6_agg, transform_config.transform_type + ) + transformed_ensemble = transform_config.apply_func( + ensemble_for_transform, + gaussian_params, + target_params, + target_dist=transform_config.transform_type, + ) results[agg] = transformed_ensemble else: @@ -1613,22 +1630,28 @@ def _generate_gridded( elif verbose: # pragma: no cover print(" → Using shared stochastic PC realizations (gridded)") - # Resolve transform and (only for transform variables, e.g. pr) load the - # CMIP6 reference field and fit the per-gridpoint target distribution once. + # Resolve transform. For variables with a distribution transform (e.g. + # pr) we generate the full prediction window in one shot and apply the + # seasonal (per-month-of-year, per-gridpoint) transform once on the + # whole window. Doing it per year-slice (the old path) re-fitted the + # Gaussian on just 12 months at each gridpoint, which normalised every + # year to its own local mean and wiped the climate-change trend. transform_config = self._get_transform_config(variable) + pre_transformed_ensemble = None target_params = None pr_baseline_field = None if transform_config and transform_config.transform_type: - ssp_data, pr_baseline_field = self._load_transform_reference( - variable, start_year, end_year, verbose=verbose - ) - if verbose: # pragma: no cover - print( - f" → Fitting per-gridpoint {transform_config.transform_type} " - "target distribution..." - ) - target_params = transform_config.fit_3d_func( - ssp_data, transform_config.transform_type + pre_transformed_ensemble = self._build_full_window_transformed_ensemble( + variable, + monthly_prediction, + monthly_warming, + noise_model, + stochastic_pcs, + start_year, + end_year, + transform_config, + include_noise, + verbose=verbose, ) # Extract requested time slices @@ -1642,6 +1665,27 @@ def _generate_gridded( def year_to_month_idx(year): return (year - start_year) * 12 + def _get_ensemble_slice(s_idx, e_idx, reduce_time): + """Slice the pre-transformed ensemble, or generate+transform per slice.""" + if pre_transformed_ensemble is not None: + sliced = pre_transformed_ensemble.isel(month=slice(s_idx, e_idx)) + if reduce_time: + sliced = sliced.mean(dim="month") + return sliced + return self._generate_gridded_slice( + monthly_prediction, + monthly_warming, + noise_model, + stochastic_pcs, + s_idx, + e_idx, + include_noise, + reduce_time=reduce_time, + transform_config=transform_config, + target_params=target_params, + pr_baseline_field=pr_baseline_field, + ) + # Annual means (12-month average of each year) if "annual" in gridded_spec: if verbose: # pragma: no cover @@ -1653,18 +1697,8 @@ def year_to_month_idx(year): start_idx = year_to_month_idx(year) end_idx = start_idx + 12 if start_idx >= 0 and end_idx <= n_months: - annual_fields[year] = self._generate_gridded_slice( - monthly_prediction, - monthly_warming, - noise_model, - stochastic_pcs, - start_idx, - end_idx, - include_noise, - reduce_time=True, - transform_config=transform_config, - target_params=target_params, - pr_baseline_field=pr_baseline_field, + annual_fields[year] = _get_ensemble_slice( + start_idx, end_idx, reduce_time=True ) else: if verbose: # pragma: no cover @@ -1684,18 +1718,8 @@ def year_to_month_idx(year): start_idx = year_to_month_idx(year) end_idx = start_idx + 12 if start_idx >= 0 and end_idx <= n_months: - monthly_fields[year] = self._generate_gridded_slice( - monthly_prediction, - monthly_warming, - noise_model, - stochastic_pcs, - start_idx, - end_idx, - include_noise, - reduce_time=False, - transform_config=transform_config, - target_params=target_params, - pr_baseline_field=pr_baseline_field, + monthly_fields[year] = _get_ensemble_slice( + start_idx, end_idx, reduce_time=False ) else: if verbose: # pragma: no cover @@ -1718,19 +1742,7 @@ def year_to_month_idx(year): end_idx = year_to_month_idx(clim_end + 1) # +1 to include end year if start_idx >= 0 and end_idx <= n_months: climatology_fields[f"{clim_start}-{clim_end}"] = ( - self._generate_gridded_slice( - monthly_prediction, - monthly_warming, - noise_model, - stochastic_pcs, - start_idx, - end_idx, - include_noise, - reduce_time=True, - transform_config=transform_config, - target_params=target_params, - pr_baseline_field=pr_baseline_field, - ) + _get_ensemble_slice(start_idx, end_idx, reduce_time=True) ) else: if verbose: # pragma: no cover @@ -1883,6 +1895,97 @@ def _apply_gridded_transform( ) return xr.DataArray(transformed, coords=data.coords, dims=data.dims) + def _build_full_window_transformed_ensemble( + self, + variable, + monthly_prediction, + monthly_warming, + noise_model, + stochastic_pcs, + start_year, + end_year, + transform_config, + include_noise, + verbose=True, + ): + """ + Generate the gridded ensemble over the FULL prediction window and apply + the seasonal (per-month-of-year, per-gridpoint) distribution transform + once. + + Fitting the Gaussian half of the quantile map over the full window — and + per month-of-year rather than across all months at once — keeps the + long-term trend in the variance the map carries through (so the gridded + trend is preserved) while removing the seasonal cycle from σ (so the + inter-realization noise band is preserved at every gridpoint). + + Returns + ------- + xr.DataArray + Transformed monthly ensemble (realization, month, lat, lon) spanning + ``start_year..end_year`` inclusive. Callers slice this once per + requested annual / monthly / climatology field. + """ + if verbose: # pragma: no cover + print( + " → Generating full-window gridded ensemble for " + f"{transform_config.transform_type} transform" + ) + + ssp_data, pr_baseline_field = self._load_transform_reference( + variable, start_year, end_year, verbose=verbose + ) + + if verbose: # pragma: no cover + print( + f" → Fitting per-month-of-year per-gridpoint " + f"{transform_config.transform_type} target distribution..." + ) + target_params = transform_config.fit_3d_seasonal_func( + ssp_data, transform_config.transform_type + ) + + if include_noise: + realizations = noise_model.generate_realization( + monthly_warming, + noise_only=True, + add_base=monthly_prediction, + stochastic_pcs=stochastic_pcs, + ) + if not isinstance(realizations, list): + realizations = [realizations] + else: + realizations = [monthly_prediction] + ensemble = _stack_realizations(realizations) + + data = ensemble + if pr_baseline_field is not None: + extra_dims = [ + d for d in pr_baseline_field.dims if d not in ("lat", "lon") + ] + if extra_dims: + pr_baseline_field = pr_baseline_field.isel( + {d: 0 for d in extra_dims}, drop=True + ) + data = data + pr_baseline_field + + if verbose: # pragma: no cover + print( + " → Fitting per-month-of-year per-gridpoint Gaussian on " + "full-window generated ensemble..." + ) + gaussian_params = transform_config.fit_3d_seasonal_func( + data.values, "gaussian" + ) + + transformed = transform_config.apply_seasonal_func( + data.values, + gaussian_params, + target_params, + target_dist=transform_config.transform_type, + ) + return xr.DataArray(transformed, coords=data.coords, dims=data.dims) + def _apply_impacts( self, var_output, variable, impact_configs, custom_regions=None, verbose=True ): diff --git a/src/meteor/precipitation_transform.py b/src/meteor/precipitation_transform.py index e0c9eb5..e56f149 100644 --- a/src/meteor/precipitation_transform.py +++ b/src/meteor/precipitation_transform.py @@ -368,6 +368,199 @@ def apply_distribution_transform( return transformed +# ============================================================================= +# Seasonal (per-month-of-year) variants +# ============================================================================= +# +# Fitting/applying one quantile map across all months conflates the seasonal +# cycle with internal variability: the Gaussian std becomes dominated by the +# seasonal swing (very large for variables like precipitation), so the CDF +# squashes inter-realization noise into a narrow quantile band and the gamma +# PPF then maps it to a narrow output band. Fitting per month-of-year removes +# the seasonal contribution from the variance the quantile map sees, so the +# transform actually preserves the within-month internal variability. + + +def _month_of_year_indices(n_time): + """Return a list of 12 index arrays selecting each month-of-year from a + contiguous monthly time axis of length ``n_time``. Assumes the series + starts in January; partial trailing years are fine (the last month-of-year + bins will just have one fewer sample).""" + return [np.arange(m, n_time, 12) for m in range(12)] + + +def fit_distribution_parameters_1d_seasonal(timeseries_data, distribution="gamma"): + """ + Fit a separate distribution per month-of-year (12 fits) to a 1D-time series. + + Parameters + ---------- + timeseries_data : np.ndarray or xr.DataArray + Monthly data. Last axis is time and must be a multiple of 12. Earlier + axes (realization / ensemble) are pooled into each per-month fit. + distribution : str + Distribution to fit at each month-of-year. + + Returns + ------- + dict + Each parameter key maps to a 1D array of length 12 (Jan..Dec). + """ + if isinstance(timeseries_data, xr.DataArray): + data = timeseries_data.values + else: + data = timeseries_data + + n_time = data.shape[-1] + month_idx = _month_of_year_indices(n_time) + + per_month = [ + fit_distribution_parameters_1d(data[..., idx], distribution=distribution) + for idx in month_idx + ] + + keys = per_month[0].keys() + return {k: np.array([p[k] for p in per_month]) for k in keys} + + +def fit_distribution_parameters_3d_seasonal(spatial_data, distribution="gamma"): + """ + Fit per-gridpoint distributions separately for each month-of-year. + + Parameters + ---------- + spatial_data : np.ndarray or xr.DataArray + Shape (n_time, n_lat, n_lon) or (n_ensemble, n_time, n_lat, n_lon). + ``n_time`` must be a multiple of 12. + distribution : str + 'gaussian' or 'gamma'. + + Returns + ------- + dict + Each parameter key maps to an array of shape (12, n_lat, n_lon). + """ + if isinstance(spatial_data, xr.DataArray): + data = spatial_data.values + else: + data = spatial_data + + if data.ndim == 3: + time_axis = 0 + n_time = data.shape[0] + elif data.ndim == 4: + time_axis = 1 + n_time = data.shape[1] + else: + raise ValueError( + "Data must be 3D (n_time, n_lat, n_lon) or " + f"4D (n_ensemble, n_time, n_lat, n_lon), got shape {data.shape}" + ) + + month_idx = _month_of_year_indices(n_time) + + per_month = [ + fit_distribution_parameters_3d( + np.take(data, idx, axis=time_axis), distribution=distribution + ) + for idx in month_idx + ] + + keys = per_month[0].keys() + return {k: np.stack([p[k] for p in per_month], axis=0) for k in keys} + + +def apply_distribution_transform_seasonal( + gaussian_data, gaussian_params, target_params, target_dist="gamma" +): + """ + Apply a per-month-of-year quantile transform. + + Automatically detects 1D vs 3D based on the rank of the parameter arrays + (``ndim == 1`` -> 1D scalar-per-month; ``ndim == 3`` -> 3D per-gridpoint). + + Parameters + ---------- + gaussian_data : np.ndarray or xr.DataArray + Shape (n_realizations, n_time) for 1D or + (n_realizations, n_time, n_lat, n_lon) / (n_time, n_lat, n_lon) for 3D. + ``n_time`` must be a multiple of 12. + gaussian_params, target_params : dict + Per-month-of-year parameter arrays. For 1D: shape (12,) each. + For 3D: shape (12, n_lat, n_lon) each. + target_dist : str + Target distribution name. + + Returns + ------- + Same type/shape as input. + """ + is_xarray = isinstance(gaussian_data, xr.DataArray) + if is_xarray: + coords = gaussian_data.coords + dims = gaussian_data.dims + data = gaussian_data.values + else: + data = gaussian_data + + sample_param = gaussian_params["mean"] + if sample_param.ndim == 1: + # 1D scalar-per-month + if data.ndim != 2: + raise ValueError( + "For 1D seasonal transform, expected 2D data (n_real, n_time), " + f"got shape {data.shape}" + ) + time_axis = 1 + elif sample_param.ndim == 3: + # 3D per-gridpoint-per-month: params shape (12, n_lat, n_lon) + if data.ndim == 3: + time_axis = 0 + elif data.ndim == 4: + time_axis = 1 + else: + raise ValueError( + "For 3D seasonal transform, expected 3D or 4D data, " + f"got shape {data.shape}" + ) + else: + raise ValueError( + "Seasonal params must have ndim 1 (scalar per month) or 3 " + f"(per gridpoint per month), got ndim={sample_param.ndim}" + ) + + n_time = data.shape[time_axis] + month_idx = _month_of_year_indices(n_time) + + transformed = np.empty_like(data, dtype=np.float64) + for m, idx in enumerate(month_idx): + data_m = np.take(data, idx, axis=time_axis) + + if sample_param.ndim == 1: + # Extract Python scalars so the underlying apply_distribution_transform + # takes the 1D-scalar path (which uses np.isscalar). + gp_m = {k: float(v[m]) for k, v in gaussian_params.items()} + tp_m = {k: float(v[m]) for k, v in target_params.items()} + else: + gp_m = {k: v[m] for k, v in gaussian_params.items()} + tp_m = {k: v[m] for k, v in target_params.items()} + + transformed_m = apply_distribution_transform( + data_m, gp_m, tp_m, target_dist=target_dist + ) + + # Write back into the appropriate slice of ``transformed``. + if time_axis == 0: + transformed[idx] = transformed_m + else: + # time_axis == 1 covers both 2D and 4D layouts. + transformed[:, idx] = transformed_m + + if is_xarray: + return xr.DataArray(transformed, coords=coords, dims=dims) + return transformed + + # ============================================================================= # Empirical Quantile Mapping (non-parametric alternative) # ============================================================================= diff --git a/src/meteor/variable_transforms.py b/src/meteor/variable_transforms.py index 3620a14..f543870 100644 --- a/src/meteor/variable_transforms.py +++ b/src/meteor/variable_transforms.py @@ -7,8 +7,11 @@ from meteor.precipitation_transform import ( apply_distribution_transform, + apply_distribution_transform_seasonal, fit_distribution_parameters_1d, + fit_distribution_parameters_1d_seasonal, fit_distribution_parameters_3d, + fit_distribution_parameters_3d_seasonal, ) @@ -30,6 +33,10 @@ class VariableTransformConfig: # pylint: disable=too-few-public-methods Function to fit 3D transform parameters apply_func : callable, optional Function to apply the transform + fit_1d_seasonal_func, fit_3d_seasonal_func, apply_seasonal_func : callable, optional + Per-month-of-year variants. When present, the generator path uses these + in place of the pooled variants so the quantile map sees variability, + not the seasonal cycle. """ def __init__( @@ -40,6 +47,9 @@ def __init__( fit_1d_func=None, fit_3d_func=None, apply_func=None, + fit_1d_seasonal_func=None, + fit_3d_seasonal_func=None, + apply_seasonal_func=None, ): self.name = name self.transform_type = transform_type @@ -47,6 +57,9 @@ def __init__( self.fit_1d_func = fit_1d_func self.fit_3d_func = fit_3d_func self.apply_func = apply_func + self.fit_1d_seasonal_func = fit_1d_seasonal_func + self.fit_3d_seasonal_func = fit_3d_seasonal_func + self.apply_seasonal_func = apply_seasonal_func def __repr__(self): """Return string representation of VariableTransformConfig.""" @@ -64,6 +77,9 @@ def __repr__(self): fit_1d_func=fit_distribution_parameters_1d, fit_3d_func=fit_distribution_parameters_3d, apply_func=apply_distribution_transform, + fit_1d_seasonal_func=fit_distribution_parameters_1d_seasonal, + fit_3d_seasonal_func=fit_distribution_parameters_3d_seasonal, + apply_seasonal_func=apply_distribution_transform_seasonal, ), "tas": VariableTransformConfig( name=None, From 8007a211a2ca4cd58edd37c425af4254a0922847 Mon Sep 17 00:00:00 2001 From: maritsandstad Date: Mon, 29 Jun 2026 09:15:45 +0200 Subject: [PATCH 09/11] linting (cherry picked from commit 2ad4f85, dropping the GCAM_predict.ipynb re-execution that was bundled with the lint fix.) --- src/meteor/meteor_interface.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/meteor/meteor_interface.py b/src/meteor/meteor_interface.py index f637bac..3da2183 100644 --- a/src/meteor/meteor_interface.py +++ b/src/meteor/meteor_interface.py @@ -1960,9 +1960,7 @@ def _build_full_window_transformed_ensemble( data = ensemble if pr_baseline_field is not None: - extra_dims = [ - d for d in pr_baseline_field.dims if d not in ("lat", "lon") - ] + extra_dims = [d for d in pr_baseline_field.dims if d not in ("lat", "lon")] if extra_dims: pr_baseline_field = pr_baseline_field.isel( {d: 0 for d in extra_dims}, drop=True @@ -1974,9 +1972,7 @@ def _build_full_window_transformed_ensemble( " → Fitting per-month-of-year per-gridpoint Gaussian on " "full-window generated ensemble..." ) - gaussian_params = transform_config.fit_3d_seasonal_func( - data.values, "gaussian" - ) + gaussian_params = transform_config.fit_3d_seasonal_func(data.values, "gaussian") transformed = transform_config.apply_seasonal_func( data.values, From 7c8e83a2e796f1ddf5751f2bd2cfe80cb631d1aa Mon Sep 17 00:00:00 2001 From: maritsandstad Date: Mon, 29 Jun 2026 15:33:53 +0200 Subject: [PATCH 10/11] Fix offending docstring (cherry picked from commit a11f2693648a6eff7d30b80e60d7741e44eab1cf) --- src/meteor/precipitation_transform.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/src/meteor/precipitation_transform.py b/src/meteor/precipitation_transform.py index e56f149..f9a68f6 100644 --- a/src/meteor/precipitation_transform.py +++ b/src/meteor/precipitation_transform.py @@ -382,10 +382,24 @@ def apply_distribution_transform( def _month_of_year_indices(n_time): - """Return a list of 12 index arrays selecting each month-of-year from a + """ + Return a list of 12 index arrays selecting each month-of-year + + Return a list of 12 index arrays selecting each month-of-year from a contiguous monthly time axis of length ``n_time``. Assumes the series starts in January; partial trailing years are fine (the last month-of-year - bins will just have one fewer sample).""" + bins will just have one fewer sample). + + Parameters + ---------- + n_time : int + Length of the time axis (must be a multiple of 12). + + Returns + ------- + list of np.ndarray + Each element is a 1D array of indices selecting the corresponding month-of-year. + """ return [np.arange(m, n_time, 12) for m in range(12)] From 4f22d611d47db76d4f9f519f3d13b2bf2c2c0f5d Mon Sep 17 00:00:00 2001 From: Ben Sanderson Date: Tue, 30 Jun 2026 21:10:58 +0200 Subject: [PATCH 11/11] Address review: drop debug print, flesh out module docstring - Remove the leftover `print(n_lat, n_lon, n_time)` debug call in noise_generator._generate_realizations_unified_implementation. - Expand the one-line geo_data_utils module docstring into a grouped inventory of the helpers it exposes. --- src/meteor/geo_data_utils.py | 19 ++++++++++++++++++- src/meteor/noise_generator.py | 1 - 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/src/meteor/geo_data_utils.py b/src/meteor/geo_data_utils.py index 743c849..5d0eb7a 100644 --- a/src/meteor/geo_data_utils.py +++ b/src/meteor/geo_data_utils.py @@ -1,4 +1,21 @@ -"""Utility functions for handling geographic data.""" +""" +Geographic data utilities for METEOR's gridded climate fields. + +This module groups the helpers METEOR uses to navigate, weight, and +spatially aggregate xarray datasets on (lat, lon, time) grids: + +- Coordinate discovery: :func:`get_time_name`, :func:`get_lat_name`, + :func:`get_lon_name`, :func:`get_year_series`. +- Area weighting: :func:`get_weights_for_ds`, + :func:`apply_weights_and_do_spatial_mean`. +- Spatial aggregation: :func:`global_mean`, :func:`regional_mean` (AR6 + reference regions via ``regionmask``), :func:`extract_point`, + :func:`create_region_mask` (custom bounding boxes or precomputed masks), + :func:`list_ar6_regions`. +- Time-axis helpers used by the pattern-scaling pipeline: + :func:`extend_temeperature_anomaly_timeseries_for_scaling`, + :func:`find_time_dim_and_cut`. +""" import logging diff --git a/src/meteor/noise_generator.py b/src/meteor/noise_generator.py index 3462cc7..6b9cfe8 100644 --- a/src/meteor/noise_generator.py +++ b/src/meteor/noise_generator.py @@ -572,7 +572,6 @@ def generate_realization( # Now reshape to (n_time, n_lat, n_lon) # If the time dimension doesn't match, select the first n_time steps - print(n_lat, n_lon, n_time) base_values = find_time_dim_and_cut(base_values, n_time, n_lat, n_lon) base_clim_np = base_values.reshape(n_time, n_lat, n_lon)