From fe9f44e0b92b4c24ea78972ac05385e89ae0ab26 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 4 Sep 2026 17:07:52 -0600 Subject: [PATCH 1/5] Settle meridional and zonal re-runs from the file names Re-running the ADF with redo_plot false was doing far more than deciding there was nothing to do. Measured on a two-case, ten-variable run whose plots were all already present, the plotting stage took 13.6 s and rewrote five files. zonal_mean built the log-pressure file name two different ways: the pre-flight pass looked for '{var}_{season}_Zonal_logp_Mean.png' while the plotting pass wrote '{var}_logp_{season}_Zonal_Mean.png'. The name it looked for never existed, so every log-pressure plot was redrawn on every run whatever redo_plot said -- the five rewritten files, and the 3.6 s in matplotlib's _update_ticks and 3.1 s in mathtext.parse that dominated the profile. It also computed both seasonal means before finding out whether either plot was wanted. meridional_mean now follows zonal_mean, which is what #425 asks for: a pre-flight pass that registers the plots that already exist, and data read through AdfData rather than a hand-rolled glob of cam_regrid_loc. Because none of its file names depend on whether the variable has a 'lev' dimension, a variable or a case whose plots are all present is settled from the names alone, without opening a file. Reading through AdfData also stops the plots being scaled twice. The regridding stage applies the variable-defaults conversion when it writes and stamps 'transformed' on the result; this script applied scale_factor and add_offset again on top. PRECT was drawn 86400000 times too large and PS a hundred times too small. AdfData's _regrid_converters exists to make exactly that decision. The no-op run now takes 1.3 s and rewrites nothing. Of the thirty meridional plots, the fifteen for variables with no unit conversion are byte-identical; PRECT and PS change because they were wrong; RESTOM changes only its units label, from the variable defaults' "W m$^{-2}$" to the file's own "W/m2", because AdfData applies new_unit only alongside a conversion -- which is what zonal_mean has always shown. Fixes #425. Co-Authored-By: Claude Opus 5 (1M context) --- scripts/plotting/meridional_mean.py | 420 +++++++++++++++++----------- scripts/plotting/zonal_mean.py | 21 +- 2 files changed, 264 insertions(+), 177 deletions(-) diff --git a/scripts/plotting/meridional_mean.py b/scripts/plotting/meridional_mean.py index dc47854e2..9314e91f2 100644 --- a/scripts/plotting/meridional_mean.py +++ b/scripts/plotting/meridional_mean.py @@ -1,76 +1,66 @@ from pathlib import Path import numpy as np -import xarray as xr import plotting_functions as pf import adf_utils as utils import warnings # use to warn user about missing files. warnings.formatwarning = utils.my_formatwarning + def meridional_mean(adfobj): """ - This script plots meridional averages. - Follows the old AMWG convention of plotting 5S to 5N. - **Note:** the constraint of 5S to 5N is easily changed; - the function that calculates the average can take any range of latitudes. - Compare CAM climatologies against - other climatological data (observations or baseline runs). + Plots meridional average from climatological files (annual and seasonal). + + Follows the old AMWG convention of averaging 5S to 5N. **Note:** that + constraint is easily changed; the function that calculates the average + takes any range of latitudes. Compares CAM climatologies against other + climatological data (observations or baseline runs). + + Parameters + ---------- + adfobj : AdfDiag + The diagnostics object that contains all the configuration information + + Returns + ------- + None + Does not return value, produces files. + + Notes + ----- + Uses AdfData for loading data described by adfobj. + + Directly uses adfobj for the following: + plot_var_list, plot_location, climo_yrs, variable_defaults, + read_config_var, get_basic_info, add_website_data, debug_log + + Every plot this makes is named the same way whether or not the variable + has a `lev` dimension, so a case whose plots are all present can be + settled from the file names alone -- see the pre-flight pass below, which + is what keeps re-running the ADF cheap. """ #Notify user that script has started: msg = "\n Generating meridional mean plots..." print(f"{msg}\n {'-' * (len(msg)-3)}") - #Extract needed quantities from ADF object: - #----------------------------------------- var_list = adfobj.plot_var_list - model_rgrid_loc = adfobj.get_basic_info("cam_regrid_loc", required=True) #Special ADF variable which contains the output paths for #all generated plots and tables: plot_locations = adfobj.plot_location - #CAM simulation variables (this is always assumed to be a list): - case_names = adfobj.get_cam_info("cam_case_name", required=True) - #Grab case years syear_cases = adfobj.climo_yrs["syears"] eyear_cases = adfobj.climo_yrs["eyears"] - # CAUTION: - # "data" here refers to either obs or a baseline simulation, - # Until those are both treated the same (via intake-esm or similar) - # we will do a simple check and switch options as needed: - if adfobj.get_basic_info("compare_obs"): - #Set obs call for observation details for plot titles - obs = True - - #Extract variable-obs dictionary: - var_obs_dict = adfobj.var_obs_dict - - #If dictionary is empty, then there are no observations to regrid to, - #so quit here: - if not var_obs_dict: - print("\t No observations found to plot against, so no meridional-mean maps will be generated.") - return - else: - obs = False - data_name = adfobj.get_baseline_info("cam_case_name", required=True) # does not get used, is just here as a placemarker - data_list = [data_name] # gets used as just the name to search for climo files HAS TO BE LIST - data_loc = model_rgrid_loc #Just use the re-gridded model data path - #End if - #Grab baseline years (which may be empty strings if using Obs): syear_baseline = adfobj.climo_yrs["syear_baseline"] eyear_baseline = adfobj.climo_yrs["eyear_baseline"] - #Grab all case nickname(s) - test_nicknames = adfobj.case_nicknames["test_nicknames"] - base_nickname = adfobj.case_nicknames["base_nickname"] - res = adfobj.variable_defaults # will be dict of variable-specific plot preferences - # or an empty dictionary if use_defaults was not specified in the config YAML file. + # or an empty dictionary if use_defaults was not specified in YAML. #Set plot file type: # -- this should be set in basic_info_dict, but is not required @@ -84,13 +74,6 @@ def meridional_mean(adfobj): print(f"\t NOTE: redo_plot is set to {redo_plot}") #----------------------------------------- - #Set data path variables: - #----------------------- - mclimo_rg_loc = Path(model_rgrid_loc) - if not adfobj.compare_obs: - dclimo_loc = Path(data_loc) - #----------------------- - #Set seasonal ranges: seasons = {"ANN": np.arange(1,13,1), "DJF": [12, 1, 2], @@ -98,152 +81,211 @@ def meridional_mean(adfobj): "MAM": [3, 4, 5], "SON": [9, 10, 11]} + # Check if plots already exist and redo_plot boolean + # If redo_plot is false and file exists, keep track and skip the calculation + # entirely, which is what makes re-running the ADF fast. A set, because + # every season of every variable is looked up in it once per case. + merid_skip = set() + + # Loop over model cases: + for case_idx, case_name in enumerate(adfobj.data.case_names): + # Set output plot location: + plot_loc = Path(plot_locations[case_idx]) + + # Check if plot output directory exists, and if not, then create it: + if not plot_loc.is_dir(): + print(f" {plot_loc} not found, making new directory") + plot_loc.mkdir(parents=True) + # End if + + # Loop over the variables for each season + for var in var_list: + for s in seasons: + plot_name = plot_loc / f"{var}_{s}_Meridional_Mean.{plot_type}" + + # Check redo_plot. If set to True: remove old plot, if it already exists: + if (not redo_plot) and plot_name.is_file(): + merid_skip.add(plot_name) + # Add already-existing plot to website (if enabled): + adfobj.add_website_data( + plot_name, var, case_name, season=s, plot_type="Meridional" + ) + continue + elif (redo_plot) and plot_name.is_file(): + plot_name.unlink() + # End if + # End for (seasons) + # End for (variables) + # End for (cases) + # + # End redo plots check + # + + # + # Setup Plotting + # #Loop over variables: for var in var_list: #Notify user of variable being plotted: print(f"\t - meridional mean plots for {var}") - if adfobj.compare_obs: - #Check if obs exist for the variable: - if var in var_obs_dict: - #Note: In the future these may all be lists, but for - #now just convert the target_list. - #Extract target file: - dclimo_loc = var_obs_dict[var]["obs_file"] - #Extract target list (eventually will be a list, for now need to convert): - data_list = [var_obs_dict[var]["obs_name"]] - #Extract target variable name: - data_var = var_obs_dict[var]["obs_var"] - else: - dmsg = f"No obs found for variable `{var}`, meridional mean plotting skipped." - adfobj.debug_log(dmsg) - continue - #End if - else: - #Set "data_var" for consistent use below: - data_var = var + # Nothing to draw for this variable at all: settle it from the file + # names, before opening any data + if _all_plots_present( + plot_locations, adfobj.data.case_names, var, seasons, plot_type, merid_skip + ): + continue + # End if + + if var not in adfobj.data.ref_var_nam: + dmsg = f"\t WARNING: No reference data found for variable `{var}`, meridional mean plotting skipped." + adfobj.debug_log(dmsg) + print(dmsg) + continue #End if # Check res for any variable specific options that need to be used BEFORE going to the plot: if var in res: vres = res[var] - #If found then notify user, assuming debug log is enabled: - adfobj.debug_log(f"meridional_mean: Found variable defaults for {var}") - + # If found then notify user, assuming debug log is enabled: + adfobj.debug_log( + f"\t INFO: meridional_mean: Found variable defaults for {var}" + ) else: vres = {} #End if - #loop over different data sets to plot model against: - for data_src in data_list: - # load data (observational) comparison files - # (we should explore intake as an alternative to having this kind of repeated code): - if adfobj.compare_obs: - #For now, only grab one file (but convert to list for use below) - oclim_fils = [dclimo_loc] - else: - oclim_fils = sorted(dclimo_loc.glob(f"{data_src}_{var}_baseline.nc")) + # load reference data (observational or baseline) + if not adfobj.compare_obs: + base_name = adfobj.data.ref_case_label + else: + base_name = adfobj.data.ref_labels[var] + # End if + + # Gather reference variable data + odata = adfobj.data.load_reference_regrid_da(base_name, var) + + # Check if regridded file exists, if not skip meridional plot for this var + if odata is None: + dmsg = f"\t WARNING: No regridded baseline file for {base_name} for variable `{var}`, meridional mean plotting skipped." + adfobj.debug_log(dmsg) + continue + # End if + + # Check meridional mean dimensions + has_dims_ref = utils.validate_dims(odata, ["lat", "lev"]) + + # check if there is a lat dimension: + # if not, skip test cases and move to next variable + if not has_dims_ref["has_lat"]: + print( + f"\t WARNING: Variable {var} is missing a lat dimension for '{base_name}', cannot continue to plot." + ) + continue + # End if + + # Loop over model cases: + for case_idx, case_name in enumerate(adfobj.data.case_names): + + # Set case nickname: + case_nickname = adfobj.data.test_nicknames[case_idx] + + # Set output plot location: + plot_loc = Path(plot_locations[case_idx]) + + # Nothing to draw for this case: settle it before opening the file + if all( + (plot_loc / f"{var}_{s}_Meridional_Mean.{plot_type}") in merid_skip + for s in seasons + ): + continue #End if - oclim_ds = utils.load_dataset(oclim_fils) - #Loop over model cases: - for case_idx, case_name in enumerate(case_names): + # load re-gridded model files: + mdata = adfobj.data.load_regrid_da(case_name, var) + + if mdata is None: + dmsg = f"\t WARNING: No regridded test file for {case_name} for variable `{var}`, meridional mean plotting skipped." + adfobj.debug_log(dmsg) + continue + # End if - #Set case nickname: - case_nickname = test_nicknames[case_idx] + # determine whether it's 2D or 3D + # 3D triggers search for surface pressure + has_dims = utils.validate_dims(mdata, ["lat", "lev"]) - #Set output plot location: - plot_loc = Path(plot_locations[case_idx]) + # check if there is a lat dimension: + if not has_dims["has_lat"]: + print( + f"\t WARNING: Variable {var} is missing a lat dimension for '{case_name}', cannot continue to plot." + ) + continue + # End if - #Check if plot output directory exists, and if not, then create it: - if not plot_loc.is_dir(): - print(f" {plot_loc} not found, making new directory") - plot_loc.mkdir(parents=True) + has_lev = has_dims["has_lev"] - # load re-gridded model files: - mclim_fils = sorted(mclimo_rg_loc.glob(f"{data_src}_{case_name}_{var}_*.nc")) - mclim_ds = utils.load_dataset(mclim_fils) + # Notify user of level dimension: + if has_lev: + print(f"\t INFO: {var} has lev dimension.") + # End if - # stop if data is invalid: - if (oclim_ds is None) or (mclim_ds is None): - warnings.warn(f"invalid data, skipping meridional mean plot of {var}") + # Check to make sure each case has vertical levels if one of the cases does + if has_lev != has_dims_ref["has_lev"]: + print( + f"\t WARNING: expecting lev boolean for both case: {has_lev} and ref: {has_dims_ref['has_lev']}" + ) + continue + # End if + + # + # Seasonal Averages + # Note: xarray can do seasonal averaging, but depends on having time accessor, + # which these prototype climo files don't. + # + + # Create new dictionaries: + mseasons = {} + oseasons = {} + + # Loop over season dictionary: + for s in seasons: + # Set the file name + plot_name = plot_loc / f"{var}_{s}_Meridional_Mean.{plot_type}" + + # Found above and redo_plot is false, so there is nothing to + # draw. The seasonal averages below are the expensive part of + # a re-run, so work that out before doing them. + if plot_name in merid_skip: continue + # End if - #Extract variable of interest - odata = oclim_ds[data_var].squeeze() # squeeze in case of degenerate dimensions - mdata = mclim_ds[var].squeeze() - - # APPLY UNITS TRANSFORMATION IF SPECIFIED: - # NOTE: looks like our climo files don't have all their metadata - mdata = mdata * vres.get("scale_factor",1) + vres.get("add_offset", 0) - # update units - mdata.attrs['units'] = vres.get("new_unit", mdata.attrs.get('units', 'none')) - - # Do the same for the baseline case if need be: - if not adfobj.compare_obs: - odata = odata * vres.get("scale_factor",1) + vres.get("add_offset", 0) - # update units - odata.attrs['units'] = vres.get("new_unit", odata.attrs.get('units', 'none')) - # Or for observations - else: - odata = odata * vres.get("obs_scale_factor",1) + vres.get("obs_add_offset", 0) - # Note: we are going to assume that the specification ensures the conversion makes the units the same. - # Doesn't make sense to add a different unit. - - # determine whether it's 2D or 3D - # 3D triggers search for surface pressure - validate_lat_lev = utils.validate_dims(mdata, ['lat', 'lev']) # keys=> 'has_lat', 'has_lev', with T/F values - - #Notify user of level dimension: - if validate_lat_lev['has_lev']: - print(f"\t INFO: {var} has lev dimension.") - has_lev = True - else: - has_lev = False - - # # Seasonal Averages - # Note: xarray can do seasonal averaging, but depends on having time accessor, - # which these prototype climo files don't. - # - - #Create new dictionaries: - mseasons = {} - oseasons = {} - - #Loop over season dictionary: - for s in seasons: - plot_name = plot_loc / f"{var}_{s}_Meridional_Mean.{plot_type}" - - # Check redo_plot. If set to True: remove old plot, if it already exists: - if (not redo_plot) and plot_name.is_file(): - #Add already-existing plot to website (if enabled): - adfobj.debug_log(f"'{plot_name}' exists and clobber is false.") - adfobj.add_website_data(plot_name, var, case_name, season=s, - plot_type="Meridional") - #Continue to next iteration: - continue - elif (redo_plot) and plot_name.is_file(): - plot_name.unlink() - - mseasons[s] = utils.seasonal_mean(mdata, season=s, is_climo=True) - oseasons[s] = utils.seasonal_mean(odata, season=s, is_climo=True) - - - #Create new plot: - pf.plot_meridional_mean_and_save(plot_name, case_nickname, base_nickname, - [syear_cases[case_idx],eyear_cases[case_idx]], - [syear_baseline,eyear_baseline], - mseasons[s], oseasons[s], has_lev, latbounds=slice(-5,5), obs=obs, **vres) - - #Add plot to website (if enabled): - adfobj.add_website_data(plot_name, var, case_name, season=s, - plot_type="Meridional") - - #End for (seasons loop) - #End for (case names loop) - #End for (obs/baseline loop) - #End for (variables loop) + mseasons[s] = utils.seasonal_mean(mdata, season=s, is_climo=True) + oseasons[s] = utils.seasonal_mean(odata, season=s, is_climo=True) + + # Create new plot: + pf.plot_meridional_mean_and_save( + plot_name, + case_nickname, + adfobj.data.ref_nickname, + [syear_cases[case_idx], eyear_cases[case_idx]], + [syear_baseline, eyear_baseline], + mseasons[s], + oseasons[s], + has_lev, + latbounds=slice(-5, 5), + obs=adfobj.compare_obs, + **vres, + ) + + # Add plot to website (if enabled): + adfobj.add_website_data( + plot_name, var, case_name, season=s, plot_type="Meridional" + ) + + # End for (seasons loop) + # End for (case names loop) + # End for (variables loop) #Notify user that script has ended: print(" ...Meridional mean plots have been generated successfully.") @@ -254,5 +296,41 @@ def meridional_mean(adfobj): ######### +def _all_plots_present(plot_locations, case_names, var, seasons, plot_type, merid_skip): + """ + Report whether every plot this variable would produce is already there. + + Parameters + ---------- + plot_locations : list + output plot directory for each case + case_names : list + names of the test cases + var : str + ADF name of the variable being plotted + seasons : dict + the seasons being plotted, keyed by name + plot_type : str + file extension the plots are written with, e.g. "png" + merid_skip : set + plots found by the pre-flight pass above + + Returns + ------- + bool + ``True`` when nothing is left to draw for this variable, so the + reference and test data never have to be opened. + """ + for case_idx in range(len(case_names)): + plot_loc = Path(plot_locations[case_idx]) + for s in seasons: + if (plot_loc / f"{var}_{s}_Meridional_Mean.{plot_type}") not in merid_skip: + return False + # End if + # End for + # End for + return True + + ############## -#END OF SCRIPT \ No newline at end of file +# END OF FILE diff --git a/scripts/plotting/zonal_mean.py b/scripts/plotting/zonal_mean.py index 2f8630fb3..97a19928f 100644 --- a/scripts/plotting/zonal_mean.py +++ b/scripts/plotting/zonal_mean.py @@ -100,7 +100,7 @@ def zonal_mean(adfobj): for var in var_list: for s in seasons: #Check zonal log-p: - plot_name_log = plot_loc / f"{var}_{s}_Zonal_logp_Mean.{plot_type}" + plot_name_log = plot_loc / f"{var}_logp_{s}_Zonal_Mean.{plot_type}" # Check redo_plot. If set to True: remove old plot, if it already exists: if (not redo_plot) and plot_name_log.is_file(): @@ -252,10 +252,6 @@ def zonal_mean(adfobj): # because we can let any pressure-level interpolation happen there # This could be re-visited for efficiency or improved code structure. - #Seasonal Averages - mseasons[s] = utils.seasonal_mean(mdata, season=s, is_climo=True) - oseasons[s] = utils.seasonal_mean(odata, season=s, is_climo=True) - #Set the file name plot_name = plot_loc / f"{var}_{s}_Zonal_Mean.{plot_type}" plot_name_log = None @@ -265,6 +261,19 @@ def zonal_mean(adfobj): plot_name_log = plot_loc / f"{var}_logp_{s}_Zonal_Mean.{plot_type}" #End if + # Nothing to draw for this season: both plots were found above + # and redo_plot is false. The seasonal averages below are the + # expensive part of a re-run, so work that out before doing them. + if (plot_name in zonal_skip) and ( + (plot_name_log is None) or (plot_name_log in logp_zonal_skip) + ): + continue + # End if + + # Seasonal Averages + mseasons[s] = utils.seasonal_mean(mdata, season=s, is_climo=True) + oseasons[s] = utils.seasonal_mean(odata, season=s, is_climo=True) + #Create plots if plot_name not in zonal_skip: @@ -299,4 +308,4 @@ def zonal_mean(adfobj): ############## -#END OF SCRIPT \ No newline at end of file +#END OF SCRIPT From 4e07a1f8a958edd4bb8075636d0ab0a5dbe5650d Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 4 Sep 2026 19:04:49 -0600 Subject: [PATCH 2/5] Decide unit conversions the same way everywhere they are loaded MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The ADF converts units with the scale_factor, add_offset and new_unit in the variable defaults, and it has to apply that conversion exactly once. Whether it had already been applied was being decided seven different ways. load_da and the time series loaders applied it unconditionally. load_climo_ds applied it unconditionally and stamped 'transformed' on the result. load_reference_climo_ds checked that stamp, and also compared the file's units to new_unit -- with ==, which almost never matched, because the two strings come from different hands: CAM writes 'W/m2' where the defaults say 'Wm$^{-2}$'. load_reference_climo_da checked nothing, on the same files. Only the regrid loaders had a working guard, and it read the stamp alone. There is now one decision, AdfData.already_converted, used by all of them: the stamp if it is there, and otherwise whether the file's units already say what the conversion is converting to. adf_units carries that comparison. It reduces a unit to a canonical form first, so the many spellings of one unit agree: 'W/m2', 'W m-2', 'Wm^-2', 'W m**-2', 'W m⁻²' and the LaTeX 'W m$^{-2}$' and 'Wm$^{-2}$' -- both of which the shipped defaults contain -- are one unit written seven ways. It folds a solidus into negative exponents, splits a run-together factor into known symbols, reads superscripts, and treats a unit whose factors cancel as dimensionless, so 'kg/kg', 'kg kg-1' and 'fraction' agree too. It imports nothing but re, so CI runs its tests. Two plotting scripts scaled by hand on data the regridding stage had already converted, the same bug meridional_mean had: global_latlon_vect_map on regridded climatologies, and tem. Both now go through AdfData.apply_conversion, which makes the same decision the loaders do. global_latlon_vect_map also relabelled units from the defaults whether or not it had converted, claiming units the data was not in. Nothing that was right changes: the three AMWG tables are byte-identical to what main produces, and a plotting re-run still rewrites nothing. Co-Authored-By: Claude Opus 5 (1M context) --- lib/adf_dataset.py | 137 +++++++++- lib/adf_units.py | 289 +++++++++++++++++++++ lib/adf_utils.py | 4 + lib/test/unit_tests/test_adf_units.py | 171 ++++++++++++ lib/test/unit_tests/test_regrid_scaling.py | 25 +- scripts/plotting/global_latlon_vect_map.py | 26 +- scripts/plotting/tem.py | 18 +- 7 files changed, 638 insertions(+), 32 deletions(-) create mode 100644 lib/adf_units.py create mode 100644 lib/test/unit_tests/test_adf_units.py diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py index 6caedee6f..194ffe9ab 100644 --- a/lib/adf_dataset.py +++ b/lib/adf_dataset.py @@ -5,6 +5,7 @@ import xarray as xr import adf_utils as utils +from adf_units import units_equivalent warnings.formatwarning = utils.my_formatwarning # "reference data" @@ -294,6 +295,16 @@ def load_climo_ds(self, case, variablename): ds = self.load_dataset(fils) if ds is None: return None + if (scale_factor != 1 or add_offset != 0) and self.already_converted( + ds[variablename].attrs, variablename + ): + self.adf.debug_log( + f"\t INFO: '{variablename}' climo file is already" + " in converted units; not converting again." + ) + scale_factor = 1 + add_offset = 0 + # End if # xarray arithmetic drops attrs, so carry them across by hand -- otherwise # the regridded files lose 'units' and the plotting scripts KeyError on it. attrs = ds[variablename].attrs.copy() @@ -368,10 +379,10 @@ def load_reference_climo_ds(self, case, variablename, apply_scaling=True): if ds is None: return None vname = self.ref_var_nam[variablename] # name of variable in the reference data - # Check if already transformed (via attribute or units) - new_unit = self.adf.variable_defaults.get(variablename, {}).get('new_unit') - unit_match = new_unit is not None and ds[vname].attrs.get('units') == new_unit - if ds[vname].attrs.get('transformed', False) or unit_match: + # Check if already transformed. The units comparison was a literal + # string match, which almost never fired: an observation file says + # "W/m2" where the defaults say "Wm$^{-2}$". + if self.already_converted(ds[vname].attrs, variablename): apply_scaling = False if not apply_scaling: add_offset = 0 @@ -394,7 +405,13 @@ def load_reference_climo_da(self, case, variablename, apply_scaling=True): scale_factor = 1 else: add_offset, scale_factor = self.get_value_converters(case, variablename) - return self.load_da(fils, vname, add_offset=add_offset, scale_factor=scale_factor) + return self.load_da( + fils, + vname, + field=variablename, + add_offset=add_offset, + scale_factor=scale_factor, + ) def get_reference_climo_file(self, var): """Return a list of files to be used as reference (aka baseline) for variable var.""" @@ -532,7 +549,7 @@ def _regrid_converters(self, fils, file_field, case, field, apply_scaling): if apply_scaling or (scale_factor == 1 and add_offset == 0): return add_offset, scale_factor ds = self.load_dataset(fils) - if ds is not None and ds[file_field].attrs.get('transformed', 0): + if ds is not None and self.already_converted(ds[file_field].attrs, field): return 0, 1 return add_offset, scale_factor @@ -569,10 +586,101 @@ def load_dataset(self, fils, use_time_bounds=False): # End if return ds - def load_da(self, fils, variablename, use_time_bounds=False, **kwargs): + def already_converted(self, attrs, field): + """Report whether the conversion for `field` has already been applied. + + Parameters + ---------- + attrs : dict + attributes of the variable as it was read from the file + field : str + ADF name of the variable, i.e. the key into the variable defaults + + Returns + ------- + bool + ``True`` when the file already holds converted values. + + Notes + ----- + Two things say so. A file the ADF wrote carries ``transformed``, + stamped by whichever load applied the conversion. A file it did not + write -- an observation file, or a regridded file from an older ADF -- + carries no stamp, and then the only evidence is that its units are + already the units the defaults are converting to. That comparison is + made with :func:`units_equivalent`, because the two strings come from + different hands: CAM writes ``W/m2`` where the defaults say + ``Wm$^{-2}$``, and comparing them literally answers the wrong + question. + """ + if attrs.get("transformed", False): + return True + # End if + new_unit = self.adf.variable_defaults.get(field, {}).get("new_unit") + return units_equivalent(attrs.get("units"), new_unit) + + def apply_conversion(self, data, field, case=None): + """Apply the variable-defaults conversion to `data`, once. + + For a script that opens a file itself rather than through one of the + loaders here -- the TEM files and the vector plots do -- so that it + makes the same decision they do. + + Parameters + ---------- + data : xarray.DataArray + the variable as it was read from the file + field : str + ADF name of the variable, i.e. the key into the variable defaults + case : str, optional + case the data belongs to, which decides whether the observation + converters are used. Defaults to `field`'s test-case converters. + + Returns + ------- + xarray.DataArray + The converted data, carrying its attributes, with ``units`` set + to the defaults' ``new_unit`` and ``transformed`` stamped on when + a conversion was applied. Returned unchanged when the file + already holds converted values. + """ + if data is None: + return None + # End if + add_offset, scale_factor = self.get_value_converters( + case if case is not None else field, field + ) + if scale_factor == 1 and add_offset == 0: + return data + # End if + if self.already_converted(data.attrs, field): + dmsg = f"\t INFO: '{field}' is already in converted units in the" + dmsg += " file, so the conversion is not applied again." + self.adf.debug_log(dmsg) + return data + # End if + attrs = data.attrs.copy() + converted = data * scale_factor + add_offset + converted.attrs = attrs + new_unit = self.adf.variable_defaults.get(field, {}).get("new_unit") + if new_unit: + converted.attrs["units"] = new_unit + # End if + # int, not bool: netCDF4 cannot store a Python bool as an attribute + converted.attrs["transformed"] = 1 + return converted + + def load_da(self, fils, variablename, use_time_bounds=False, field=None, **kwargs): """Return xarray DataArray from file(s) w/ optional scale factor, offset, new units. `use_time_bounds` is passed to `load_dataset`; see there. + + `field` is the ADF name of the variable when it differs from its name + in the file, which is the case for observations. The variable + defaults are keyed by the ADF name. + + A conversion that has already been applied to the file is not applied + again; see :meth:`already_converted`. """ ds = self.load_dataset(fils, use_time_bounds=use_time_bounds) if ds is None: @@ -581,12 +689,25 @@ def load_da(self, fils, variablename, use_time_bounds=False, **kwargs): da = ds[variablename].squeeze() scale_factor = kwargs.get('scale_factor', 1) add_offset = kwargs.get('add_offset', 0) + + if (scale_factor != 1 or add_offset != 0) and self.already_converted( + da.attrs, field if field is not None else variablename + ): + dmsg = f"\t INFO: '{variablename}' is already in converted units" + dmsg += " in the file, so the conversion is not applied again." + self.adf.debug_log(dmsg) + scale_factor = 1 + add_offset = 0 + # End if + attrs = da.attrs.copy() da = da * scale_factor + add_offset da.attrs = attrs if scale_factor != 1 or add_offset != 0: - new_unit = self.adf.variable_defaults.get(variablename, {}).get("new_unit") + new_unit = self.adf.variable_defaults.get( + field if field is not None else variablename, {} + ).get("new_unit") if new_unit: da.attrs['units'] = new_unit # Stamp on any conversion, not only one that renames the units -- diff --git a/lib/adf_units.py b/lib/adf_units.py new file mode 100644 index 000000000..7285993be --- /dev/null +++ b/lib/adf_units.py @@ -0,0 +1,289 @@ +""" +Unit string handling. + +Kept apart from adf_utils so it can be unit tested without the scientific +stack: the ADF unit test workflow installs only PyYAML and pytest, so anything +that imports xarray/geocat at module level cannot be exercised in CI. This +module imports nothing but re. + +Functions +--------- +normalize_units(units) + Reduce a unit string to a canonical form for comparison. +units_equivalent(first, second) + Report whether two unit strings mean the same thing. + +Notes +----- +The ADF compares units to decide whether the conversion named in the variable +defaults has already been applied to a file. The two strings being compared +come from different places -- one written by CAM, the other by whoever edited +`adf_variable_defaults.yaml` -- so they agree on the physics far more often +than they agree character for character. `W/m2`, `W m-2`, `Wm^-2` and the +LaTeX `W m$^{-2}$` are one unit written four ways, and the shipped defaults +alone contain both `Wm$^{-2}$` and `W m$^{-2}$`. Comparing the raw strings +answers "were these typed the same way", which is not the question being +asked, and getting it wrong scales the data twice. + +Rendering a unit for somewhere with no LaTeX renderer, such as a table cell, +is `adf_utils.plain_text_units`; this module only compares. +""" + +import re + +# Unit names that mean the same thing. Keys and values are compared after +# normalization, so only real synonyms belong here, not spelling variants: +_ALIASES = { + "fraction": "1", + "frac": "1", + "unitless": "1", + "dimensionless": "1", + "none": "1", + "-": "1", + "[-]": "1", + "percent": "%", + "pct": "%", + "degrees_kelvin": "k", + "deg_k": "k", + "kelvin": "k", + "degrees_celsius": "degc", + "celsius": "degc", + "deg_c": "degc", + "degrees_east": "degrees_east", + "meters": "m", + "meter": "m", + "metre": "m", + "seconds": "s", + "second": "s", + "sec": "s", + "days": "d", + "day": "d", + "grams": "g", + "gram": "g", + "micron": "um", + "microns": "um", +} + +# Unit symbols this knows how to recognise. Used to split a run-together +# factor such as "Wm-2" into "W" and "m-2": climate files write units both +# ways, and the shipped variable defaults contain both spellings. Note that +# "ms" is therefore read as metre-second rather than millisecond -- no CAM +# field is reported in milliseconds, and "ms$^{-1}$" for wind speed is in the +# defaults today: +_SYMBOLS = { + "w", + "m", + "s", + "k", + "g", + "kg", + "pa", + "hpa", + "j", + "n", + "mol", + "l", + "d", + "h", + "hr", + "yr", + "cm", + "mm", + "um", + "nm", + "km", + "rad", + "sr", + "ppb", + "ppm", + "ppt", + "ppbv", + "ppmv", + "pptv", + "%", + "c", + "v", + "a", + "1", +} + +# LaTeX and unicode fragments that carry no meaning for a comparison: +_LATEX = ( + (r"\mathrm", ""), + (r"\text", ""), + (r"\mu", "u"), + (r"\,", " "), + (r"\;", " "), + (r"\ ", " "), + ("µ", "u"), # micro sign + ("μ", "u"), # greek small letter mu + ("·", " "), # middle dot, used as a multiplication sign + ("**", "^"), +) + +# Superscript digits, which appear in units copied out of documents: +_SUPERSCRIPTS = str.maketrans( + { + "⁰": "0", + "¹": "1", + "²": "2", + "³": "3", + "⁴": "4", + "⁵": "5", + "⁶": "6", + "⁷": "7", + "⁸": "8", + "⁹": "9", + "⁻": "-", + } +) + + +def _strip_markup(units): + """Return `units` with LaTeX, superscripts and stray braces removed.""" + text = units.translate(_SUPERSCRIPTS) + for old, new in _LATEX: + text = text.replace(old, new) + # A LaTeX exponent is written ^{-2}; the braces are noise once the + # superscript digits are gone: + text = re.sub(r"\^\s*\{([^}]*)\}", r"^\1", text) + text = text.replace("$", "").replace("{", "").replace("}", "") + return text.strip() + + +def _split_symbols(name): + """ + Return `name` as a list of unit symbols, splitting a run-together factor. + + "wm" is watt-metre written without a space; "hpa" is a symbol in its own + right and must not become hecto-pascal. A name that cannot be covered + exactly by known symbols is left alone, so an unfamiliar unit still + compares equal to itself. + """ + if name in _SYMBOLS or name in _ALIASES: + return [name] + # End if + parts = [] + rest = name + while rest: + # Longest match first, so "kg" wins over "k": + for size in range(min(len(rest), 4), 0, -1): + if rest[:size] in _SYMBOLS: + parts.append(rest[:size]) + rest = rest[size:] + break + # End if + else: + return [name] + # End for + # End while + return parts if len(parts) > 1 else [name] + + +def _tokenize(text, sign): + """Return (name, exponent) pairs for one side of a solidus.""" + tokens = [] + for chunk in re.split(r"[\s*.]+", text): + if not chunk: + continue + # A trailing exponent, written either "m2", "m^2", "m-2" or "m^-2": + match = re.fullmatch(r"([a-z%_]+)\^?([+-]?\d+)?", chunk) + if match is None: + # Not something this understands; keep it whole so that two + # identical odd strings still compare equal: + tokens.append((chunk, sign)) + continue + name, exponent = match.group(1), match.group(2) + power = sign * int(exponent if exponent else 1) + # Only the last symbol of a run-together factor carries the exponent: + # "Wm-2" is watt per metre squared, not per watt per metre squared. + symbols = _split_symbols(name) + for symbol in symbols[:-1]: + tokens.append((symbol, sign)) + # End for + tokens.append((symbols[-1], power)) + return tokens + + +def normalize_units(units): + """ + Reduce a unit string to a canonical form for comparison. + + Parameters + ---------- + units : str or None + a unit string as it appears in a file or in the variable defaults + + Returns + ------- + str + A canonical form: lower case, no LaTeX, every factor written + ``name^exponent`` and sorted, so that the many spellings of one unit + reduce to a single string. An empty string for ``None`` or for a + string that holds nothing. + + Examples + -------- + ``W/m2``, ``W m-2``, ``Wm^-2`` and ``W m$^{-2}$`` all give ``m^-2 w^1``. + """ + if units is None: + return "" + text = _strip_markup(str(units)).lower() + if not text: + return "" + text = _ALIASES.get(text, text) + # Split on the solidus: everything after the first one is a denominator. + # "a/b/c" is read the way it is written, left to right. + parts = text.split("/") + tokens = _tokenize(parts[0], 1) + for part in parts[1:]: + tokens += _tokenize(part, -1) + + # Fold the aliases once more, now that the factors are separated, so that + # "m/s" and "meters/second" agree: + folded = {} + for name, exponent in tokens: + name = _ALIASES.get(name, name) + folded[name] = folded.get(name, 0) + exponent + # A factor that cancels out carries no information, and a unit whose + # factors all cancel is dimensionless -- "kg/kg" and "kg kg-1" are the same + # thing, and both are the same thing as "fraction": + remaining = { + name: exponent + for name, exponent in folded.items() + if exponent != 0 and name != "1" + } + if not remaining: + return "1" + # End if + return " ".join( + f"{name}^{exponent}" for name, exponent in sorted(remaining.items()) + ) + + +def units_equivalent(first, second): + """ + Report whether two unit strings mean the same thing. + + Parameters + ---------- + first, second : str or None + the unit strings to compare + + Returns + ------- + bool + ``True`` when the two describe the same unit, however they are + spelled. Two strings that are both empty or ``None`` are not + equivalent to anything, including each other: nothing is known about + a variable that does not say what its units are. + """ + left = normalize_units(first) + right = normalize_units(second) + if not left or not right: + return False + return left == right + + +############## +# END OF FILE diff --git a/lib/adf_utils.py b/lib/adf_utils.py index 858124578..7a01c50f3 100644 --- a/lib/adf_utils.py +++ b/lib/adf_utils.py @@ -7,6 +7,9 @@ ts_file_span(), as_hist_str_list(), pick_hist_str() re-exported from adf_file_utils; time series file discovery +normalize_units(), units_equivalent() + re-exported from adf_units; compare unit strings that are spelled + differently but mean the same thing plain_text_units() render a unit string from the variable defaults as plain text use_time_bounds_midpoint() @@ -63,6 +66,7 @@ #tested without importing the scientific stack (see adf_file_utils). Re-export #here so `utils.find_ts_files(...)` keeps working for every existing caller: # pylint: disable=unused-import +from adf_units import normalize_units, units_equivalent from adf_file_utils import ( as_hist_str_list, describe_dir_problem, diff --git a/lib/test/unit_tests/test_adf_units.py b/lib/test/unit_tests/test_adf_units.py new file mode 100644 index 000000000..4e798e200 --- /dev/null +++ b/lib/test/unit_tests/test_adf_units.py @@ -0,0 +1,171 @@ +""" +Collection of python unit tests for the "adf_units" unit string helpers. +""" + +# +++++++++++++++++++++++ +# Import required modules +# +++++++++++++++++++++++ + +import unittest +import sys +import os +import os.path + +# Set relevant path variables: +_CURRDIR = os.path.abspath(os.path.dirname(__file__)) +_ADF_LIB_DIR = os.path.join(_CURRDIR, os.pardir, os.pardir) + +# Add ADF "lib" directory to python path: +sys.path.append(_ADF_LIB_DIR) + +# adf_units imports nothing but re, so these run in CI, where only PyYAML and +# pytest are installed: +from adf_units import normalize_units, units_equivalent + +# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +# Main adf_units testing routine, used when script is run directly +# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + +class AdfUnitsTestRoutine(unittest.TestCase): + """ + Unit tests for comparing unit strings, which the ADF does to decide + whether a conversion has already been applied to a file. + """ + + def test_one_unit_written_many_ways(self): + """ + Check that the spellings of watts per square metre all agree. + + These are not hypothetical: CAM writes 'W/m2', and + adf_variable_defaults.yaml contains both 'Wm$^{-2}$' and + 'W m$^{-2}$'. + """ + + spellings = [ + "W/m2", + "W m-2", + "Wm-2", + "Wm^-2", + "W m^-2", + "W/m^2", + "W m**-2", + "W m$^{-2}$", + "Wm$^{-2}$", + "W m⁻²", + ] + + canonical = {normalize_units(u) for u in spellings} + + self.assertEqual( + len(canonical), 1, f"expected one canonical form, got {canonical}" + ) + for unit in spellings: + self.assertTrue(units_equivalent(unit, "W/m2"), unit) + + def test_defaults_spellings_match_file_spellings(self): + """ + Check the pairs that actually occur: a unit as CAM writes it against + the same unit as the variable defaults write it. + """ + + pairs = [ + ("W/m2", "Wm$^{-2}$"), + ("m/s", "ms$^{-1}$"), + ("mm/day", "mm d$^{-1}$"), + ("ug/m3", "$\\mu$g/m3"), + ("mol/mol", "mol mol$^{-1}$"), + ("fraction", "Fraction"), + ("%", "Percent"), + ] + + for from_file, from_defaults in pairs: + with self.subTest(units=from_file): + self.assertTrue(units_equivalent(from_file, from_defaults)) + + def test_different_units_stay_different(self): + """ + Check that units which differ are not reported as equivalent. + + This is the direction that matters for correctness: a false match + means a conversion is skipped and the data is plotted in the wrong + units. + """ + + pairs = [ + ("W/m2", "W/m3"), + ("m/s", "s/m"), + ("K", "W/m2"), + ("mm/day", "mm/s"), + ("Pa", "hPa"), + ("kg/m2", "kg/m3"), + ("ppbv", "ppmv"), + ("degrees_east", "degrees_north"), + ] + + for first, second in pairs: + with self.subTest(units=(first, second)): + self.assertFalse(units_equivalent(first, second)) + + def test_dimensionless_forms_agree(self): + """ + Check that a unit whose factors cancel is dimensionless. + """ + + for unit in [ + "kg/kg", + "kg kg-1", + "fraction", + "Fraction", + "1", + "unitless", + "none", + ]: + with self.subTest(units=unit): + self.assertEqual(normalize_units(unit), "1") + + def test_missing_units_match_nothing(self): + """ + Check that an absent unit is not equivalent to anything. + + Nothing is known about a variable that does not say what its units + are, so the ADF must not conclude a conversion was already applied. + """ + + self.assertEqual(normalize_units(None), "") + for missing in [None, "", " "]: + with self.subTest(units=missing): + self.assertFalse(units_equivalent(missing, "K")) + self.assertFalse(units_equivalent("K", missing)) + self.assertFalse(units_equivalent(missing, missing)) + + def test_unfamiliar_unit_matches_itself(self): + """ + Check that a unit this does not understand still compares equal to + itself, and unequal to something else. + """ + + self.assertTrue(units_equivalent("ppbv", "ppbv")) + self.assertTrue(units_equivalent("kg m-2 s-1", "kg/m2/s")) + self.assertFalse(units_equivalent("frobnicate", "ppbv")) + + def test_hpa_is_not_hecto_pascal_split(self): + """ + Check that a known symbol is not split into smaller ones. + + 'hPa' would become hour-pascal if the run-together splitting were + applied to a symbol that stands on its own. + """ + + self.assertEqual(normalize_units("hPa"), "hpa^1") + self.assertFalse(units_equivalent("hPa", "Pa")) + + +# ++++++++++++++++++ + +# Run unit tests if this script is called directly: +if __name__ == "__main__": + unittest.main() + +############# +# End of file diff --git a/lib/test/unit_tests/test_regrid_scaling.py b/lib/test/unit_tests/test_regrid_scaling.py index 08f85860c..576f9f9ee 100644 --- a/lib/test/unit_tests/test_regrid_scaling.py +++ b/lib/test/unit_tests/test_regrid_scaling.py @@ -25,11 +25,19 @@ NO_CONVERSION = (0, 1) -def _data(stamped, converters=CONVERTED): +def _data(stamped, converters=CONVERTED, units=None, new_unit=None): """An AdfData stand-in holding one variable, with or without the stamp.""" attrs = {'transformed': 1} if stamped else {} + if units is not None: + attrs["units"] = units + defaults = {"TAUX": {"new_unit": new_unit}} if new_unit else {} obj = SimpleNamespace( get_value_converters=lambda case, field: converters, + adf=SimpleNamespace(variable_defaults=defaults), + ) + # The real decision, so that this pins down what the loaders actually do: + obj.already_converted = lambda attrs_, field: AdfData.already_converted( + obj, attrs_, field ) obj.reads = 0 @@ -63,3 +71,18 @@ def test_variable_without_a_conversion_never_opens_the_file(): def test_explicit_override_wins_over_the_stamp(): assert _call(_data(stamped=False), apply_scaling=False) == NO_CONVERSION assert _call(_data(stamped=True), apply_scaling=True) == CONVERTED + + +def test_units_already_the_converted_ones_are_not_converted_again(): + """A file an older ADF wrote carries no stamp, but its units give it away. + + The comparison has to survive the two being spelled differently: CAM + writes "W/m2" where the variable defaults say "Wm$^{-2}$". + """ + obj = _data(stamped=False, units="W/m2", new_unit="Wm$^{-2}$") + assert _call(obj) == NO_CONVERSION + + +def test_units_that_differ_still_convert(): + obj = _data(stamped=False, units="m/s", new_unit="Wm$^{-2}$") + assert _call(obj) == CONVERTED diff --git a/scripts/plotting/global_latlon_vect_map.py b/scripts/plotting/global_latlon_vect_map.py index 403a0d8a8..d8f24a2cb 100644 --- a/scripts/plotting/global_latlon_vect_map.py +++ b/scripts/plotting/global_latlon_vect_map.py @@ -272,9 +272,12 @@ def global_latlon_vect_map(adfobj): uodata = uoclim_ds[data_var[0]].squeeze() # squeeze in case of degenerate dimensions vodata = voclim_ds[data_var[1]].squeeze() # squeeze in case of degenerate dimensions - #Convert units if requested (assumes units between model and data are the same): - uodata = uodata * vres.get("scale_factor",1) + vres.get("add_offset", 0) - vodata = vodata * vres.get("scale_factor",1) + vres.get("add_offset", 0) + # Convert units if requested (assumes units between model and data are the same). + # Through the ADF's data layer, so that a file already holding + # converted values -- which is what the regridding stage writes -- + # is not scaled a second time: + uodata = adfobj.data.apply_conversion(uodata, var) + vodata = adfobj.data.apply_conversion(vodata, var_pair) #Check zonal mean dimensions has_lat_ref, has_lev_ref = utils.zm_validate_dims(uodata) @@ -331,9 +334,9 @@ def global_latlon_vect_map(adfobj): umdata = umclim_ds[var].squeeze() vmdata = vmclim_ds[var_pair].squeeze() - #Convert units if requested: - umdata = umdata * vres.get("scale_factor",1) + vres.get("add_offset", 0) - vmdata = vmdata * vres.get("scale_factor",1) + vres.get("add_offset", 0) + # Convert units if requested, once -- see above: + umdata = adfobj.data.apply_conversion(umdata, var) + vmdata = adfobj.data.apply_conversion(vmdata, var_pair) #Check dimensions: has_lat, has_lev = utils.zm_validate_dims(umdata) @@ -346,12 +349,9 @@ def global_latlon_vect_map(adfobj): continue # End if - # update units - # NOTE: looks like our climo files don't have all their metadata - uodata.attrs['units'] = vres.get("new_unit", uodata.attrs.get('units', 'none')) - vodata.attrs['units'] = vres.get("new_unit", vodata.attrs.get('units', 'none')) - umdata.attrs['units'] = vres.get("new_unit", umdata.attrs.get('units', 'none')) - vmdata.attrs['units'] = vres.get("new_unit", vmdata.attrs.get('units', 'none')) + # Units are set by apply_conversion above, which renames them + # only when it actually converts. Relabelling here regardless + # claimed the defaults' units for data still in the file's own. #Determine if observations/baseline have the correct dimensions: if has_lev: @@ -520,4 +520,4 @@ def global_latlon_vect_map(adfobj): print(" ...lat/lon vector maps have been generated successfully.") ############## -#END OF SCRIPT \ No newline at end of file +#END OF SCRIPT diff --git a/scripts/plotting/tem.py b/scripts/plotting/tem.py index 13f9c49d4..63404f56a 100644 --- a/scripts/plotting/tem.py +++ b/scripts/plotting/tem.py @@ -408,19 +408,17 @@ def tem(adf): mdata = ds[var].squeeze() odata = ds_base[var].squeeze() - #Apply the unit conversion from the variable defaults. TEM - #files carry little metadata, so the new unit is taken from the - #defaults too. Observations have their own scaling, which is - #assumed to bring them to the same units, so they keep the unit - #string they arrived with. - mdata = mdata * vres.get("scale_factor", 1) + vres.get("add_offset", 0) - mdata.attrs['units'] = vres.get("new_unit", mdata.attrs.get('units', 'none')) + # Apply the unit conversion from the variable defaults, through + # the ADF's data layer so that a TEM file already holding + # converted values is not scaled a second time. Observations + # have their own scaling, which is assumed to bring them to the + # same units, so they keep the unit string they arrived with. + mdata = adfobj.data.apply_conversion(mdata, var) if obs: odata = (odata * vres.get("obs_scale_factor", 1) + vres.get("obs_add_offset", 0)) else: - odata = odata * vres.get("scale_factor", 1) + vres.get("add_offset", 0) - odata.attrs['units'] = vres.get("new_unit", odata.attrs.get('units', 'none')) + odata = adfobj.data.apply_conversion(odata, var) #End if #Month-length weighted seasonal (or annual) mean. The weighted @@ -571,4 +569,4 @@ def _mesh(da): print(" ...TEM plots have been generated successfully.") # Helper functions -################## \ No newline at end of file +################## From a45c8ecea258942116364eececeaaaa55871539b Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 4 Sep 2026 19:07:28 -0600 Subject: [PATCH 3/5] Let polar_map settle a re-run from the file names too polar_map already worked out that every plot existed and skipped the reference data, but it opened a regridded file for each case first, only to see whether the variable has a 'lev' dimension. A 2-D variable is plotted as {var}_{season}_{hemisphere}_Mean and a 3-D one as {var}_{pressure}hpa_{season}_{hemisphere}_Mean, so which of the two it is looking at can be read off the names, and the file does not have to be opened at all when nothing is missing. The plots still reach the website: the names are registered before the variable is skipped, which is what the loop below did with them. Co-Authored-By: Claude Opus 5 (1M context) --- scripts/plotting/polar_map.py | 103 +++++++++++++++++++++++++++++++++- 1 file changed, 101 insertions(+), 2 deletions(-) diff --git a/scripts/plotting/polar_map.py b/scripts/plotting/polar_map.py index 7893fee8e..e3b62645f 100644 --- a/scripts/plotting/polar_map.py +++ b/scripts/plotting/polar_map.py @@ -138,10 +138,32 @@ def polar_map(adfobj): vres = res.get(var, {}) web_category = vres.get("category", None) + # A complete set of plots can be recognised from the file names alone, + # because a 2-D variable and a 3-D one are named differently. Doing + # that first means a re-run never opens the regridded files just to + # find out which of the two it is looking at. + settled = _existing_plot_set( + plot_locations, case_names, var, seasons, pres_levs, plot_type + ) + if settled and not redo_plot: + for path, web_name, case_name, season, hemi_type in settled: + adfobj.add_website_data( + path, + web_name, + case_name, + category=web_category, + season=season, + plot_type=hemi_type, + ) + # End for + print(f"\t Skipping {var} - all plots already exist") + continue + # End if + # Get all plot info and check existence plot_info = [] all_plots_exist = True - + for case_idx, case_name in enumerate(case_names): plot_loc = Path(plot_locations[case_idx]) @@ -266,5 +288,82 @@ def polar_map(adfobj): ############## #END OF `polar_map` function + +def _existing_plot_set(plot_locations, case_names, var, seasons, pres_levs, plot_type): + """ + Return the complete set of polar plots for `var`, if one is on disk. + + A 2-D variable is plotted as ``{var}_{season}_{hemisphere}_Mean`` and a + 3-D one as ``{var}_{pressure}hpa_{season}_{hemisphere}_Mean``, so which + of the two a variable is can be read off the file names. That is what + lets a re-run decide there is nothing to do without opening the regridded + file to look for a ``lev`` dimension. + + Parameters + ---------- + plot_locations : list + output plot directory for each case + case_names : list + names of the test cases + var : str + ADF name of the variable being plotted + seasons : list or dict + the seasons being plotted + pres_levs : list + the configured pressure levels, empty when none are set + plot_type : str + file extension the plots are written with, e.g. "png" + + Returns + ------- + list + ``(path, website name, case, season, hemisphere)`` for every plot, + when a complete set exists. An empty list when it does not, which + means the data has to be opened to find out what is missing. + """ + for pressures in ([None], pres_levs): + found = [] + complete = bool(pressures) + for case_idx, case_name in enumerate(case_names): + plot_loc = Path(plot_locations[case_idx]) + for season in seasons: + for hemi_type in ["NHPolar", "SHPolar"]: + for pres in pressures: + if pres is None: + name = f"{var}_{season}_{hemi_type}_Mean.{plot_type}" + web_name = var + else: + name = ( + f"{var}_{pres}hpa_{season}_{hemi_type}" + f"_Mean.{plot_type}" + ) + web_name = f"{var}_{pres}hpa" + # End if + path = plot_loc / name + if not path.is_file(): + complete = False + break + # End if + found.append((path, web_name, case_name, season, hemi_type)) + # End for + if not complete: + break + # End if + # End for + if not complete: + break + # End if + # End for + if not complete: + break + # End if + # End for + if complete and found: + return found + # End if + # End for + return [] + + ############## -# END OF FILE \ No newline at end of file +# END OF FILE From 6bd01aa1a73e2e2fac9c2e10bb8a9505f5854318 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 4 Sep 2026 19:21:55 -0600 Subject: [PATCH 4/5] Settle three more plotting scripts from the file names global_latlon_vect_map opened four regridded files for a vector pair, and global_mean_timeseries opened the reference time series and every case's and took global means of them, before either found out that the plot it would draw was already there. aod_latlon opened both observation files and every case's climatology up front and only then checked each panel. All three name their output from things known before any data is read -- the vector's name and the configured pressure levels, the field name, the observation source and the season -- so a re-run can be settled from the names, and the plots still reach the website. The check that a named set is complete is now plotting_utils. first_complete_plot_set, shared with polar_map: a 2-D variable and a 3-D one are named differently, so trying both spellings says both that nothing needs drawing and which kind of variable it was, which is the only thing the data was being opened for. Verified against a run with U and V added: ten vector plots and five global-mean plots on the first pass, and a second pass that draws nothing, rewrites nothing and reports every variable as already done. Co-Authored-By: Claude Opus 5 (1M context) --- lib/plotting_utils.py | 40 +++++++++- scripts/plotting/aod_latlon.py | 68 +++++++++++++++- scripts/plotting/global_latlon_vect_map.py | 92 ++++++++++++++++++++++ scripts/plotting/global_mean_timeseries.py | 23 +++++- scripts/plotting/polar_map.py | 68 +++++++--------- 5 files changed, 248 insertions(+), 43 deletions(-) diff --git a/lib/plotting_utils.py b/lib/plotting_utils.py index e83836429..aa28fc1e0 100644 --- a/lib/plotting_utils.py +++ b/lib/plotting_utils.py @@ -1,8 +1,10 @@ -""" . +""". Generic plotting helper functions Functions --------- +first_complete_plot_set(layouts) + pick the set of plots that is already on disk, if there is one use_this_norm() switches matplotlib color normalization method get_difference_colors(values) @@ -27,6 +29,8 @@ """ #import statements: +from pathlib import Path + import numpy as np import xarray as xr import matplotlib as mpl @@ -72,6 +76,40 @@ def load_dataset(fils): ####### + +def first_complete_plot_set(layouts): + """ + Return the first set of plots whose files all exist. + + A plotting script can usually tell what it would produce from the + variable defaults and the configured seasons and pressure levels, without + opening any data -- except that a 2-D variable and a 3-D one are named + differently, and which one it has is a property of the data. Both + spellings can be checked instead: if either set is complete on disk, the + script knows both that nothing needs drawing and which kind of variable it + was looking at, and never has to open the file. + + Parameters + ---------- + layouts : list + candidate sets, most specific first. Each is a list of entries whose + first element is the plot's path; the rest is whatever the caller + needs to register the plot on the website. + + Returns + ------- + list + The first complete set, or an empty list when none is complete, which + means the data has to be opened to find out what is missing. + """ + for entries in layouts: + if entries and all(Path(entry[0]).is_file() for entry in entries): + return entries + # End if + # End for + return [] + + def use_this_norm(): """Just use the right normalization; avoids a deprecation warning.""" diff --git a/scripts/plotting/aod_latlon.py b/scripts/plotting/aod_latlon.py index 926dc884f..00e055516 100644 --- a/scripts/plotting/aod_latlon.py +++ b/scripts/plotting/aod_latlon.py @@ -32,7 +32,18 @@ class AODPlotConfig: def aod_latlon(adfobj): """Generate AOD comparison plots.""" config = AODPlotConfig() - + + # Every panel is named from the observation source and the season, both + # of which are known here. Settling a re-run first means neither the + # observation files nor any case's climatologies are opened, which is + # what this script spends its time on: + if _register_existing_panels(adfobj, config): + print( + "\t INFO: All AOD panels exist. Existing plots added to" " website data." + ) + return + # End if + # Load observations obs_data = load_observations(adfobj) if not obs_data: @@ -490,4 +501,57 @@ def monthly_to_seasonal(ds, obs=False): ds_seasonal['season'] = seasons ds_seasonal = ds_seasonal.transpose('lat', 'lon', 'season') - return ds_seasonal \ No newline at end of file + return ds_seasonal + + +def _register_existing_panels(adfobj, config): + """ + Put the AOD panels that are already on disk onto the website. + + Parameters + ---------- + adfobj : AdfDiag + The diagnostics object containing configuration + config : AODPlotConfig + the observation sources and seasons being plotted + + Returns + ------- + bool + ``True`` when every panel exists and ``redo_plot`` is false, so there + is nothing left to draw. ``False`` leaves the plots alone: the run + goes on to make them, and registers them as it does. + """ + if adfobj.get_basic_info("redo_plot"): + return False + # End if + file_type = adfobj.read_config_var("diag_basic_info").get("plot_type", "png") + plot_dir = Path(adfobj.plot_location[0]) + panels = [] + for obs_name in config.obs_sources: + for season in config.seasons: + label = obs_name.replace(" ", "_") + panels.append( + ( + plot_dir / f"AOD_diff_{label}_{season}_LatLon_Mean.{file_type}", + f"AOD_diff_{label}", + season, + ) + ) + # End for + # End for + if not all(path.is_file() for path, _, _ in panels): + return False + # End if + for path, web_name, season in panels: + adfobj.add_website_data( + path, + web_name, + None, + season=season, + multi_case=True, + plot_type="LatLon", + category="4-Panel AOD Diags", + ) + # End for + return True diff --git a/scripts/plotting/global_latlon_vect_map.py b/scripts/plotting/global_latlon_vect_map.py index d8f24a2cb..eb5457c1e 100644 --- a/scripts/plotting/global_latlon_vect_map.py +++ b/scripts/plotting/global_latlon_vect_map.py @@ -185,6 +185,32 @@ def global_latlon_vect_map(adfobj): # otherwise defaults to 180 vres['central_longitude'] = plot_utils.get_central_longitude(adfobj) + # A complete set of plots can be recognised from the file names alone, + # because a 2-D vector and a 3-D one are named differently. Doing that + # first means a re-run never opens the four regridded files a vector + # pair needs just to find out which of the two it is looking at. + settled = _existing_plot_set( + plot_locations, case_names, var_name, seasons, pres_levs, plot_type + ) + if settled and not redo_plot: + for path, web_name, case_name, season in settled: + adfobj.debug_log(f"'{path}' exists and clobber is false.") + adfobj.add_website_data( + path, + web_name, + case_name, + category=web_category, + season=season, + plot_type="LatLon_Vector", + ) + # End for + print( + f"\t INFO: All plots exist for {var_name}. " + f"Redo is {redo_plot}. Existing plots added to website data." + ) + continue + # End if + #Determine observations to compare against: if adfobj.compare_obs: if var not in adfobj.data.ref_var_nam: @@ -519,5 +545,71 @@ def global_latlon_vect_map(adfobj): #Notify user that script has ended: print(" ...lat/lon vector maps have been generated successfully.") + +def _existing_plot_set( + plot_locations, case_names, var_name, seasons, pres_levs, plot_type +): + """ + Return the complete set of vector plots for `var_name`, if one is on disk. + + A 2-D vector is plotted as ``{name}_{season}_LatLon_Vector_Mean`` and a + 3-D one as ``{name}_{pressure}hpa_{season}_LatLon_Vector_Mean``, so which + of the two a vector pair is can be read off the file names. + + Parameters + ---------- + plot_locations : list + output plot directory for each case + case_names : list + names of the test cases + var_name : str + the vector's name from the variable defaults, e.g. "Wind" + seasons : list or dict + the seasons being plotted + pres_levs : list + the configured pressure levels, empty when none are set + plot_type : str + file extension the plots are written with, e.g. "png" + + Returns + ------- + list + ``(path, website name, case, season)`` for every plot, when a + complete set exists; an empty list when it does not. + """ + import plotting_utils as plot_utils + + flat = [] + levelled = [] + for case_idx, case_name in enumerate(case_names): + plot_loc = Path(plot_locations[case_idx]) + for season in seasons: + flat.append( + ( + plot_loc / f"{var_name}_{season}_LatLon_Vector_Mean.{plot_type}", + var_name, + case_name, + season, + ) + ) + for lev in pres_levs: + levelled.append( + ( + plot_loc + / ( + f"{var_name}_{lev}hpa_{season}" + f"_LatLon_Vector_Mean.{plot_type}" + ), + f"{var_name}_{lev}hpa", + case_name, + season, + ) + ) + # End for + # End for + # End for + return plot_utils.first_complete_plot_set([flat, levelled]) + + ############## #END OF SCRIPT diff --git a/scripts/plotting/global_mean_timeseries.py b/scripts/plotting/global_mean_timeseries.py index 5ea5e9731..54ff1c9d3 100644 --- a/scripts/plotting/global_mean_timeseries.py +++ b/scripts/plotting/global_mean_timeseries.py @@ -34,6 +34,8 @@ def global_mean_timeseries(adfobj): # Gather ADF configurations plot_loc = get_plot_loc(adfobj) plot_type = adfobj.read_config_var("diag_basic_info").get("plot_type", "png") + redo_plot = adfobj.get_basic_info("redo_plot") + print(f"\t NOTE: redo_plot is set to {redo_plot}") res = adfobj.variable_defaults # will be dict of variable-specific plot preferences # or an empty dictionary if use_defaults was not specified in YAML. @@ -42,6 +44,24 @@ def global_mean_timeseries(adfobj): #Notify user of variable being plotted: print(f"\t - time series plot for {field}") + # This field's plot is named the same way whatever the data turns out + # to hold, so a re-run can settle it here -- before opening the + # reference time series and every case's, and before taking the + # global means of them, which is the whole cost of this script: + plot_name = plot_loc / f"{field}_GlobalMean_ANN_TimeSeries_Mean.{plot_type}" + if (not redo_plot) and plot_name.is_file(): + adfobj.debug_log(f"'{plot_name}' exists and clobber is false.") + adfobj.add_website_data( + plot_name, + f"{field}_GlobalMean", + None, + season="ANN", + multi_case=True, + plot_type="TimeSeries", + ) + continue + # End if + # Check res for any variable specific options that need to be used BEFORE going to the plot: if field in res: vres = res[field] @@ -162,7 +182,6 @@ def global_mean_timeseries(adfobj): unit = vres.get("new_unit","[-]") ax.set_ylabel(getattr(ref_ts_da,"unit", unit)) # add units - plot_name = plot_loc / f"{field}_GlobalMean_ANN_TimeSeries_Mean.{plot_type}" conditional_save(adfobj, plot_name, fig) @@ -313,4 +332,4 @@ def make_plot(case_ts, lens2, label=None, ref_ts_da=None): ############## -#END OF SCRIPT \ No newline at end of file +#END OF SCRIPT diff --git a/scripts/plotting/polar_map.py b/scripts/plotting/polar_map.py index e3b62645f..753ce659b 100644 --- a/scripts/plotting/polar_map.py +++ b/scripts/plotting/polar_map.py @@ -4,6 +4,7 @@ # ADF library import plotting_functions as pf +import plotting_utils as plot_utils import adf_utils as utils def get_hemisphere(hemi_type): @@ -318,51 +319,42 @@ def _existing_plot_set(plot_locations, case_names, var, seasons, pres_levs, plot ------- list ``(path, website name, case, season, hemisphere)`` for every plot, - when a complete set exists. An empty list when it does not, which - means the data has to be opened to find out what is missing. + when a complete set exists; an empty list when it does not. """ - for pressures in ([None], pres_levs): - found = [] - complete = bool(pressures) - for case_idx, case_name in enumerate(case_names): - plot_loc = Path(plot_locations[case_idx]) - for season in seasons: - for hemi_type in ["NHPolar", "SHPolar"]: - for pres in pressures: - if pres is None: - name = f"{var}_{season}_{hemi_type}_Mean.{plot_type}" - web_name = var - else: - name = ( + flat = [] + levelled = [] + for case_idx, case_name in enumerate(case_names): + plot_loc = Path(plot_locations[case_idx]) + for season in seasons: + for hemi_type in ["NHPolar", "SHPolar"]: + flat.append( + ( + plot_loc / f"{var}_{season}_{hemi_type}_Mean.{plot_type}", + var, + case_name, + season, + hemi_type, + ) + ) + for pres in pres_levs: + levelled.append( + ( + plot_loc + / ( f"{var}_{pres}hpa_{season}_{hemi_type}" f"_Mean.{plot_type}" - ) - web_name = f"{var}_{pres}hpa" - # End if - path = plot_loc / name - if not path.is_file(): - complete = False - break - # End if - found.append((path, web_name, case_name, season, hemi_type)) - # End for - if not complete: - break - # End if + ), + f"{var}_{pres}hpa", + case_name, + season, + hemi_type, + ) + ) # End for - if not complete: - break - # End if # End for - if not complete: - break - # End if # End for - if complete and found: - return found - # End if # End for - return [] + return plot_utils.first_complete_plot_set([flat, levelled]) ############## From 99d3ac8d1226c0d9dcdeb639b40c36cd71f81462 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 4 Sep 2026 20:07:10 -0600 Subject: [PATCH 5/5] Fix defects found reviewing the pre-flight and units work tem.py called adfobj.data.apply_conversion, but the entry point's argument is named adf, so the first TEM variable that needed drawing raised NameError and took the rest of the plots with it. A re-run whose plots all existed skipped the line, which is how the end-to-end testing missed it. normalize_units split a run-together factor by greedily matching known symbols, after lower-casing. That reported 'Nm-2' and 'nm-2' -- newton per metre squared and per square nanometre -- as the same unit, which is exactly the false match the module exists to prevent, and it shredded unit names that happen to spell other symbols: 'Sv' became siemens-volt, 'cal' became c-a-l. The splitting is now an explicit table of the forms that occur, applied while the case is still intact, and unit names it does not know are left whole. Multi-factor units written with several solidi -- 'kg/m2/s', 'W/m2/K' -- work either way round. A ratio whose factors cancel now keeps what cancelled, so 'kg/kg' and 'kg kg-1' still agree but 'kg/kg' and 'mol/mol' do not: a mass mixing ratio and a volume mixing ratio are both dimensionless and are not the same number. Also from the review: - load_reference_regrid_da did not pass 'field', so for observations the guard in load_da looked the obs variable name up in the variable defaults, found nothing, and could never fire. - load_reference_climo_ds stamped 'transformed' but never set the new units, so a converted reference kept the units it arrived with: the regridded baseline files on disk say 'm/s' over values in mm/day. - polar_map and global_latlon_vect_map passed the layouts to first_complete_plot_set least specific first, against what its docstring says, so a 3-D variable with stale 2-D-named plots could be declared complete and its levelled plots never drawn. - adf_units imports nothing but re and scores 10.00/10, so it goes in the set of files CI lints. - plot_meridional_mean_and_save assigned the reference's time mean to adata, throwing away the latitude weighting just applied to the test case. Latent, because seasonal_mean removes 'time' first, and fixed here because #425 is about the meridional plots. The AMWG tables are still byte-identical to main, a plotting re-run still rewrites nothing, and the thirty meridional plots are unchanged. Co-Authored-By: Claude Opus 5 (1M context) --- .github/scripts/pr_mod_file_tests.py | 21 ++-- lib/adf_dataset.py | 17 +++- lib/adf_units.py | 111 ++++++--------------- lib/plotting_functions.py | 10 +- lib/plotting_utils.py | 3 +- lib/test/unit_tests/test_adf_units.py | 68 ++++++++++--- scripts/plotting/global_latlon_vect_map.py | 2 +- scripts/plotting/meridional_mean.py | 3 +- scripts/plotting/polar_map.py | 2 +- scripts/plotting/tem.py | 4 +- 10 files changed, 127 insertions(+), 114 deletions(-) diff --git a/.github/scripts/pr_mod_file_tests.py b/.github/scripts/pr_mod_file_tests.py index 6b1ba1375..150febff8 100755 --- a/.github/scripts/pr_mod_file_tests.py +++ b/.github/scripts/pr_mod_file_tests.py @@ -140,15 +140,18 @@ def _main_prog(): print("Generating list of modified files...") - #This should eventually be passed in via a command-line - #argument, and include everything inside the "lib" directory -JN: - testable_files = {"lib/adf_base.py", - "lib/adf_config.py", - "lib/adf_file_utils.py", - "lib/adf_info.py", - "lib/adf_obs.py", - "lib/adf_web.py", - "lib/adf_diag.py"} + # This should eventually be passed in via a command-line + # argument, and include everything inside the "lib" directory -JN: + testable_files = { + "lib/adf_base.py", + "lib/adf_config.py", + "lib/adf_file_utils.py", + "lib/adf_info.py", + "lib/adf_obs.py", + "lib/adf_units.py", + "lib/adf_web.py", + "lib/adf_diag.py", + } #+++++++++++++++++++++++ #Read in input arguments diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py index 194ffe9ab..227e3d0c5 100644 --- a/lib/adf_dataset.py +++ b/lib/adf_dataset.py @@ -392,6 +392,13 @@ def load_reference_climo_ds(self, case, variablename, apply_scaling=True): ds[vname] = ds[vname] * scale_factor + add_offset ds[vname].attrs = attrs if scale_factor != 1 or add_offset != 0: + # Rename the units as load_climo_ds does for the test cases. Without + # this the reference kept the units it arrived with while holding + # converted values, so the file said 'm/s' over data in mm/day. + new_unit = self.adf.variable_defaults.get(variablename, {}).get("new_unit") + if new_unit: + ds[vname].attrs["units"] = new_unit + # End if # int, not bool: netCDF4 cannot store a Python bool as an attribute ds[vname].attrs['transformed'] = 1 return ds @@ -527,7 +534,15 @@ def load_reference_regrid_da(self, case, field, apply_scaling=None): file_field = self.ref_var_nam[field] if self.adf.compare_obs else field add_offset, scale_factor = self._regrid_converters(fils, file_field, case, field, apply_scaling) - return self.load_da(fils, file_field, add_offset=add_offset, scale_factor=scale_factor) + # field, not file_field: an observation file names the variable its own + # way, and the variable defaults are keyed by the ADF name + return self.load_da( + fils, + file_field, + field=field, + add_offset=add_offset, + scale_factor=scale_factor, + ) def _regrid_converters(self, fils, file_field, case, field, apply_scaling): """Return the (add_offset, scale_factor) to use for a regridded file. diff --git a/lib/adf_units.py b/lib/adf_units.py index 7285993be..46308d474 100644 --- a/lib/adf_units.py +++ b/lib/adf_units.py @@ -64,47 +64,23 @@ "microns": "um", } -# Unit symbols this knows how to recognise. Used to split a run-together -# factor such as "Wm-2" into "W" and "m-2": climate files write units both -# ways, and the shipped variable defaults contain both spellings. Note that -# "ms" is therefore read as metre-second rather than millisecond -- no CAM -# field is reported in milliseconds, and "ms$^{-1}$" for wind speed is in the -# defaults today: -_SYMBOLS = { - "w", - "m", - "s", - "k", - "g", - "kg", - "pa", - "hpa", - "j", - "n", - "mol", - "l", - "d", - "h", - "hr", - "yr", - "cm", - "mm", - "um", - "nm", - "km", - "rad", - "sr", - "ppb", - "ppm", - "ppt", - "ppbv", - "ppmv", - "pptv", - "%", - "c", - "v", - "a", - "1", +# Factors written without a space between them. Case matters and is still +# intact here, which is the point: "Nm" is newton-metre while "nm" is +# nanometre, and lower-casing first makes them the same string. Only forms +# that actually occur are listed -- a general splitter guessing where to cut a +# run of letters reads "Sv" as siemens-volt and "nm" as newton-metre, the +# second of which reports two different units as equal. Note that "ms" is +# read as metre-second: no CAM field is reported in milliseconds, and +# "ms$^{-1}$" for wind speed is in the shipped variable defaults today. +_RUN_TOGETHER = { + "Wm": "W m", + "Nm": "N m", + "Jm": "J m", + "Km": "K m", + "kgm": "kg m", + "gm": "g m", + "ms": "m s", + "Pas": "Pa s", } # LaTeX and unicode fragments that carry no meaning for a comparison: @@ -148,38 +124,14 @@ def _strip_markup(units): # superscript digits are gone: text = re.sub(r"\^\s*\{([^}]*)\}", r"^\1", text) text = text.replace("$", "").replace("{", "").replace("}", "") + # Split the run-together factors before anything lower-cases the string, + # because which factors they are depends on their case: + for joined, apart in _RUN_TOGETHER.items(): + text = re.sub(rf"(? 1 else [name] - - def _tokenize(text, sign): """Return (name, exponent) pairs for one side of a solidus.""" tokens = [] @@ -194,14 +146,7 @@ def _tokenize(text, sign): tokens.append((chunk, sign)) continue name, exponent = match.group(1), match.group(2) - power = sign * int(exponent if exponent else 1) - # Only the last symbol of a run-together factor carries the exponent: - # "Wm-2" is watt per metre squared, not per watt per metre squared. - symbols = _split_symbols(name) - for symbol in symbols[:-1]: - tokens.append((symbol, sign)) - # End for - tokens.append((symbols[-1], power)) + tokens.append((name, sign * int(exponent if exponent else 1))) return tokens @@ -245,16 +190,18 @@ def normalize_units(units): for name, exponent in tokens: name = _ALIASES.get(name, name) folded[name] = folded.get(name, 0) + exponent - # A factor that cancels out carries no information, and a unit whose - # factors all cancel is dimensionless -- "kg/kg" and "kg kg-1" are the same - # thing, and both are the same thing as "fraction": + # A factor that cancels out carries no dimension, but it does carry an + # identity: "kg/kg" and "kg kg-1" are the same thing, and neither is + # "mol/mol". A mass mixing ratio and a volume mixing ratio are both + # dimensionless and are not the same number, so what cancelled is kept. remaining = { name: exponent for name, exponent in folded.items() if exponent != 0 and name != "1" } if not remaining: - return "1" + cancelled = sorted(name for name in folded if name != "1") + return f"1[{' '.join(cancelled)}]" if cancelled else "1" # End if return " ".join( f"{name}^{exponent}" for name, exponent in sorted(remaining.items()) diff --git a/lib/plotting_functions.py b/lib/plotting_functions.py index 86911759e..88b9987e2 100644 --- a/lib/plotting_functions.py +++ b/lib/plotting_functions.py @@ -1210,8 +1210,12 @@ def plot_meridional_mean_and_save(wks, case_nickname, base_nickname, latweight = np.cos(np.radians(adata.lat)) adata = adata.weighted(latweight).mean(dim='lat', keep_attrs=True) if 'time' in bdata.dims: - adata = bdata.mean(dim='time', keep_attrs=True) - if 'lat' in bdata.dims: + # bdata, not adata: this took the reference's time mean and put it in + # the test case, throwing away the latitude weighting just applied to + # it. Latent, because seasonal_mean(is_climo=True) has already + # removed 'time' by the time the ADF calls this. + bdata = bdata.mean(dim="time", keep_attrs=True) + if "lat" in bdata.dims: latweight = np.cos(np.radians(bdata.lat)) bdata = bdata.weighted(latweight).mean(dim='lat', keep_attrs=True) # If there are other dimensions, they are still going to be there: @@ -1511,4 +1515,4 @@ def square_contour_difference(fld1, fld2, **kwargs): return fig ##################### -#END HELPER FUNCTIONS \ No newline at end of file +#END HELPER FUNCTIONS diff --git a/lib/plotting_utils.py b/lib/plotting_utils.py index aa28fc1e0..18e08ceb9 100644 --- a/lib/plotting_utils.py +++ b/lib/plotting_utils.py @@ -1,5 +1,4 @@ -""". -Generic plotting helper functions +"""Generic plotting helper functions. Functions --------- diff --git a/lib/test/unit_tests/test_adf_units.py b/lib/test/unit_tests/test_adf_units.py index 4e798e200..982e21d43 100644 --- a/lib/test/unit_tests/test_adf_units.py +++ b/lib/test/unit_tests/test_adf_units.py @@ -109,21 +109,65 @@ def test_different_units_stay_different(self): def test_dimensionless_forms_agree(self): """ - Check that a unit whose factors cancel is dimensionless. - """ - - for unit in [ - "kg/kg", - "kg kg-1", - "fraction", - "Fraction", - "1", - "unitless", - "none", - ]: + Check that the words for "no units" all mean the same thing. + + A ratio that cancels is dimensionless too, but keeps what cancelled -- + see test_dimensionless_ratios_keep_what_cancelled. + """ + + for unit in ["fraction", "Fraction", "1", "unitless", "none"]: with self.subTest(units=unit): self.assertEqual(normalize_units(unit), "1") + def test_case_decides_a_run_together_factor(self): + """ + Check that "Nm" and "nm" are not read as the same unit. + + Newton-metre and nanometre differ only in case, so the run-together + factors have to be split before anything lower-cases the string. + Reporting these two as equal would let a conversion be skipped on a + file whose units are nothing like the ones being converted to. + """ + + self.assertTrue(units_equivalent("N/m2", "Nm-2")) + self.assertFalse(units_equivalent("Nm-2", "nm-2")) + self.assertNotEqual(normalize_units("nm"), normalize_units("N m")) + + def test_unit_names_are_not_chopped_into_letters(self): + """ + Check that a unit name is left whole rather than read as a product of + whatever symbols happen to spell it. + + "Sv" is sieverts, not siemens-volt; "cal" is calories. + """ + + for unit in ["Sv", "cal", "dam", "DU", "molec"]: + with self.subTest(units=unit): + self.assertEqual(normalize_units(unit), unit.lower() + "^1") + self.assertFalse(units_equivalent("Sv", "vs")) + + def test_dimensionless_ratios_keep_what_cancelled(self): + """ + Check that two dimensionless ratios of different things differ. + + A mass mixing ratio and a volume mixing ratio are both dimensionless + and are not the same number. + """ + + self.assertTrue(units_equivalent("kg/kg", "kg kg-1")) + self.assertFalse(units_equivalent("kg/kg", "mol/mol")) + self.assertFalse(units_equivalent("kg/kg", "m3/m3")) + + def test_units_with_several_factors(self): + """ + Check units written with more than one solidus or exponent. + """ + + self.assertTrue(units_equivalent("kg/m2/s", "kg m-2 s-1")) + self.assertTrue(units_equivalent("W/m2/K", "W m-2 K-1")) + self.assertTrue(units_equivalent("1/s", "s-1")) + self.assertFalse(units_equivalent("kg/m2/s", "kg/m2")) + def test_missing_units_match_nothing(self): """ Check that an absent unit is not equivalent to anything. diff --git a/scripts/plotting/global_latlon_vect_map.py b/scripts/plotting/global_latlon_vect_map.py index eb5457c1e..c38004aa4 100644 --- a/scripts/plotting/global_latlon_vect_map.py +++ b/scripts/plotting/global_latlon_vect_map.py @@ -608,7 +608,7 @@ def _existing_plot_set( # End for # End for # End for - return plot_utils.first_complete_plot_set([flat, levelled]) + return plot_utils.first_complete_plot_set([levelled, flat]) ############## diff --git a/scripts/plotting/meridional_mean.py b/scripts/plotting/meridional_mean.py index 9314e91f2..9fa9df3e9 100644 --- a/scripts/plotting/meridional_mean.py +++ b/scripts/plotting/meridional_mean.py @@ -33,7 +33,8 @@ def meridional_mean(adfobj): Directly uses adfobj for the following: plot_var_list, plot_location, climo_yrs, variable_defaults, - read_config_var, get_basic_info, add_website_data, debug_log + read_config_var, get_basic_info, add_website_data, debug_log, + compare_obs Every plot this makes is named the same way whether or not the variable has a `lev` dimension, so a case whose plots are all present can be diff --git a/scripts/plotting/polar_map.py b/scripts/plotting/polar_map.py index 753ce659b..5d634b49a 100644 --- a/scripts/plotting/polar_map.py +++ b/scripts/plotting/polar_map.py @@ -354,7 +354,7 @@ def _existing_plot_set(plot_locations, case_names, var, seasons, pres_levs, plot # End for # End for # End for - return plot_utils.first_complete_plot_set([flat, levelled]) + return plot_utils.first_complete_plot_set([levelled, flat]) ############## diff --git a/scripts/plotting/tem.py b/scripts/plotting/tem.py index 63404f56a..a72c153ec 100644 --- a/scripts/plotting/tem.py +++ b/scripts/plotting/tem.py @@ -413,12 +413,12 @@ def tem(adf): # converted values is not scaled a second time. Observations # have their own scaling, which is assumed to bring them to the # same units, so they keep the unit string they arrived with. - mdata = adfobj.data.apply_conversion(mdata, var) + mdata = adf.data.apply_conversion(mdata, var) if obs: odata = (odata * vres.get("obs_scale_factor", 1) + vres.get("obs_add_offset", 0)) else: - odata = adfobj.data.apply_conversion(odata, var) + odata = adf.data.apply_conversion(odata, var) #End if #Month-length weighted seasonal (or annual) mean. The weighted