Taylor Diagram with observations and multi-case - #435
Conversation
|
I should also mention that I have a test configuration file here: The web page generation does not work with this configuration... part of the reasoning for making it a draft PR for now. |
Resolved conflicts in adf_dataset.py (kept apply_scaling kwargs on top of main's hist_str-aware file lookup; factored the per-case regrid location into a helper) and create_climo_files.py (main already carries the spawn-context multiprocessing fix).
- regrid_and_vert_interp: source_grid was undefined for unstructured input (NameError); load_climo_da(...)['PS'] indexed a DataArray; the PS fallback passed a DataArray where a grid Dataset was expected; cam_regrid_loc as a plain string was indexed character-wise; weights dir needed parents=True; dropped two no-op 'final_ds = final_ds' blocks. - _create_clean_grid: derive cell edges from centers when the file has no bounds, so conservative regridding works on CAM climo files. - cam_taylor_diagram: reference data was only recognized when the label was literally 'Obs', so model-vs-model looked for the baseline in the test-case regrid files; route all reference loads through _is_ref/_load_ref_da, which falls back to the baseline climo when there is no regridded file. Obs weight directory was built from a case *name* rather than a path. - adf_config: collapse the two near-identical recursion helpers into one, and cover nested-list keyword expansion in the unit tests. - Remove leftover debug prints (create_climo_files, adf_info) and the regrid_and_vert_interp.py.OLD copy; repair the mangled comment at the end of hist_str_to_list.
vert_remap uses np.interp, which requires increasing x. Pressure-level data often arrives surface-first (descending), which silently produced a constant field equal to the topmost level.
load_regrid_da returns a DataArray, whose 'in' operator tests values, not variable names, so 'hyam' not in dset was always True and the hybrid branch could never run. Regridded files are already on pressure levels, so the pressure path is the correct one; collapse the two identical branches and drop the now-unused geocat/warnings imports. Also pass the Dataset (not the DataArray) to _determine_vertical_coord_type so its hyam fallback can fire.
load_climo_ds and load_reference_climo_ds indexed the result of load_dataset without checking for None, and load_reference_climo_ds looked up variable_defaults[var] unguarded (KeyError for a variable with no defaults entry). Also require new_unit to exist before treating a units match as evidence the data were already transformed.
…posed
The rewritten regridder only wrote the test-case files, dropping the
{base}_{var}_baseline.nc products the old driver produced for every variable
in non-obs mode. zonal_mean, global_latlon_map and polar_map read them via
load_reference_regrid_da; meridional_mean and global_latlon_vect_map glob them
directly. All five degrade with a warning rather than crashing, so a
model-vs-model run silently lost the reference and difference panels from
every plot.
- Add _write_reference_files, run once before the test-case loop, reusing
_determine_vertical_coord_type / _handle_vertical_interpolation /
_apply_ocean_mask / save_to_nc. No horizontal regrid: the reference defines
the target grid. Gated on 'not compare_obs', as the old driver was.
- Revert cam_regrid_loc to a single directory. The per-case list broke
meridional_mean and global_latlon_vect_map (Path() of a list) and
_regrid_loc_for_case (the baseline is not in case_names); regridded
filenames already carry the case name.
- _find_surface_pressure: prefer the PS carried inside the climo file, as the
old driver did, before a standalone PS climo. PS is often not in
diag_var_list, which otherwise skipped every 3-D variable.
Bugs the run surfaced, all on paths shared with the test case:
- attrs['transformed'] = True could not be written; netCDF4 rejects a Python
bool. Now 1.
- DEFAULT_PLEVS_Pa was a list; geocat's interp_hybrid_to_pressure needs .size
and utils.vert_remap needs .shape, so every hybrid- and height-coordinate
variable crashed.
- _create_clean_grid clipped the target grid CENTERS to +/-89.999999. Those
become the output coordinate, so regridded files disagreed with the
reference at the poles and xarray alignment silently dropped both polar rows
from every difference field. Only the bounds are clamped now, to +/-90.
- load_climo_ds dropped attrs (xarray arithmetic does not keep them), so every
regridded file lost 'units' and polar_map raised KeyError on it.
- cam_taylor_diagram advertised 'PRECT or PRECC+PRECL' but raised KeyError for
either when PRECT was absent. Added _load_field/_load_field_ds, which return
None for a variable that is not part of the run, and routed the loaders
through them.
Verified by a full model-vs-model run (b.e23_alpha17f.BLT1850.ne30_t232.098 vs
.093, yrs 10-11): completes successfully, 16/16 baseline files written, 3-D
references on the same 21 pressure levels as the test case, difference fields
keep all 192 latitudes, and every plotting script produces its reference and
difference panels.
The ERA5 land-sea mask this PR adds as the LANDFRAC observation is a single-timestamp file, not the 12-month climatology (time = 1..12) every other file in ADF_obs uses. load_da's squeeze() collapses its length-1 time dimension to a scalar coordinate, so data.time.dt.days_in_month succeeds but the subsequent .sel(time=...) raises KeyError: no index found for coordinate 'time', which aborted global_latlon_map for every model-vs-obs run. A time-invariant field has the same value in every season, so return it unchanged (dropping the stale scalar time coordinate, matching what the averaging path produces). Guarded on 'time' being a coordinate but not a dimension, so the existing month-dimension and bare-numpy paths are untouched.
The Taylor diagram numbered its points with a counter that advanced only for points actually drawn, while the figure legend numbers the full variable list. Any variable without data therefore shifted every later label onto the wrong name. In the model-vs-obs run just completed, three of ten variables have no observational counterpart, so four of the seven plotted points were labelled with the wrong variable -- the plot read as if TropicalLandPrecip were present when point 4 was actually Land2mTemperature. Number by position in the frame instead, which leaves gaps for missing variables and keeps labels correct; it also removes the old requirement that all cases share the same row order. Dropped the 'annos' list, which was accumulated and never used. cam_taylor_diagram carried its own copy of _create_clean_grid that still clipped the grid CENTERS to +/-89.999999 and never inferred bounds -- so the regridding behind its derived variables ran without bounds (hence the 'Latitude is outside of [-90, 90]' warnings) and would reintroduce the polar misalignment fixed in the regridder. Hoisted one implementation into adf_utils.create_clean_grid, used by both, and deleted both copies.
The Taylor diagram's EquatorialPacificStress entry could never resolve against observations because TAUX had no obs_file. The TAUX climatologies already existed in ADF_obs (built from ERA5 iews) but were referenced only by the era5-1deg defaults; the default file referenced neither TAUX nor TAUY, and TAUY had no observation at all. TAUY has now been built from the sibling ERA5 field (inss, monthly means from RDA ds633.1, 1979-2022) at 0.25 degree native and conservatively regridded to 1 degree, by AMP_toolbox/make_era5_tau_obs.py. Point both defaults files at the matching resolution: the default file at the native climatologies, as its other ERA5 entries are, and the era5-1deg file at the 1 degree pair. scale_factor: -1 is deliberately left untouched. ERA5 stores stress positive eastward/northward on the surface, CAM uses the opposite sign, and ADF applies that factor to the *model* when its climo is loaded for regridding while observations are read unscaled -- verified over 5S-5N: ERA5 -0.0228 N/m2, CAM raw +0.0283, CAM regridded -0.0283. Adding obs_scale_factor would double-flip.
The Taylor diagram accepts 'PRECT or PRECC+PRECL' for its two tropical precipitation entries. That is right against a baseline simulation, but PRECC and PRECL have no observational counterpart, so under compare_obs the check passed and then both entries silently vanished from the diagram -- which is how the obs run came out at 8 of 10 variables with no indication why. Warn explicitly in that case, and note that ADF derives PRECT from PRECC + PRECL when it builds the time series, so adding it costs the user nothing. Verified: with PRECT in diag_var_list the derivation runs, the derived field is exactly PRECC + PRECL, and it regrids onto ERAI's grid in mm/day matching the observations. The obs Taylor diagram now shows all 10 variables.
|
This comment authored by Claude, and so is very verbose. But it is comprehensive about what was changed in this latest update. At this point, the code seems to run and is improved. I might do another round of refactoring, but at this time I think it is ready for review/testing. Branch The blocker: the rewritten regridder stopped producing the referenceThe old driver wrote two products per variable; the rewrite wrote only one.
The baseline file is not limited to 3-D fields — in the old driver a 2-D Fixed by adding Bugs the end-to-end run exposed — all on paths shared with the test caseThese were latent in the PR; the reference pass just hit them first.
Other defects fixed while reviewingRegridder:
Housekeeping: removed leftover Verified by the run
Worth deciding before merge
Model-vs-obs run — the PR's headline featureAlso run end to end ( Three further defects it exposed, all fixed:
Two obs-mode limitations that are not bugs but constrain the feature:
Not exercisedUnstructured ( Memory note: regridding a 3-D CAM field to ERA5 0.25° needs ~3.2 GB per |
Ran the Taylor diagram against observations with two unrelated cases on
different vertical grids -- ceresmip_amip02 (cam.h0, 32 levels, 1995-2004) and
b.e30_alpha09d_m.B1850C_MTso_Gris_Marbl.ne30_t233_wgx3.377 (cam.h0a, 93 levels,
0138-0147). That combination had never been run, and it found three problems.
Surface pressure was being read in hPa and used as if it were Pa. This PR
inverted main's precedence for PS: main prefers the copy carried inside the
variable's own climo file, which is in Pa, and only falls back to the regridded
file; this branch checked "{target}_{case}_PS_regridded.nc" first. That file has
the variable defaults applied, so it is in hPa. Handed to
interp_hybrid_to_pressure with P0=100000 it squeezes the whole model column into
3.6-178 hPa, and every requested level outside that band comes out NaN -- the
entire troposphere, silently. U300 disappeared from the diagram and the two
column-integral entries still produced plausible-looking numbers computed from
stratosphere-only columns. Normalize PS to Pa at both consumption sites instead,
falling back to a magnitude test when the units attribute is missing.
Every earlier verification missed this because the example model-vs-model config
has no PS in diag_var_list, so no regridded PS file existed and the code fell
through to the correct path. Adding PS is what triggers it.
Multi-case website generation died with FileNotFoundError at the last step.
_write_run_info_to_web wrote the run-info markdown only to
plot_location[0]/website, but create_website reads it back from the directory
belonging to each web_data's case -- including the synthetic "multi-case" entry,
which points at main_website. Write it into every website directory the run will
build. Single-case behaviour is unchanged, since case 0's directory is exactly
the path that was used before.
Restore the DJF/JJA/MAM/SON entries in the Taylor script's seasons dict. They
were commented out in January's "refining taylor diagram" commit, apparently to
shorten a debug loop, so the branch produced only the annual diagram where main
produces five.
Verified: the rerun completes end to end. NaN in the regridded 3-D fields now
appears only below topography and above each model's own top -- the 32-level
case cuts off above 3.64 hPa, the 93-level case stays valid to 1 hPa, matching a
native-grid reference computation. All five seasons and all ten variables are
plotted for both cases, and the multi-case website builds with the Taylor pages
linked. adf_web.py is in the CI-gated lint set and holds at 9.70.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The masked entries on the diagram -- Land2mTemperature, TropicalLandPrecip, TropicalOceanPrecip -- carry different masks in the model and in the observations, because the two disagree about coastlines. taylor_stats_single then took each moment over a different population: the case mean and variance over the case's valid points, the reference mean and variance over the reference's, and the covariance over the intersection. That is not a correlation, and it does not have to land in [-1, 1]. Land2mTemperature came out at 1.0401, which never showed up as an error because plot_taylor_data clips the angle to corr <= 1, quietly parking the point on the REF ray. Mask both fields to the points where both are valid, before anything is computed, which is what NCL's taylor_stats effectively does. Doing it once at the top fixes the means, both variances, the covariance and the bias together rather than patching weighted_correlation alone. Verified against the two-case observation run: Land2mTemperature goes to 0.9965 and 0.9964 for the two cases, its bias from -0.26% to -0.02% now that both means cover the same points, and the tropical precipitation entries shift slightly. Every unmasked variable is unchanged to four decimal places, which is the control -- with nothing masked the intersection is everything. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Where a run carries PMID (or PINT for data on interfaces), that pressure field now takes precedence over rebuilding pressure from PS and the hybrid coefficients. It is the pressure the model actually used rather than a reconstruction, and it is the only option for vertical coordinates that are not hybrid-sigma, so treating it as the preferred source is where the rest of ADF has to end up. PS plus hyam/hybm remains the fallback. Three things the implementation has to get right: _find_pressure_field deliberately ignores "*_PMID_regridded.nc". PMID is a 3-D field, so when it is in diag_var_list that file has already been interpolated onto the output pressure levels and is useless as a source pressure -- the same shape of trap as reading PS back in hPa. It reads the climatology instead and horizontally regrids that, and normalizes units through _pressure_in_pa (renamed from _surface_pressure_in_pa, since it now serves both). Out-of-range levels are masked per column. vert_remap interpolates with np.interp, which clamps to the end values, while the geocat hybrid path returns NaN outside the column. Without the mask, changing the preferred pressure source would have quietly started reporting the lowest model level's value at 1000 hPa over Tibet instead of a NaN. Interface data is handled by renaming ilev to lev on the way in, because utils.pmid_to_plev stacks on a dimension named literally "lev". That path could not have worked for PINT before. The old "height" branch reached for source_ds['PMID'] on the native grid while its input had already been horizontally regridded, a latent grid mismatch. It is now unreachable, and reaching it raises rather than silently misaligning. Verified two ways. A unit test exploits vert_remap being exact for a field linear in log(p): exact recovery inside the column, NaN rather than clamping outside it, per-column masking, and the ilev path. And an end-to-end run with PMID in diag_var_list for a config where only one of the two cases writes PMID, so one case took the new path and the other the PS fallback: for the PMID case the regridded U agrees with the PS-based result to identical NaN masks at all 21 levels and RMS differences of 0.002-0.1 m/s against RMS winds of 4-33 m/s. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
create_website copies each case's website into the combined "main_website" with shutil.copytree and no dirs_exist_ok, so a second run into the same cam_diag_plot_loc dies on the directory the first run left behind: FileExistsError: .../plots/main_website/multi-case The line is identical in main; it was simply unreachable there, because multi-case website generation always died earlier on the missing run-info markdown. Fixing that exposed this. Verified by re-running the two-case observation config into a populated plot location: it now completes, refreshes both per-case sub-sites and the multi-case one, and every season's plot page is reachable from the season navigation. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
|
Fixed the seasonal plots. 
|
_determine_vertical_coord_type inspected the Dataset's dimensions, so any variable whose file happened to carry a 'lev' dimension was treated as a model-level field. _handle_vertical_interpolation then did da['lev'] on a 2-D DataArray and raised KeyError: 'lev', ending the run. This never fires with ADF's own ncrcat time series, which copy hyam/hybm only into files for variables that are actually on model levels. It fires for every GenTS run, because GenTS copies every non-time-varying variable - lev, ilev, hyam, hybm - into every file, so TS, PS and ICEFRAC all arrive in datasets carrying a lev dimension they do not use. GenTS TS dataset-only: 'hybrid' -> with the variable: 'none' GenTS ICEFRAC dataset-only: 'hybrid' -> with the variable: 'none' GenTS PS dataset-only: 'hybrid' -> with the variable: 'none' GenTS Q dataset-only: 'hybrid' -> with the variable: 'hybrid' The Dataset is still passed, because the hyam/hybm fallback check needs the other variables in the file; the DataArray now decides whether the variable needs interpolating at all. This is a fix to the regridding rewrite, found by running the two branches together - neither shows it alone. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
2 participants
The main objective of this PR is to provide Taylor Diagrams with observation/reanalysis reference data, and thus to address #426.
Since Taylor Diagrams are especially useful for comparing multiple simulations to a common reference, I also prioritized being able to run with multiple test cases.
Refactoring the original
cam_taylor_diagram.pymostly needed moving to use more ofadf_dataset. The other key part was allowing for the potential need to regrid obs-to-obs for derived variables that the Taylor diagram uses. I also refactored theregrid_and_vert_interp.pymodule. Both use xESMF for regridding, defaulting to conservative methods. The regridding weights are currently being stored in a subdirectory of the regrid location.For cases when we need to regrid from pressure levels to pressure levels, I added a plev_to_plev method in
adf_utils. Elsewhere in that module I made some modifications to try to keep metadata.I was running in to a hang when creating climo files, and I tracked it down to our use of multiprocessing. I think I've implemented a fix for that (which I think we've had to work around before). The main thing seems to be
mp.get_context('spawn')but there's also some dask stuff to try to prevent zombie processes being made. This seems to work, but isn't crucial for this PR.There are some updates to
adf_info.pythat help to deal with the multi-case configurations. This has to do with handling the specification of history files from the YAML file. This could break previously working configurations because it might make some assumptions about how those history files are specified.To read nested lists correctly (i.e., [['h0a', 'h1a'],['h0a', 'h1a']]), I also needed to modify
adf_config.py.I added the ERA5 land fraction to the observations repository and included it in the variable defaults file.
A few modifications are included in
adf_dataset, but mainly just trying to apply unit transformations more carefully. These might not be super necessary, but I don't think they will hurt anything.Added: vertical interpolation is now decided by the variable, not its file
Found while running this branch together with #467 (GenTS time series) on real
data. Neither branch shows it alone, so it is worth explaining.
_determine_vertical_coord_typeinspected the Dataset's dimensions. ADF's ownncrcattime series copyhyam/hybmonly into files for variables that aregenuinely on model levels, so the dataset dims and the variable's dims agree and
nothing goes wrong. GenTS copies every non-time-varying variable —
lev,ilev,hyam,hybm— into every file, so a 2-D field such asTSarrives in a datasetcarrying a
levdimension it does not use:levdimlevTS'none'TS'hybrid'ICEFRAC'hybrid'PS'hybrid'Q'hybrid'The mis-detection sent 2-D fields into hybrid interpolation, where
_handle_vertical_interpolationdoesda[lev_coord_name]on a DataArray that has nolev, raisingKeyError: 'lev'out ofcreate_time_seriesand ending the run with atraceback — which AGENTS.md §4.5 rules out for framework code.
The Dataset is still passed, because the
hyam/hybmfallback check needs the othervariables in the file; the DataArray now decides whether the variable needs
interpolating at all. With ncrcat output the answer is unchanged for every variable,
so this is a no-op for the existing back end.
Verified on real output from both back ends:
29 unit tests pass on this branch.
Note on the regridding fix in this PR
The same runs confirmed the problem this PR set out to fix. On
main, the test-caseclimo lookup globbed
{case}_{var}_*.ncwhile files are named{case}_{hist_str}_{var}_climo.nc, so nothing matched — three full ADF runs producedzero plots and still printed "ADF diagnostics has completed successfully", because
with no variables surviving regridding the plotters' "all plots already exist" check is
vacuously true. With this branch merged in, the same configs produce 280 plots each.