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 |
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.