Conversation
# Conflicts: # tests/unit/test_config.py
resolve conflicts in config, snakefile and schema (additive: scoremaps + spectra side by side). drop the eckit_grids input from the spectra rules, mirroring 3e7660b which deleted data_download_eckit_geo_grids.
jonasbhend
left a comment
jonasbhend
left a comment
There was a problem hiding this comment.
@frazane Thanks for pushing this. I have a few comments in particular around organization and generalization that should be easy to address. Good to go after that.
| variables: List[str] = Field( | ||
| default=["T_2M", "WIND_KE", "TOT_PREC"], | ||
| description="Spectra variables. Supported: T_2M, WIND_KE (from U/V_10M), TOT_PREC.", |
There was a problem hiding this comment.
This now uses explicit aggregations, so this should probably be either TOT_PREC1 and/or TOT_PREC6
| raise ValueError( | ||
| "`variables` must list at least one variable when spectra is configured." | ||
| ) | ||
| allowed = {"T_2M", "WIND_KE", "TOT_PREC"} |
There was a problem hiding this comment.
same as above for TOT_PREC explicit aggregations
There was a problem hiding this comment.
Also, spectra for all the parameters that can be accessed with load_forecast/truth_data should be able to be computed, or does the computation require variable-specific configuration?
| def compute_source_spectra(ds: xr.Dataset, variables, lead_times, method, label): | ||
| """Compute spectra for one source (already-loaded native dataset). | ||
|
|
||
| Returns an x.Dataset with dims (variable, leadtime, wavenumber) and a shared |
There was a problem hiding this comment.
| Returns an x.Dataset with dims (variable, leadtime, wavenumber) and a shared | |
| Returns an xr.Dataset with dims (variable, leadtime, wavenumber) and a shared |
| def aggregate_spectra(spectra_files) -> xr.Dataset: | ||
| """Average power over init times (nanmean). All inputs share one grid.""" | ||
| datasets = [xr.open_dataset(f) for f in spectra_files] | ||
| try: | ||
| stacked = xr.concat(datasets, dim="init") | ||
| agg = stacked.mean(dim="init", skipna=True) | ||
| agg.attrs = datasets[0].attrs | ||
| agg["wavelength"] = datasets[0]["wavelength"].copy() | ||
| agg.load() | ||
| finally: | ||
| for ds in datasets: | ||
| ds.close() | ||
| return agg |
There was a problem hiding this comment.
In verification_aggregation.py we verify stratified by season (if configured). Could this be leveraged here, i.e. would it be easy to harmonize how we aggregate by initialization results to average results across the different sets of scores and metrics?
| def plot_experiment_spectra( | ||
| truth_file, participant_files, out_dir, variables, lead_times | ||
| ): |
There was a problem hiding this comment.
Not sure I agree on the way this module is organized. To me plotting definitely doesn't belong in compute, but would be better put in a separate sub-module. But as evalml is quite chaotic in this regard, I would be happy to defer this to a more general clean-up during the planned refactor.
There was a problem hiding this comment.
seeing that plot_power_spectra is defined in core, could this also go there?
There was a problem hiding this comment.
I suggest to move all plotting to a separate sub-module for clarity.
| VARIABLE_COMPONENTS: dict[str, tuple[list[str], float]] = { | ||
| "T_2M": (["T_2M"], 1.0), | ||
| "WIND_KE": (["U_10M", "V_10M"], 0.5), | ||
| "TOT_PREC": (["TOT_PREC"], 1.0), | ||
| } |
There was a problem hiding this comment.
This could be generalized to all available parameters (if computed directly) and only WIND_KE needs special treatment.
There was a problem hiding this comment.
Given that the regridding will be useful for other metrics, would it make sense to move into a separate module?
| def test_aggregate_spectra_nanmean_ignores_missing(tmp_path): | ||
| wl = np.array([100.0, 50.0]) | ||
|
|
||
| def mk(power, init): | ||
| ds = xr.Dataset( | ||
| {"power": (("variable", "leadtime", "wavenumber"), np.array([[power]]))}, | ||
| coords={ | ||
| "variable": ["T_2M"], | ||
| "leadtime": [6], | ||
| "wavenumber": np.arange(2), | ||
| "wavelength": ("wavenumber", wl), | ||
| }, | ||
| attrs={"dx_km": 1.1, "npoints": 10, "label": "m"}, | ||
| ) | ||
| p = tmp_path / f"s_{init}.nc" | ||
| ds.to_netcdf(p) | ||
| return p | ||
|
|
||
| a = mk([np.nan, 2.0], 0) | ||
| b = mk([4.0, 4.0], 1) | ||
| agg = compute.aggregate_spectra([a, b]) | ||
| np.testing.assert_allclose( | ||
| agg["power"].sel(variable="T_2M", leadtime=6).values, [4.0, 3.0] | ||
| ) |
There was a problem hiding this comment.
how could missing values, that we would want to exclude from aggregation, occur? Aren't all missing values real missing values that should result in a missing value in the aggregation? I am worried that we may get aggregated results based on very few non-missing samples without noticing.
| def spectra_participants(): | ||
| """participant key -> aggregated spectra path (runs + baselines).""" | ||
| out = {} | ||
| for base in BASELINE_CONFIGS: | ||
| out[base] = OUT_ROOT / f"data/baselines/{base}/spectra_aggregated.nc" | ||
| for run_id, cfg in RUN_CONFIGS.items(): | ||
| if cfg.get("_is_candidate", False): | ||
| out[run_id] = OUT_ROOT / f"data/runs/{run_id}/spectra_aggregated.nc" | ||
| return out |
There was a problem hiding this comment.
are spectra participants different from experiment participants? Or should we consolidate and use only one?
2 participants
Adds an optional power-spectra diagnostic to the experiment pipeline, so we can compare the effective resolution of model outputs against the truth and spot over-smoothing or spurious small-scale noise that point metrics like RMSE do not reveal.
The diagnostic computes 2D variance spectra (DCT by default, FFT optional) per participant and lead time, averages over initialisation times, and overlays each model against the truth with grid-resolution, effective-resolution and Nyquist reference lines.
What changed
Notes
Example output, T_2M at +24h: