From b28b6ac4e8633cacb2884c27afd6345244295c04 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Thu, 26 Oct 2023 21:14:20 -0600 Subject: [PATCH 01/46] sketching out a path for AdfData; adf_dataset.py and zonal_Mean_B.py --- lib/adf_dataset.py | 162 +++++++++++++++++++ scripts/plotting/zonal_mean_B.py | 258 +++++++++++++++++++++++++++++++ 2 files changed, 420 insertions(+) create mode 100644 lib/adf_dataset.py create mode 100644 scripts/plotting/zonal_mean_B.py diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py new file mode 100644 index 000000000..4813ca734 --- /dev/null +++ b/lib/adf_dataset.py @@ -0,0 +1,162 @@ +from pathlib import Path +import xarray as xr + +import warnings # use to warn user about missing files + +def my_formatwarning(msg, *args, **kwargs): + # ignore everything except the message + return str(msg) + '\n' +warnings.formatwarning = my_formatwarning + +# "reference data" +# It is often just a "baseline case", +# but could be some totally external data (reanalysis or observation or other model) +# When it is another simulation, we treat it basically like another "case" +# When it is external data expect: +# - "climo" files (12 monthly climos in the file) +# - one variable per "climo" +# - source can differ for each variable, requires label +# - resolution can differ for each variable, requires regridded file(s) +# - the variable name and units in the file may differ from CAM; use defaults.yaml to set conversion +# - there could be multiple instances of a variable from different sources (e.g. different observations) + +# NOTE: the last item (multiple instances of a variable) is not allowed in AdfObs.var_obs_dict +# Since we are not able to handle this case, for now we exclude it from the AdfData class. + +# NOTE: To make the "baseline case" vs "external data" cases as similar as possible, +# below construct the "baseline case" version to be similar to "external data" +# To do that, need to provide a dictionary of variables as keys with paths to the +# file(s) as values. For external data, we get that from AdfObs.var_obs_dict, +# which provides a dict of all the available variables. +# When we have a reference simulation, we just need to look for +# files that match the diag_var_list + +# NOTE: There is currently a "base_nickname" allowed from AdfInfo. I will set AdfData.ref_nickname to that. +# Could be altered from "Obs" to be the data source label. + +class AdfData: + """A class instantiated with an AdfDiag object. + Methods provide means to load data. + This class does not interact with plotting, + just provides access to data locations and loading data. + + We will probably need to add some kind of frequency/sampling + parameters to allow for non-h0 files. + + """ + def __init__(self, adfobj): + self.adf = adfobj # provides quick access to the AdfDiag object + # paths + self.model_rgrid_loc = adfobj.get_basic_info("cam_regrid_loc", required=True) + + # variables (and info for unit transform) + self.var_list = adfobj.diag_var_list + self.res = adfobj.variable_defaults + + # case names and nicknames + self.case_names = adfobj.get_cam_info("cam_case_name", required=True) + self.test_nicknames = adfobj.case_nicknames["test_nicknames"] + self.base_nickname = adfobj.case_nicknames["base_nickname"] + self.reference_is_obs = adfobj.get_basic_info("compare_obs") + self.set_reference() # specify "ref_labels" -> called "data_list" in zonal_mean (name of data source) + self.ref_nickname = self.base_nickname + + def set_reference(self): + """Set attributes for reference (aka baseline) data location, names, and variables.""" + if self.reference_is_obs: + self.ref_var_loc = {v: self.adf.var_obs_dict[v]['obs_file'] for v in self.adf.var_obs_dict} + self.ref_labels = {v: self.adf.var_obs_dict[v]['obs_name'] for v in self.adf.var_obs_dict} + self.ref_var_nam = {v: self.adf.var_obs_dict[v]['obs_var'] for v in self.adf.var_obs_dict} + if not self.adf.var_obs_dict: + print("\t WARNING: reference is observations, but no observations found to plot against.") + else: + self.ref_var_loc = {} + self.ref_var_nam = {} + self.ref_labels = {} + for v in self.var_list: + f = self.get_reference_climo_file(v) + if f is None: + continue + else: + self.ref_var_loc[v] = f + self.ref_var_nam[v] = v + self.ref_labels[v] = self.adf.get_baseline_info("cam_case_name", required=True) + + def get_reference_climo_file(self, var): + """Return a list of files to be used as reference (aka baseline) for variable var.""" + if self.reference_is_obs: + return [self.ref_var_loc[var]] + else: + fils = sorted(Path(self.ref_loc).glob(f"{self.ref_label}_{var}_baseline.nc")) + if fils: + return fils + else: + return None + + def load_reference_dataset(self, var): + fils = self.get_reference_climo_file(var) + if not fils: + print(f"WARNING: Did not find any reference files for variable: {var}. Will try to skip.") + return None + return self.load_dataset(fils) + + def load_reference_da(self, variablename): + da = self.load_reference_dataset(variablename)[self.ref_var_nam[variablename]] + if variablename in self.res: + vres = self.res[variablename] + if self.reference_is_obs: + scale_factor = vres.get("obs_scale_factor",1) + add_offset = vres.get("obs_add_offset", 0) + else: + scale_factor = vres.get("scale_factor",1) + add_offset = vres.get("add_offset", 0) + da = da * scale_factor + add_offset + da.attrs['units'] = vres.get("new_unit", da.attrs.get('units', 'none')) + + def get_climo_file(self, variablename): + pass + + def get_timeseries_file(self, variablename): + pass + + def get_regrid_file(self, case, field): + model_rg_loc = Path(self.adf.get_basic_info("cam_regrid_loc", required=True)) + rlbl = self.ref_labels[field] # rlbl = "reference label" = the name of the reference data that defines target grid + return sorted(model_rg_loc.glob(f"{rlbl}_{case}_{field}_*.nc")) + + def load_regrid_dataset(self, case, field): + fils = self.get_regrid_file(case, field) + if not fils: + print(f"ERROR: Did not find regrid file(s) for case: {case}, variable: {field}") + return None + return self.load_dataset(fils) + + def load_regrid_da(self, case, field): + fils = self.get_regrid_file(case, field) + if not fils: + print(f"ERROR: Did not find regrid file(s) for case: {case}, variable: {field}") + return None + return self.load_da(fils, field) + + def get_file_list(): + pass + + def load_dataset(self, fils): + if (len(fils) == 0): + warnings.warn("Input file list is empty.") + return None + elif (len(fils) > 1): + ds = xr.open_mfdataset(fils, combine='by_coords') + else: + sfil = str(fils[0]) + ds = xr.open_dataset(sfil) + if ds is None: + warnings.warn(f"invalid data on load_dataset") + return ds + + def load_da(self, fils, variablename): + da = (self.load_dataset(fils)[variablename]).squeeze() + if variablename in self.res: + vres = self.res[variablename] + da = da * vres.get("scale_factor",1) + vres.get("add_offset", 0) + da.attrs['units'] = vres.get("new_unit", da.attrs.get('units', 'none')) diff --git a/scripts/plotting/zonal_mean_B.py b/scripts/plotting/zonal_mean_B.py new file mode 100644 index 000000000..c73f6e222 --- /dev/null +++ b/scripts/plotting/zonal_mean_B.py @@ -0,0 +1,258 @@ +from pathlib import Path +import numpy as np +import xarray as xr +import plotting_functions as pf +import warnings # use to warn user about missing files. + +def my_formatwarning(msg, *args, **kwargs): + # ignore everything except the message + return str(msg) + '\n' + +warnings.formatwarning = my_formatwarning + +def zonal_mean(adfobj): + + """ + This script plots zonal averages. + Compare CAM climatologies against + other climatological data (observations or baseline runs). + Description of needed inputs from ADF: + case_name -> Name of CAM case provided by "cam_case_name". + model_rgrid_loc -> Location of re-gridded CAM climo files provided by "cam_regrid_loc". + data_name -> Name of data set CAM case is being compared against, + which is always either "obs" or the baseline CAM case name, + depending on whether "compare_obs" is true or false. + data_loc -> Location of comparison data, which is either "obs_climo_loc" + or "cam_baseline_climo_loc", depending on whether + "compare_obs" is true or false. + var_list -> List of CAM output variables provided by "diag_var_list" + data_list -> List of data sets CAM will be compared against, which + is simply the baseline case name in situations when + "compare_obs" is false. + climo_yrs -> Dictionary containing the start and end years of the test + and baseline model data (if applicable). + plot_location -> Location where plot files will be written to, which is + specified by "cam_diag_plot_loc". + variable_defaults -> optional, + Dict that has keys that are variable names and values that are plotting preferences/defaults. + Notes: + The script produces plots of 2-D and 3-D variables, + but needs to determine which type along the way. + For 3-D variables, the default behavior is to interpolate + climo files to pressure levels, which requires the hybrid-sigma + coefficients and surface pressure. That ASSUMES that the climo + files are using native hybrid-sigma levels rather than being + transformed to pressure levels. + """ + + #Notify user that script has started: + print("\n Generating zonal mean plots...") + + # Start by instantiating the AdfData object + # and Extract needed quantities from ADF object: + + from adf_dataset import AdfData + data = AdfData(adfobj) + + var_list = adfobj.diag_var_list + + #Special ADF variable which contains the output paths for + #all generated plots and tables: + plot_locations = adfobj.plot_location + + #Grab case years + syear_cases = adfobj.climo_yrs["syears"] + eyear_cases = adfobj.climo_yrs["eyears"] + + #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"] + + res = adfobj.variable_defaults # will be dict of variable-specific plot preferences + # 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 + # -- So check for it, and default to png + basic_info_dict = adfobj.read_config_var("diag_basic_info") + plot_type = basic_info_dict.get('plot_type', 'png') + print(f"\t NOTE: Plot type is set to {plot_type}") + + # check if existing plots need to be redone + redo_plot = adfobj.get_basic_info('redo_plot') + print(f"\t NOTE: redo_plot is set to {redo_plot}") + #----------------------------------------- + + + #Set seasonal ranges: + seasons = {"ANN": np.arange(1,13,1), + "DJF": [12, 1, 2], + "JJA": [6, 7, 8], + "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 attempt to skip calcs to + #speed up preformance a bit if re-running the ADF + zonal_skip = [] + logp_zonal_skip = [] + + #Loop over model cases: + for case_idx, case_name in enumerate(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: + #Check zonal log-p: + plot_name_log = plot_loc / f"{var}_{s}_Zonal_logp_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(): + logp_zonal_skip.append(plot_name_log) + #Continue to next iteration: + adfobj.add_website_data(plot_name_log, f"{var}_logp", case_name, season=s, + plot_type="Zonal", category="Log-P") + pass + + elif (redo_plot) and plot_name_log.is_file(): + plot_name_log.unlink() + #End if + + #Check regular zonal + plot_name = plot_loc / f"{var}_{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.is_file(): + zonal_skip.append(plot_name) + #Add already-existing plot to website (if enabled): + adfobj.add_website_data(plot_name, var, case_name, season=s, + plot_type="Zonal") + + 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: + if var not in data.ref_var_nam: + dmsg = f"No obs found for variable `{var}`, zonal mean plotting skipped." + adfobj.debug_log(dmsg) + continue + + #Notify user of variable being plotted: + print(f"\t - zonal mean plots for {var}") + + # 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"zonal_mean: Found variable defaults for {var}") + + else: + vres = {} + #End if + + # load reference data (observational or baseline) + odata = data.load_reference_da(var) + + #Loop over model cases: + for case_idx, case_name in enumerate(data.case_names): + + #Set case nickname: + case_nickname = data.test_nicknames[case_idx] + + #Set output plot location: + plot_loc = Path(plot_locations[case_idx]) + + # load re-gridded model files: + mdata = data.load_regrid_da(case_name, var) + + # determine whether it's 2D or 3D + # 3D triggers search for surface pressure + has_lat, has_lev = pf.zm_validate_dims(mdata) # assumes will work for both mdata & odata + + #Notify user of level dimension: + if has_lev: + print(f"\t {var} has lev dimension.") + + # + # Seasonal Averages + # + + #Create new dictionaries: + mseasons = {} + oseasons = {} + + #Loop over season dictionary: + for s in seasons: + + # time to make plot; here we'd probably loop over whatever plots we want for this variable + # I'll just call this one "Zonal_Mean" ... would this work as a pattern [operation]_[AxesDescription] ? + # NOTE: Up to this point, nothing really differs from global_latlon_map, + # so we could have made one script instead of two. + # Merging would make overall timing better because looping twice will double I/O steps. + # + plot_name = plot_loc / f"{var}_{s}_Zonal_Mean.{plot_type}" + + if plot_name not in zonal_skip: + + #Seasonal Averages + mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) + oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) + + # difference: each entry should be (lat, lon) or (plev, lat, lon) + # dseasons[s] = mseasons[s] - oseasons[s] + # difference will be calculated in plot_zonal_mean_and_save; + # because we can let any pressure-level interpolation happen there + # This could be re-visited for efficiency or improved code structure. + + #Create new plot: + pf.plot_zonal_mean_and_save(plot_name, case_nickname, data.ref_nickname, + [syear_cases[case_idx],eyear_cases[case_idx]], + [syear_baseline,eyear_baseline], + mseasons[s], oseasons[s], has_lev, log_p=False, obs=data.reference_is_obs, **vres) + + #Add plot to website (if enabled): + adfobj.add_website_data(plot_name, var, case_name, season=s, plot_type="Zonal") + + #Create new plot with log-p: + if has_lev: + plot_name_log = plot_loc / f"{var}_{s}_Zonal_logp_Mean.{plot_type}" + + if plot_name_log not in logp_zonal_skip: + pf.plot_zonal_mean_and_save(plot_name_log, case_nickname, data.ref_nickname, + [syear_cases[case_idx],eyear_cases[case_idx]], + [syear_baseline,eyear_baseline], + mseasons[s], oseasons[s], has_lev, log_p=True, obs=data.reference_is_obs, **vres) + + #Add plot to website (if enabled): + adfobj.add_website_data(plot_name_log, f"{var}_logp", case_name, season=s, plot_type="Zonal", category="Log-P") + + #End for (seasons loop) + #End for (case names loop) + #End for (variables loop) + + #Notify user that script has ended: + print(" ...Zonal mean plots have been generated successfully.") + + +############## +#END OF SCRIPT + From 796a8c021fdb8416276fece708274c5264fcfbd5 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 27 Oct 2023 09:20:46 -0600 Subject: [PATCH 02/46] working on making it run --- scripts/averaging/create_climo_files.py | 15 +++++++++++++-- scripts/plotting/zonal_mean_B.py | 2 +- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/scripts/averaging/create_climo_files.py b/scripts/averaging/create_climo_files.py index 776c3846e..78c6d257b 100644 --- a/scripts/averaging/create_climo_files.py +++ b/scripts/averaging/create_climo_files.py @@ -8,11 +8,20 @@ def my_formatwarning(msg, *args, **kwargs): return str(msg) + '\n' warnings.formatwarning = my_formatwarning - +import numpy as np import xarray as xr # module-level import so all functions can get to it. import multiprocessing as mp +def get_time_slice_by_year(time, startyear, endyear): + if not hasattr(time, 'dt'): + print("Warning: get_time_slice_by_year requires the `time` parameter to be an xarray time coordinate with a dt accessor. Returning generic slice (which will probably fail).") + return slice(startyear, endyear) + start_time_index = np.argwhere((time.dt.year >= startyear).values).flatten().min() + end_time_index = np.argwhere((time.dt.year <= endyear).values).flatten().max() + return slice(start_time_index, end_time_index+1) + + ############## #Main function @@ -205,7 +214,9 @@ def process_variable(ts_files, syr, eyr, output_file): cam_ts_data.assign_coords(time=time) cam_ts_data = xr.decode_cf(cam_ts_data) #Extract data subset using provided year bounds: - cam_ts_data = cam_ts_data.sel(time=slice(syr, eyr)) + #cam_ts_data = cam_ts_data.sel(time=slice(syr, eyr)) + tslice = get_time_slice_by_year(cam_ts_data.time, int(syr), int(eyr)) + cam_ts_data = cam_ts_data.isel(time=tslice) #Group time series values by month, and average those months together: cam_climo_data = cam_ts_data.groupby('time.month').mean(dim='time') #Rename "months" to "time": diff --git a/scripts/plotting/zonal_mean_B.py b/scripts/plotting/zonal_mean_B.py index c73f6e222..eef15bde8 100644 --- a/scripts/plotting/zonal_mean_B.py +++ b/scripts/plotting/zonal_mean_B.py @@ -10,7 +10,7 @@ def my_formatwarning(msg, *args, **kwargs): warnings.formatwarning = my_formatwarning -def zonal_mean(adfobj): +def zonal_mean_B(adfobj): """ This script plots zonal averages. From 18706b3fe429f44b4ae4c0ce881dd60acbde8e86 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 27 Oct 2023 09:21:50 -0600 Subject: [PATCH 03/46] fix getting reference case path --- lib/adf_dataset.py | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py index 4813ca734..5a3585feb 100644 --- a/lib/adf_dataset.py +++ b/lib/adf_dataset.py @@ -87,6 +87,7 @@ def get_reference_climo_file(self, var): if self.reference_is_obs: return [self.ref_var_loc[var]] else: + self.ref_loc = self.adf.get_baseline_info("cam_climo_loc") fils = sorted(Path(self.ref_loc).glob(f"{self.ref_label}_{var}_baseline.nc")) if fils: return fils From efdfba5b6d007fd1da02714060aaee745ea50cb2 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 27 Oct 2023 09:34:20 -0600 Subject: [PATCH 04/46] allow reference case to have a label --- lib/adf_dataset.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py index 5a3585feb..251ceadde 100644 --- a/lib/adf_dataset.py +++ b/lib/adf_dataset.py @@ -57,9 +57,11 @@ def __init__(self, adfobj): self.case_names = adfobj.get_cam_info("cam_case_name", required=True) self.test_nicknames = adfobj.case_nicknames["test_nicknames"] self.base_nickname = adfobj.case_nicknames["base_nickname"] + self.ref_nickname = self.base_nickname + + # define reference data self.reference_is_obs = adfobj.get_basic_info("compare_obs") self.set_reference() # specify "ref_labels" -> called "data_list" in zonal_mean (name of data source) - self.ref_nickname = self.base_nickname def set_reference(self): """Set attributes for reference (aka baseline) data location, names, and variables.""" @@ -73,6 +75,8 @@ def set_reference(self): self.ref_var_loc = {} self.ref_var_nam = {} self.ref_labels = {} + # when using a reference simulation, allow a "special" attribute with the case name: + self.ref_case_label = self.adf.get_baseline_info("cam_case_name", required=True) for v in self.var_list: f = self.get_reference_climo_file(v) if f is None: @@ -88,7 +92,7 @@ def get_reference_climo_file(self, var): return [self.ref_var_loc[var]] else: self.ref_loc = self.adf.get_baseline_info("cam_climo_loc") - fils = sorted(Path(self.ref_loc).glob(f"{self.ref_label}_{var}_baseline.nc")) + fils = sorted(Path(self.ref_loc).glob(f"{self.ref_case_label}_{var}_baseline.nc")) if fils: return fils else: From 6456fdc853a45378c5631a87b71a0814fef52353 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Thu, 30 Nov 2023 17:40:49 -0700 Subject: [PATCH 05/46] adfdata applied to global maps script --- lib/adf_dataset.py | 41 +++- scripts/plotting/global_latlon_map_B.py | 312 ++++++++++++++++++++++++ 2 files changed, 340 insertions(+), 13 deletions(-) create mode 100644 scripts/plotting/global_latlon_map_B.py diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py index 251ceadde..a51aba579 100644 --- a/lib/adf_dataset.py +++ b/lib/adf_dataset.py @@ -11,7 +11,7 @@ def my_formatwarning(msg, *args, **kwargs): # "reference data" # It is often just a "baseline case", # but could be some totally external data (reanalysis or observation or other model) -# When it is another simulation, we treat it basically like another "case" +# When it is another simulation, it gets treated like another "case" # When it is external data expect: # - "climo" files (12 monthly climos in the file) # - one variable per "climo" @@ -21,17 +21,17 @@ def my_formatwarning(msg, *args, **kwargs): # - there could be multiple instances of a variable from different sources (e.g. different observations) # NOTE: the last item (multiple instances of a variable) is not allowed in AdfObs.var_obs_dict -# Since we are not able to handle this case, for now we exclude it from the AdfData class. +# Since ADF is not able to handle this case, for now it is excluded the AdfData class. # NOTE: To make the "baseline case" vs "external data" cases as similar as possible, -# below construct the "baseline case" version to be similar to "external data" -# To do that, need to provide a dictionary of variables as keys with paths to the -# file(s) as values. For external data, we get that from AdfObs.var_obs_dict, -# which provides a dict of all the available variables. -# When we have a reference simulation, we just need to look for -# files that match the diag_var_list - -# NOTE: There is currently a "base_nickname" allowed from AdfInfo. I will set AdfData.ref_nickname to that. +# below construct the "baseline case" version to be similar to "external data". +# - provide a dictionary of (variable: file-path) +# + For external data, that dictionay is from AdfObs.var_obs_dict, +# which provides a dict of all the available variables. +# + For reference simulation, look for files that match the diag_var_list + +# NOTE: There is currently a "base_nickname" allowed from AdfInfo. +# Set AdfData.ref_nickname to that. # Could be altered from "Obs" to be the data source label. class AdfData: @@ -40,7 +40,7 @@ class AdfData: This class does not interact with plotting, just provides access to data locations and loading data. - We will probably need to add some kind of frequency/sampling + A future need is to add some kind of frequency/sampling parameters to allow for non-h0 files. """ @@ -121,8 +121,23 @@ def load_reference_da(self, variablename): def get_climo_file(self, variablename): pass - def get_timeseries_file(self, variablename): - pass + def get_timeseries_file(self, case, field): + ts_locs = self.adf.get_cam_info("cam_ts_loc", required=True) + ts_loc = Path(ts_locs[case]) + ts_filenames = f'{case}.*.{field}.*nc' + ts_files = sorted(ts_loc.glob(ts_filenames)) + return ts_files + + def get_ref_timeseries_file(self, field): + if self.reference_is_obs: + print("Where do we get observation time series?") + return None + else: + ts_loc = Path(self.adf.get_baseline_info("cam_ts_loc", required=True)) + ts_filenames = f'{self.ref_case_label}.*.{field}.*nc' + ts_files = sorted(ts_loc.glob(ts_filenames)) + return ts_files + def get_regrid_file(self, case, field): model_rg_loc = Path(self.adf.get_basic_info("cam_regrid_loc", required=True)) diff --git a/scripts/plotting/global_latlon_map_B.py b/scripts/plotting/global_latlon_map_B.py new file mode 100644 index 000000000..2b2daa9ac --- /dev/null +++ b/scripts/plotting/global_latlon_map_B.py @@ -0,0 +1,312 @@ +""" +Generate global maps of 2-D fields + +Functions +--------- +global_latlon_map(adfobj) + use ADF object to make maps +my_formatwarning(msg, *args, **kwargs) + format warning messages + (private method) +_load_dataset(fils) + load files into dataset + (private methd) +""" +#Import standard modules: +from pathlib import Path +import numpy as np +import xarray as xr +import warnings # use to warn user about missing files. + +import plotting_functions as pf +from adf_dataset import AdfData + +#Format warning messages: +def my_formatwarning(msg, *args, **kwargs): + """Issue `msg` as warning.""" + return str(msg) + '\n' +warnings.formatwarning = my_formatwarning + +def global_latlon_map(adfobj): + """ + This script/function is designed to generate global + 2-D lat/lon maps of model fields with continental overlays. + + Parameters + ---------- + adfobj : AdfDiag + The diagnostics object that contains all the configuration information + + Returns + ------- + Does not return a value; produces plots and saves files. + + Notes + ----- + This function imports `pandas` and `plotting_functions` + + It uses the AdfDiag object's methods to get necessary information. + Specificially: + adfobj.diag_var_list + List of variables + adfobj.get_basic_info + Regrid data path, checks `compare_obs`, checks `redo_plot`, checks `plot_press_levels` + adfobj.plot_location + output plot path + adfobj.get_cam_info + Get `cam_case_name` and `case_nickname` + adfobj.climo_yrs + start and end climo years of the case(s), `syears` & `eyears` + start and end climo years of the reference, `syear_baseline` & `eyear_baseline` + adfobj.var_obs_dict + reference data (conditional) + adfobj.get_baseline_info + get reference case, `cam_case_name` + adfobj.variable_defaults + dict of variable-specific plot preferences + adfobj.read_config_var + dict of basic info, `diag_basic_info` + Then use to check `plot_type` + adfobj.compare_obs + Used to set data path + adfobj.debug_log + Issues debug message + adfobj.add_website_data + Communicates information to the website generator + + + The `plotting_functions` module is needed for: + pf.get_central_longitude() + determine central longitude for global plots + pf.lat_lon_validate_dims() + makes sure latitude and longitude are valid + pf.seasonal_mean() + calculate seasonal mean + pf.plot_map_and_save() + send information to make the plot and save the file + pf.zm_validate_dims() + Checks on pressure level dimension + """ + + #Notify user that script has started: + print("\n Generating lat/lon maps...") + + # + # Use ADF api to get all necessary information + # + data = AdfData(adfobj) + var_list = adfobj.diag_var_list + + #Special ADF variable which contains the output paths for + #all generated plots and tables for each case: + plot_locations = adfobj.plot_location + + #Grab case years + syear_cases = adfobj.climo_yrs["syears"] + eyear_cases = adfobj.climo_yrs["eyears"] + + #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"] + + res = adfobj.variable_defaults # will be dict of variable-specific plot preferences + # 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 + # -- So check for it, and default to png + basic_info_dict = adfobj.read_config_var("diag_basic_info") + plot_type = basic_info_dict.get('plot_type', 'png') + print(f"\t NOTE: Plot type is set to {plot_type}") + + # check if existing plots need to be redone + redo_plot = adfobj.get_basic_info('redo_plot') + print(f"\t NOTE: redo_plot is set to {redo_plot}") + #----------------------------------------- + + #Determine if user wants to plot 3-D variables on + #pressure levels: + pres_levs = adfobj.get_basic_info("plot_press_levels") + + weight_season = True #always do seasonal weighting + + #Set seasonal ranges: + seasons = {"ANN": np.arange(1,13,1), + "DJF": [12, 1, 2], + "JJA": [6, 7, 8], + "MAM": [3, 4, 5], + "SON": [9, 10, 11] + } + + # probably want to do this one variable at a time: + for var in var_list: + if var not in data.ref_var_nam: + dmsg = f"No reference data found for variable `{var}`, zonal mean plotting skipped." + adfobj.debug_log(dmsg) + continue + + #Notify user of variable being plotted: + print("\t - lat/lon maps for {}".format(var)) + + # 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"global_latlon_map: Found variable defaults for {var}") + + #Extract category (if available): + web_category = vres.get("category", None) + + else: + vres = {} + web_category = None + #End if + + # For global maps, also set the central longitude: + # can be specified in adfobj basic info as 'central_longitude' or supplied as a number, + # otherwise defaults to 180 + vres['central_longitude'] = pf.get_central_longitude(adfobj) + + # load reference data (observational or baseline) + odata = data.load_reference_da(var) + if odata is None: + continue + has_dims = pf.lat_lon_validate_dims(odata) # T iff dims are (lat,lon) -- can't plot unless we have both + if not has_dims: + print(f"\t = skipping global map for {var} as REFERENCE does not have both lat and lon") + continue + + #Loop over model cases: + for case_idx, case_name in enumerate(data.case_names): + + #Set case nickname: + case_nickname = data.test_nicknames[case_idx] + + #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(" {} not found, making new directory".format(plot_loc)) + plot_loc.mkdir(parents=True) + + #Load re-gridded model files: + mdata = data.load_regrid_da(case_name, var) + + #Skip this variable/case if the regridded climo file doesn't exist: + if mdata is None: + continue + + #Determine dimensions of variable: + has_dims_cam = pf.lat_lon_validate_dims(mdata) # T iff dims are (lat,lon) -- can't plot unless we have both + _, has_lev = pf.zm_validate_dims(mdata) # has_lev T if lev in mdata + if not has_dims_cam: + print(f"\t = skipping global map for {var} for case {case_name} as it does not have both lat and lon") + continue + else: # i.e., has lat&lon + if pres_levs and (not has_lev): + print(f"\t - skipping global map for {var} as it has more than lat/lon dims, but no pressure levels were provided") + continue + + # Check output file. If file does not exist, proceed. + # If file exists: + # if redo_plot is true: delete it now and make plot + # if redo_plot is false: add to website and move on + doplot = {} + if not pres_levs: + for s in seasons: + plot_name = plot_loc / f"{var}_{s}_LatLon_Mean.{plot_type}" + doplot[plot_name] = plot_file_op(adfobj, plot_name, var, case_name, s, web_category, redo_plot, "LatLon") + else: + for pres in pres_levs: + for s in seasons: + plot_name = plot_loc / f"{var}_{pres}hpa_{s}_LatLon_Mean.{plot_type}" + doplot[plot_name] = plot_file_op(adfobj, plot_name, f"{var}_{pres}hpa", case_name, s, web_category, redo_plot, "LatLon") + if all(value is None for value in doplot.values()): + print(f"All plots exist for {var}. Redo is {redo_plot}. Existing plots added to website data. Continue.") + continue + + #Create new dictionaries: + mseasons = {} + oseasons = {} + dseasons = {} # hold the differences + + if weight_season: + mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) + oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) + else: + #Just average months as-is: + mseasons[s] = mdata.sel(time=seasons[s]).mean(dim='time') + oseasons[s] = odata.sel(time=seasons[s]).mean(dim='time') + #End if + + # difference: each entry should be (lat, lon) + dseasons[s] = mseasons[s] - oseasons[s] + + if not pres_levs: + + #Loop over season dictionary: + for s in seasons: + plot_name = plot_loc / f"{var}_{s}_LatLon_Mean.{plot_type}" + if doplot[plot_name] is None: + continue + pf.plot_map_and_save(plot_name, case_nickname, data.ref_nickname, + [syear_cases[case_idx],eyear_cases[case_idx]], + [syear_baseline,eyear_baseline], + mseasons[s], oseasons[s], dseasons[s], + obs=data.reference_is_obs, **vres) + + #Add plot to website (if enabled): + adfobj.add_website_data(plot_name, var, case_name, category=web_category, + season=s, plot_type="LatLon") + + else: # => pres_levs has values, & we already checked that lev is in mdata (has_lev) + + for pres in pres_levs: + + #Check that the user-requested pressure level + #exists in the model data, which should already + #have been interpolated to the standard reference + #pressure levels: + if not (pres in mdata['lev']): + print(f"plot_press_levels value '{pres}' not present in {var}, so skipping.") + continue + + #Loop over seasons: + for s in seasons: + plot_name = plot_loc / f"{var}_{pres}hpa_{s}_LatLon_Mean.{plot_type}" + if doplot[plot_name] is None: + continue + pf.plot_map_and_save(plot_name, case_nickname, data.ref_nickname, + [syear_cases[case_idx],eyear_cases[case_idx]], + [syear_baseline,eyear_baseline], + mseasons[s].sel(lev=pres), oseasons[s].sel(lev=pres), dseasons[s].sel(lev=pres), + obs=data.reference_is_obs, **vres) + + #Add plot to website (if enabled): + adfobj.add_website_data(plot_name, f"{var}_{pres}hpa", case_name, category=web_category, + season=s, plot_type="LatLon") + #End for (seasons) + #End for (pressure levels) + #End if (plotting pressure levels) + #End for (case loop) + #End for (variable loop) + + #Notify user that script has ended: + print(" ...lat/lon maps have been generated successfully.") + + +def plot_file_op(adfobj, plot_name, var, case_name, season, web_category, redo_plot, plot_type): + """Check if output plot needs to be made or remade.""" + # 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.add_website_data(plot_name, var, case_name, category=web_category, + season=season, plot_type=plot_type) + return None # None tells caller that file exists and not to overwrite + elif (redo_plot) and plot_name.is_file(): + plot_name.unlink() + return 1 + +############## +#END OF SCRIPT From aababb277cca88ddbf0998c8560fd16ae3cf0bf8 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Thu, 30 Nov 2023 17:43:34 -0700 Subject: [PATCH 06/46] adfdata applied to global maps script -- rename function --- scripts/plotting/global_latlon_map_B.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/plotting/global_latlon_map_B.py b/scripts/plotting/global_latlon_map_B.py index 2b2daa9ac..4d8b7e0fd 100644 --- a/scripts/plotting/global_latlon_map_B.py +++ b/scripts/plotting/global_latlon_map_B.py @@ -27,7 +27,7 @@ def my_formatwarning(msg, *args, **kwargs): return str(msg) + '\n' warnings.formatwarning = my_formatwarning -def global_latlon_map(adfobj): +def global_latlon_map_B(adfobj): """ This script/function is designed to generate global 2-D lat/lon maps of model fields with continental overlays. From 988011b945487e25256a3db5345e880dbd815c24 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Thu, 30 Nov 2023 17:46:51 -0700 Subject: [PATCH 07/46] adfdata applied to global maps script -- debugging --- scripts/plotting/global_latlon_map_B.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/plotting/global_latlon_map_B.py b/scripts/plotting/global_latlon_map_B.py index 4d8b7e0fd..a6249d3ea 100644 --- a/scripts/plotting/global_latlon_map_B.py +++ b/scripts/plotting/global_latlon_map_B.py @@ -96,7 +96,7 @@ def global_latlon_map_B(adfobj): # data = AdfData(adfobj) var_list = adfobj.diag_var_list - + print(f"DEBUG: var_list = {var_list}") #Special ADF variable which contains the output paths for #all generated plots and tables for each case: plot_locations = adfobj.plot_location From 91bce2a5eb1efbce6a499f5cc960ca95d48f4827 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Thu, 30 Nov 2023 17:49:00 -0700 Subject: [PATCH 08/46] adfdata applied to global maps script -- debugging --- scripts/plotting/global_latlon_map_B.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/scripts/plotting/global_latlon_map_B.py b/scripts/plotting/global_latlon_map_B.py index a6249d3ea..8e2cb7dc5 100644 --- a/scripts/plotting/global_latlon_map_B.py +++ b/scripts/plotting/global_latlon_map_B.py @@ -140,9 +140,12 @@ def global_latlon_map_B(adfobj): # probably want to do this one variable at a time: for var in var_list: + print(f"DEBUG -- {var = }") if var not in data.ref_var_nam: dmsg = f"No reference data found for variable `{var}`, zonal mean plotting skipped." adfobj.debug_log(dmsg) + print(dmsg) + print(f"DEBUG: {data.ref_var_nam = }") continue #Notify user of variable being plotted: From bf02a1f414f44db7c16d13abd4cf2dd58ce27a2b Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Thu, 30 Nov 2023 17:52:01 -0700 Subject: [PATCH 09/46] adfdata applied to global maps script -- debugging --- lib/adf_dataset.py | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py index a51aba579..b28886920 100644 --- a/lib/adf_dataset.py +++ b/lib/adf_dataset.py @@ -80,6 +80,7 @@ def set_reference(self): for v in self.var_list: f = self.get_reference_climo_file(v) if f is None: + print(f"DEBUG ADFDATA -- no ref climo file for {v}") continue else: self.ref_var_loc[v] = f From 71578ef582ce4964f357179ec440a06f604852a8 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 09:01:14 -0700 Subject: [PATCH 10/46] debugging --- scripts/plotting/global_latlon_map_B.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/scripts/plotting/global_latlon_map_B.py b/scripts/plotting/global_latlon_map_B.py index 8e2cb7dc5..25a809713 100644 --- a/scripts/plotting/global_latlon_map_B.py +++ b/scripts/plotting/global_latlon_map_B.py @@ -145,7 +145,8 @@ def global_latlon_map_B(adfobj): dmsg = f"No reference data found for variable `{var}`, zonal mean plotting skipped." adfobj.debug_log(dmsg) print(dmsg) - print(f"DEBUG: {data.ref_var_nam = }") + print(f"DEBUG: {data.ref_var_nam = } [check if obs: {data.reference_is_obs}]") + print(f"DEBUG2: {data.ref_var_loc}") continue #Notify user of variable being plotted: From 06d3ea17fbfc5495898db3e2a0e91a6f962f3d43 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 09:09:20 -0700 Subject: [PATCH 11/46] debugging --- lib/adf_dataset.py | 1 + scripts/plotting/global_latlon_map_B.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py index b28886920..09aa0b1bc 100644 --- a/lib/adf_dataset.py +++ b/lib/adf_dataset.py @@ -93,6 +93,7 @@ def get_reference_climo_file(self, var): return [self.ref_var_loc[var]] else: self.ref_loc = self.adf.get_baseline_info("cam_climo_loc") + print(f"DEBGUG/get_reference_climo_file: {self.ref_loc}") fils = sorted(Path(self.ref_loc).glob(f"{self.ref_case_label}_{var}_baseline.nc")) if fils: return fils diff --git a/scripts/plotting/global_latlon_map_B.py b/scripts/plotting/global_latlon_map_B.py index 25a809713..c2d2b2d4c 100644 --- a/scripts/plotting/global_latlon_map_B.py +++ b/scripts/plotting/global_latlon_map_B.py @@ -146,7 +146,7 @@ def global_latlon_map_B(adfobj): adfobj.debug_log(dmsg) print(dmsg) print(f"DEBUG: {data.ref_var_nam = } [check if obs: {data.reference_is_obs}]") - print(f"DEBUG2: {data.ref_var_loc}") + print(f"DEBUG2: Label: {data.ref_case_label}") continue #Notify user of variable being plotted: From 11f1477b7506bb5dec8b0e990fffeff397c0560a Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 09:21:50 -0700 Subject: [PATCH 12/46] debugging --- lib/adf_dataset.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py index 09aa0b1bc..85f83fa55 100644 --- a/lib/adf_dataset.py +++ b/lib/adf_dataset.py @@ -94,7 +94,8 @@ def get_reference_climo_file(self, var): else: self.ref_loc = self.adf.get_baseline_info("cam_climo_loc") print(f"DEBGUG/get_reference_climo_file: {self.ref_loc}") - fils = sorted(Path(self.ref_loc).glob(f"{self.ref_case_label}_{var}_baseline.nc")) + # NOTE: originally had this looking for *_baseline.nc + fils = sorted(Path(self.ref_loc).glob(f"{self.ref_case_label}_{var}_climo.nc")) if fils: return fils else: From c7286b74462ab96b4da442aff0d0281166aab9f2 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 09:26:12 -0700 Subject: [PATCH 13/46] debugging --- lib/adf_dataset.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py index 85f83fa55..4d4a7230c 100644 --- a/lib/adf_dataset.py +++ b/lib/adf_dataset.py @@ -159,6 +159,8 @@ def load_regrid_da(self, case, field): if not fils: print(f"ERROR: Did not find regrid file(s) for case: {case}, variable: {field}") return None + else: + print(f"DEBUG/load_regrid_da: {fils}") return self.load_da(fils, field) def get_file_list(): From 97e1cbdff73eba78d86eddec4ce72680164c19d9 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 09:28:19 -0700 Subject: [PATCH 14/46] debugging --- scripts/plotting/zonal_mean_B.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/scripts/plotting/zonal_mean_B.py b/scripts/plotting/zonal_mean_B.py index eef15bde8..4ed41b424 100644 --- a/scripts/plotting/zonal_mean_B.py +++ b/scripts/plotting/zonal_mean_B.py @@ -183,6 +183,8 @@ def zonal_mean_B(adfobj): # load re-gridded model files: mdata = data.load_regrid_da(case_name, var) + if mdata is None: + print(f"MAJOR PROBLEM -- mdata is none -- {case_name = }, {var = }") # determine whether it's 2D or 3D # 3D triggers search for surface pressure From 9068ec901f0e8410f7d3bedb6dd081be6bd16321 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 09:32:06 -0700 Subject: [PATCH 15/46] debugging -- fix load_da; actually return the DataArray --- lib/adf_dataset.py | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py index 4d4a7230c..72b4738d5 100644 --- a/lib/adf_dataset.py +++ b/lib/adf_dataset.py @@ -185,3 +185,4 @@ def load_da(self, fils, variablename): vres = self.res[variablename] da = da * vres.get("scale_factor",1) + vres.get("add_offset", 0) da.attrs['units'] = vres.get("new_unit", da.attrs.get('units', 'none')) + return da From ceecec9ab7e53ec7ca22ad48aef319b1b3e8cbee Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 09:33:26 -0700 Subject: [PATCH 16/46] debugging --- scripts/plotting/zonal_mean_B.py | 1 + 1 file changed, 1 insertion(+) diff --git a/scripts/plotting/zonal_mean_B.py b/scripts/plotting/zonal_mean_B.py index 4ed41b424..2ee8df19f 100644 --- a/scripts/plotting/zonal_mean_B.py +++ b/scripts/plotting/zonal_mean_B.py @@ -171,6 +171,7 @@ def zonal_mean_B(adfobj): # load reference data (observational or baseline) odata = data.load_reference_da(var) + print(f"DEBUGGING ODATA: {odata}") #Loop over model cases: for case_idx, case_name in enumerate(data.case_names): From b48d01a94c07f9747f6875004c3b6d310d63bee6 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 09:34:33 -0700 Subject: [PATCH 17/46] debugging -- fix load_reference_da; actually return the DataArray --- lib/adf_dataset.py | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py index 72b4738d5..dbe92e191 100644 --- a/lib/adf_dataset.py +++ b/lib/adf_dataset.py @@ -120,6 +120,7 @@ def load_reference_da(self, variablename): add_offset = vres.get("add_offset", 0) da = da * scale_factor + add_offset da.attrs['units'] = vres.get("new_unit", da.attrs.get('units', 'none')) + return da def get_climo_file(self, variablename): pass From dd015ef5261a5c96917910369243bbfd98342359 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 09:38:07 -0700 Subject: [PATCH 18/46] debugging --- scripts/plotting/zonal_mean_B.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/scripts/plotting/zonal_mean_B.py b/scripts/plotting/zonal_mean_B.py index 2ee8df19f..76c7d9ffb 100644 --- a/scripts/plotting/zonal_mean_B.py +++ b/scripts/plotting/zonal_mean_B.py @@ -171,7 +171,6 @@ def zonal_mean_B(adfobj): # load reference data (observational or baseline) odata = data.load_reference_da(var) - print(f"DEBUGGING ODATA: {odata}") #Loop over model cases: for case_idx, case_name in enumerate(data.case_names): @@ -184,8 +183,6 @@ def zonal_mean_B(adfobj): # load re-gridded model files: mdata = data.load_regrid_da(case_name, var) - if mdata is None: - print(f"MAJOR PROBLEM -- mdata is none -- {case_name = }, {var = }") # determine whether it's 2D or 3D # 3D triggers search for surface pressure @@ -219,7 +216,10 @@ def zonal_mean_B(adfobj): #Seasonal Averages mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) - + print("CHECK mseasons[s]") + print(mseasons[s]) + print("CHECK oseasons[s]") + print(oseasons[s]) # difference: each entry should be (lat, lon) or (plev, lat, lon) # dseasons[s] = mseasons[s] - oseasons[s] # difference will be calculated in plot_zonal_mean_and_save; From 84150a4905cb4e87156fb14f7562139661255e02 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 09:41:44 -0700 Subject: [PATCH 19/46] debugging --- lib/plotting_functions.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/lib/plotting_functions.py b/lib/plotting_functions.py index fa8bcc2a7..b468740da 100644 --- a/lib/plotting_functions.py +++ b/lib/plotting_functions.py @@ -1784,7 +1784,10 @@ def prep_contour_plot(adata, bdata, diffdata, **kwargs): levelsdiff = np.arange(*kwargs['diff_contour_range']) else: # set a symmetric color bar for diff: - absmaxdif = np.max(np.abs(diffdata)) + try: + absmaxdif = np.max(np.abs(diffdata)) + except: + print(f"WHAT IS GOING ON: {np.max(diffdata) = }, {np.min(diffdata) = }") # set levels for difference plot: levelsdiff = np.linspace(-1*absmaxdif, absmaxdif, 12) #End if From 42614fb8099e182d7c2efe38a9a0b9018ae952b0 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 09:43:30 -0700 Subject: [PATCH 20/46] debugging --- lib/plotting_functions.py | 1 + scripts/plotting/zonal_mean_B.py | 5 +---- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/lib/plotting_functions.py b/lib/plotting_functions.py index b468740da..e75158392 100644 --- a/lib/plotting_functions.py +++ b/lib/plotting_functions.py @@ -1784,6 +1784,7 @@ def prep_contour_plot(adata, bdata, diffdata, **kwargs): levelsdiff = np.arange(*kwargs['diff_contour_range']) else: # set a symmetric color bar for diff: + print(f"[pre] WHAT IS GOING ON: {np.max(diffdata) = }, {np.min(diffdata) = }") try: absmaxdif = np.max(np.abs(diffdata)) except: diff --git a/scripts/plotting/zonal_mean_B.py b/scripts/plotting/zonal_mean_B.py index 76c7d9ffb..50caebe5a 100644 --- a/scripts/plotting/zonal_mean_B.py +++ b/scripts/plotting/zonal_mean_B.py @@ -216,10 +216,7 @@ def zonal_mean_B(adfobj): #Seasonal Averages mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) - print("CHECK mseasons[s]") - print(mseasons[s]) - print("CHECK oseasons[s]") - print(oseasons[s]) + # difference: each entry should be (lat, lon) or (plev, lat, lon) # dseasons[s] = mseasons[s] - oseasons[s] # difference will be calculated in plot_zonal_mean_and_save; From b7765ab281d9f57218858a0f51d8d83b6aea9c43 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 09:45:20 -0700 Subject: [PATCH 21/46] debugging --- lib/plotting_functions.py | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/plotting_functions.py b/lib/plotting_functions.py index e75158392..0c09c9d16 100644 --- a/lib/plotting_functions.py +++ b/lib/plotting_functions.py @@ -1784,6 +1784,7 @@ def prep_contour_plot(adata, bdata, diffdata, **kwargs): levelsdiff = np.arange(*kwargs['diff_contour_range']) else: # set a symmetric color bar for diff: + print(f"[pre] {type(diffdata)}, {diffdata.shape = }, {diffdata}") print(f"[pre] WHAT IS GOING ON: {np.max(diffdata) = }, {np.min(diffdata) = }") try: absmaxdif = np.max(np.abs(diffdata)) From 857a4411b04ec4c328d8208a5f2ba9cd8554da18 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 09:49:12 -0700 Subject: [PATCH 22/46] debugging --- scripts/plotting/zonal_mean_B.py | 43 ++++++++++++++++++-------------- 1 file changed, 24 insertions(+), 19 deletions(-) diff --git a/scripts/plotting/zonal_mean_B.py b/scripts/plotting/zonal_mean_B.py index 50caebe5a..b1e7e1b49 100644 --- a/scripts/plotting/zonal_mean_B.py +++ b/scripts/plotting/zonal_mean_B.py @@ -209,34 +209,39 @@ def zonal_mean_B(adfobj): # so we could have made one script instead of two. # Merging would make overall timing better because looping twice will double I/O steps. # - plot_name = plot_loc / f"{var}_{s}_Zonal_Mean.{plot_type}" + if not has_lev: + plot_name = plot_loc / f"{var}_{s}_Zonal_Mean.{plot_type}" - if plot_name not in zonal_skip: + if plot_name not in zonal_skip: - #Seasonal Averages - mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) - oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) + #Seasonal Averages + mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) + oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) - # difference: each entry should be (lat, lon) or (plev, lat, lon) - # dseasons[s] = mseasons[s] - oseasons[s] - # difference will be calculated in plot_zonal_mean_and_save; - # because we can let any pressure-level interpolation happen there - # This could be re-visited for efficiency or improved code structure. + # difference: each entry should be (lat, lon) or (plev, lat, lon) + # dseasons[s] = mseasons[s] - oseasons[s] + # difference will be calculated in plot_zonal_mean_and_save; + # because we can let any pressure-level interpolation happen there + # This could be re-visited for efficiency or improved code structure. - #Create new plot: - pf.plot_zonal_mean_and_save(plot_name, case_nickname, data.ref_nickname, - [syear_cases[case_idx],eyear_cases[case_idx]], - [syear_baseline,eyear_baseline], - mseasons[s], oseasons[s], has_lev, log_p=False, obs=data.reference_is_obs, **vres) + #Create new plot: + pf.plot_zonal_mean_and_save(plot_name, case_nickname, data.ref_nickname, + [syear_cases[case_idx],eyear_cases[case_idx]], + [syear_baseline,eyear_baseline], + mseasons[s], oseasons[s], has_lev, log_p=False, obs=data.reference_is_obs, **vres) - #Add plot to website (if enabled): - adfobj.add_website_data(plot_name, var, case_name, season=s, plot_type="Zonal") + #Add plot to website (if enabled): + adfobj.add_website_data(plot_name, var, case_name, season=s, plot_type="Zonal") #Create new plot with log-p: - if has_lev: + # NOTE: The log-p should be an option here. + else: plot_name_log = plot_loc / f"{var}_{s}_Zonal_logp_Mean.{plot_type}" - if plot_name_log not in logp_zonal_skip: + #Seasonal Averages + mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) + oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) + pf.plot_zonal_mean_and_save(plot_name_log, case_nickname, data.ref_nickname, [syear_cases[case_idx],eyear_cases[case_idx]], [syear_baseline,eyear_baseline], From 577c9e7f68f6698eeea7a2bbaa7f6501f818601a Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 12:19:00 -0700 Subject: [PATCH 23/46] debugging --- lib/plotting_functions.py | 1 + 1 file changed, 1 insertion(+) diff --git a/lib/plotting_functions.py b/lib/plotting_functions.py index 0c09c9d16..cebc3b88d 100644 --- a/lib/plotting_functions.py +++ b/lib/plotting_functions.py @@ -1895,6 +1895,7 @@ def plot_zonal_mean_and_save(wks, case_nickname, base_nickname, # calculate difference: diff = azm - bzm + print(f"[ZONAL] {azm.shape = }, {bzm.shape = }, {diff.shape = }") # generate dictionary of contour plot settings: cp_info = prep_contour_plot(azm, bzm, diff, **kwargs) From 87aebf200c493046251246cc2ab8b9e023832d37 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 12:25:54 -0700 Subject: [PATCH 24/46] debugging --- scripts/plotting/zonal_mean_B.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/scripts/plotting/zonal_mean_B.py b/scripts/plotting/zonal_mean_B.py index b1e7e1b49..e734fcb6e 100644 --- a/scripts/plotting/zonal_mean_B.py +++ b/scripts/plotting/zonal_mean_B.py @@ -171,6 +171,7 @@ def zonal_mean_B(adfobj): # load reference data (observational or baseline) odata = data.load_reference_da(var) + has_lat_ref, has_lev_ref = pf.zm_validate_dims(mdata) #Loop over model cases: for case_idx, case_name in enumerate(data.case_names): @@ -236,6 +237,12 @@ def zonal_mean_B(adfobj): #Create new plot with log-p: # NOTE: The log-p should be an option here. else: + if (not has_lev_ref) or (not has_lev): + print(f"Error: expecting lev for both case: {has_lev} and ref: {has_lev_ref}") + continue + if len(mdata['lev']) != len(odata['lev']): + print(f"Error: zonal mean contour expects `lev` dim to have same size, got {len(mdata['lev'])} and {odata['lev']}") + continue plot_name_log = plot_loc / f"{var}_{s}_Zonal_logp_Mean.{plot_type}" if plot_name_log not in logp_zonal_skip: #Seasonal Averages From c751d026ed4ec81a567f3fb4d34692a5f32fa734 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 12:27:06 -0700 Subject: [PATCH 25/46] debugging --- scripts/plotting/zonal_mean_B.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/plotting/zonal_mean_B.py b/scripts/plotting/zonal_mean_B.py index e734fcb6e..f2387ad9a 100644 --- a/scripts/plotting/zonal_mean_B.py +++ b/scripts/plotting/zonal_mean_B.py @@ -171,7 +171,7 @@ def zonal_mean_B(adfobj): # load reference data (observational or baseline) odata = data.load_reference_da(var) - has_lat_ref, has_lev_ref = pf.zm_validate_dims(mdata) + has_lat_ref, has_lev_ref = pf.zm_validate_dims(odata) #Loop over model cases: for case_idx, case_name in enumerate(data.case_names): From f9f296d8c34cb891469ae29f69d8267aad09119e Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 13:14:40 -0700 Subject: [PATCH 26/46] completed updates for AdfData and implement in zonal mean and global maps scripts --- lib/plotting_functions.py | 12 ++--- scripts/plotting/global_latlon_map_B.py | 65 +++++++++++++++++-------- scripts/plotting/zonal_mean_B.py | 57 ++++++++++------------ 3 files changed, 73 insertions(+), 61 deletions(-) diff --git a/lib/plotting_functions.py b/lib/plotting_functions.py index cebc3b88d..d8686e691 100644 --- a/lib/plotting_functions.py +++ b/lib/plotting_functions.py @@ -1784,12 +1784,7 @@ def prep_contour_plot(adata, bdata, diffdata, **kwargs): levelsdiff = np.arange(*kwargs['diff_contour_range']) else: # set a symmetric color bar for diff: - print(f"[pre] {type(diffdata)}, {diffdata.shape = }, {diffdata}") - print(f"[pre] WHAT IS GOING ON: {np.max(diffdata) = }, {np.min(diffdata) = }") - try: - absmaxdif = np.max(np.abs(diffdata)) - except: - print(f"WHAT IS GOING ON: {np.max(diffdata) = }, {np.min(diffdata) = }") + absmaxdif = np.max(np.abs(diffdata)) # set levels for difference plot: levelsdiff = np.linspace(-1*absmaxdif, absmaxdif, 12) #End if @@ -1895,7 +1890,6 @@ def plot_zonal_mean_and_save(wks, case_nickname, base_nickname, # calculate difference: diff = azm - bzm - print(f"[ZONAL] {azm.shape = }, {bzm.shape = }, {diff.shape = }") # generate dictionary of contour plot settings: cp_info = prep_contour_plot(azm, bzm, diff, **kwargs) @@ -2271,9 +2265,9 @@ def square_contour_difference(fld1, fld2, **kwargs): mnorm = mpl.colors.Normalize(mn, mx) coord1, coord2 = fld1.coords # ASSUMES xarray WITH coords AND 2-dimensions - print(f"{coord1}, {coord2}") + # print(f"{coord1}, {coord2}") xx, yy = np.meshgrid(fld1[coord2], fld1[coord1]) - print(f"shape of meshgrid: {xx.shape}") + # print(f"shape of meshgrid: {xx.shape}") img1 = ax1.contourf(xx, yy, fld1.transpose()) if (coord1 == 'month') and (fld1.shape[0] ==12): diff --git a/scripts/plotting/global_latlon_map_B.py b/scripts/plotting/global_latlon_map_B.py index c2d2b2d4c..4e49a44c3 100644 --- a/scripts/plotting/global_latlon_map_B.py +++ b/scripts/plotting/global_latlon_map_B.py @@ -8,9 +8,8 @@ my_formatwarning(msg, *args, **kwargs) format warning messages (private method) -_load_dataset(fils) - load files into dataset - (private methd) +plot_file_op + Check on status of output plot file. """ #Import standard modules: from pathlib import Path @@ -43,32 +42,22 @@ def global_latlon_map_B(adfobj): Notes ----- - This function imports `pandas` and `plotting_functions` It uses the AdfDiag object's methods to get necessary information. - Specificially: + This version passes the AdfDiag object to instatiate an AdfData object. + Explicitly accesses: adfobj.diag_var_list List of variables - adfobj.get_basic_info - Regrid data path, checks `compare_obs`, checks `redo_plot`, checks `plot_press_levels` adfobj.plot_location output plot path - adfobj.get_cam_info - Get `cam_case_name` and `case_nickname` adfobj.climo_yrs start and end climo years of the case(s), `syears` & `eyears` start and end climo years of the reference, `syear_baseline` & `eyear_baseline` - adfobj.var_obs_dict - reference data (conditional) - adfobj.get_baseline_info - get reference case, `cam_case_name` adfobj.variable_defaults dict of variable-specific plot preferences adfobj.read_config_var dict of basic info, `diag_basic_info` Then use to check `plot_type` - adfobj.compare_obs - Used to set data path adfobj.debug_log Issues debug message adfobj.add_website_data @@ -96,7 +85,6 @@ def global_latlon_map_B(adfobj): # data = AdfData(adfobj) var_list = adfobj.diag_var_list - print(f"DEBUG: var_list = {var_list}") #Special ADF variable which contains the output paths for #all generated plots and tables for each case: plot_locations = adfobj.plot_location @@ -140,13 +128,9 @@ def global_latlon_map_B(adfobj): # probably want to do this one variable at a time: for var in var_list: - print(f"DEBUG -- {var = }") if var not in data.ref_var_nam: dmsg = f"No reference data found for variable `{var}`, zonal mean plotting skipped." adfobj.debug_log(dmsg) - print(dmsg) - print(f"DEBUG: {data.ref_var_nam = } [check if obs: {data.reference_is_obs}]") - print(f"DEBUG2: Label: {data.ref_case_label}") continue #Notify user of variable being plotted: @@ -301,7 +285,46 @@ def global_latlon_map_B(adfobj): def plot_file_op(adfobj, plot_name, var, case_name, season, web_category, redo_plot, plot_type): - """Check if output plot needs to be made or remade.""" + """Check if output plot needs to be made or remade. + + Parameters + ---------- + adfobj : AdfDiag + The diagnostics object that contains all the configuration information + + plot_name : Path + path of the output plot + + var : str + name of variable + + case_name : str + case name + + season : str + season being plotted + + web_category : str + the category for this variable + + redo_plot : bool + whether to overwrite existing plot with this file name + + plot_type : str + the file type for the output plot + + Returns + ------- + int, None + Returns 1 if file is removed + Returns None if file exists by redo_plot is False + + Notes + ----- + The long list of parameters is because add_website_data is called + when the file exists and will not be overwritten. + + """ # 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): diff --git a/scripts/plotting/zonal_mean_B.py b/scripts/plotting/zonal_mean_B.py index f2387ad9a..1374dadd8 100644 --- a/scripts/plotting/zonal_mean_B.py +++ b/scripts/plotting/zonal_mean_B.py @@ -2,6 +2,8 @@ import numpy as np import xarray as xr import plotting_functions as pf +from adf_dataset import AdfData + import warnings # use to warn user about missing files. def my_formatwarning(msg, *args, **kwargs): @@ -13,45 +15,38 @@ def my_formatwarning(msg, *args, **kwargs): def zonal_mean_B(adfobj): """ - This script plots zonal averages. + Plots zonal average from climatological files (annual and seasonal). Compare CAM climatologies against other climatological data (observations or baseline runs). - Description of needed inputs from ADF: - case_name -> Name of CAM case provided by "cam_case_name". - model_rgrid_loc -> Location of re-gridded CAM climo files provided by "cam_regrid_loc". - data_name -> Name of data set CAM case is being compared against, - which is always either "obs" or the baseline CAM case name, - depending on whether "compare_obs" is true or false. - data_loc -> Location of comparison data, which is either "obs_climo_loc" - or "cam_baseline_climo_loc", depending on whether - "compare_obs" is true or false. - var_list -> List of CAM output variables provided by "diag_var_list" - data_list -> List of data sets CAM will be compared against, which - is simply the baseline case name in situations when - "compare_obs" is false. - climo_yrs -> Dictionary containing the start and end years of the test - and baseline model data (if applicable). - plot_location -> Location where plot files will be written to, which is - specified by "cam_diag_plot_loc". - variable_defaults -> optional, - Dict that has keys that are variable names and values that are plotting preferences/defaults. - Notes: - The script produces plots of 2-D and 3-D variables, - but needs to determine which type along the way. - For 3-D variables, the default behavior is to interpolate - climo files to pressure levels, which requires the hybrid-sigma - coefficients and surface pressure. That ASSUMES that the climo - files are using native hybrid-sigma levels rather than being - transformed to pressure levels. + + 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: + diag_var_list, climo_yrs, variable_defaults, read_config_var, + get_basic_info, add_website_data, debug_log + + Determines whether `lev` dimension is present. If not, makes + a line plot, but if so it makes a contour plot. + TODO: There's a flag to plot linear vs log pressure, but no + method to infer what the user wants. """ - #Notify user that script has started: print("\n Generating zonal mean plots...") # Start by instantiating the AdfData object # and Extract needed quantities from ADF object: - - from adf_dataset import AdfData data = AdfData(adfobj) var_list = adfobj.diag_var_list From 9ec949b507e0d05a161fb2bc948e117c2d79bc8b Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 13:20:41 -0700 Subject: [PATCH 27/46] remove extra debugging print statements --- lib/adf_dataset.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py index dbe92e191..6f06bfd32 100644 --- a/lib/adf_dataset.py +++ b/lib/adf_dataset.py @@ -80,7 +80,7 @@ def set_reference(self): for v in self.var_list: f = self.get_reference_climo_file(v) if f is None: - print(f"DEBUG ADFDATA -- no ref climo file for {v}") + print(f"\t WARNING: ADFData found no reference climo file for {v}") continue else: self.ref_var_loc[v] = f @@ -93,7 +93,6 @@ def get_reference_climo_file(self, var): return [self.ref_var_loc[var]] else: self.ref_loc = self.adf.get_baseline_info("cam_climo_loc") - print(f"DEBGUG/get_reference_climo_file: {self.ref_loc}") # NOTE: originally had this looking for *_baseline.nc fils = sorted(Path(self.ref_loc).glob(f"{self.ref_case_label}_{var}_climo.nc")) if fils: @@ -104,7 +103,7 @@ def get_reference_climo_file(self, var): def load_reference_dataset(self, var): fils = self.get_reference_climo_file(var) if not fils: - print(f"WARNING: Did not find any reference files for variable: {var}. Will try to skip.") + print(f"ERROR: Did not find any reference files for variable: {var}. Will try to skip.") return None return self.load_dataset(fils) @@ -134,7 +133,6 @@ def get_timeseries_file(self, case, field): def get_ref_timeseries_file(self, field): if self.reference_is_obs: - print("Where do we get observation time series?") return None else: ts_loc = Path(self.adf.get_baseline_info("cam_ts_loc", required=True)) @@ -160,8 +158,6 @@ def load_regrid_da(self, case, field): if not fils: print(f"ERROR: Did not find regrid file(s) for case: {case}, variable: {field}") return None - else: - print(f"DEBUG/load_regrid_da: {fils}") return self.load_da(fils, field) def get_file_list(): From 3adbfe044a28fb326934fc640aadc534f7446e10 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 13:30:10 -0700 Subject: [PATCH 28/46] refactor plot_file_op for easier logic --- scripts/plotting/global_latlon_map_B.py | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/scripts/plotting/global_latlon_map_B.py b/scripts/plotting/global_latlon_map_B.py index 4e49a44c3..bf177ba29 100644 --- a/scripts/plotting/global_latlon_map_B.py +++ b/scripts/plotting/global_latlon_map_B.py @@ -316,8 +316,8 @@ def plot_file_op(adfobj, plot_name, var, case_name, season, web_category, redo_p Returns ------- int, None - Returns 1 if file is removed - Returns None if file exists by redo_plot is False + Returns 1 if existing file is removed or no existing file. + Returns None if file exists and redo_plot is False Notes ----- @@ -326,13 +326,16 @@ def plot_file_op(adfobj, plot_name, var, case_name, season, web_category, redo_p """ # 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.add_website_data(plot_name, var, case_name, category=web_category, - season=season, plot_type=plot_type) - return None # None tells caller that file exists and not to overwrite - elif (redo_plot) and plot_name.is_file(): - plot_name.unlink() + if plot_name.is_file(): + if redo_plot: + plot_name.unlink() + return 1 + else: + #Add already-existing plot to website (if enabled): + adfobj.add_website_data(plot_name, var, case_name, category=web_category, + season=season, plot_type=plot_type) + return None # None tells caller that file exists and not to overwrite + else: return 1 ############## From 51394a6bdde087d7d508f956c56aebf1cbe30f54 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 13:33:35 -0700 Subject: [PATCH 29/46] seasonal averaging moved to inside season loop --- scripts/plotting/global_latlon_map_B.py | 29 +++++++++++++++++-------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/scripts/plotting/global_latlon_map_B.py b/scripts/plotting/global_latlon_map_B.py index bf177ba29..2b48eec09 100644 --- a/scripts/plotting/global_latlon_map_B.py +++ b/scripts/plotting/global_latlon_map_B.py @@ -219,15 +219,6 @@ def global_latlon_map_B(adfobj): oseasons = {} dseasons = {} # hold the differences - if weight_season: - mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) - oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) - else: - #Just average months as-is: - mseasons[s] = mdata.sel(time=seasons[s]).mean(dim='time') - oseasons[s] = odata.sel(time=seasons[s]).mean(dim='time') - #End if - # difference: each entry should be (lat, lon) dseasons[s] = mseasons[s] - oseasons[s] @@ -238,6 +229,16 @@ def global_latlon_map_B(adfobj): plot_name = plot_loc / f"{var}_{s}_LatLon_Mean.{plot_type}" if doplot[plot_name] is None: continue + + if weight_season: + mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) + oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) + else: + #Just average months as-is: + mseasons[s] = mdata.sel(time=seasons[s]).mean(dim='time') + oseasons[s] = odata.sel(time=seasons[s]).mean(dim='time') + #End if + pf.plot_map_and_save(plot_name, case_nickname, data.ref_nickname, [syear_cases[case_idx],eyear_cases[case_idx]], [syear_baseline,eyear_baseline], @@ -265,6 +266,16 @@ def global_latlon_map_B(adfobj): plot_name = plot_loc / f"{var}_{pres}hpa_{s}_LatLon_Mean.{plot_type}" if doplot[plot_name] is None: continue + + if weight_season: + mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) + oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) + else: + #Just average months as-is: + mseasons[s] = mdata.sel(time=seasons[s]).mean(dim='time') + oseasons[s] = odata.sel(time=seasons[s]).mean(dim='time') + #End if + pf.plot_map_and_save(plot_name, case_nickname, data.ref_nickname, [syear_cases[case_idx],eyear_cases[case_idx]], [syear_baseline,eyear_baseline], From cf2128c32fc554ed12b723dec50befd36c033ff4 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 13:35:32 -0700 Subject: [PATCH 30/46] seasonal averaging moved to inside season loop --- scripts/plotting/global_latlon_map_B.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/scripts/plotting/global_latlon_map_B.py b/scripts/plotting/global_latlon_map_B.py index 2b48eec09..fb3c4a4be 100644 --- a/scripts/plotting/global_latlon_map_B.py +++ b/scripts/plotting/global_latlon_map_B.py @@ -219,9 +219,6 @@ def global_latlon_map_B(adfobj): oseasons = {} dseasons = {} # hold the differences - # difference: each entry should be (lat, lon) - dseasons[s] = mseasons[s] - oseasons[s] - if not pres_levs: #Loop over season dictionary: @@ -239,6 +236,9 @@ def global_latlon_map_B(adfobj): oseasons[s] = odata.sel(time=seasons[s]).mean(dim='time') #End if + # difference: each entry should be (lat, lon) + dseasons[s] = mseasons[s] - oseasons[s] + pf.plot_map_and_save(plot_name, case_nickname, data.ref_nickname, [syear_cases[case_idx],eyear_cases[case_idx]], [syear_baseline,eyear_baseline], @@ -276,6 +276,9 @@ def global_latlon_map_B(adfobj): oseasons[s] = odata.sel(time=seasons[s]).mean(dim='time') #End if + # difference: each entry should be (lat, lon) + dseasons[s] = mseasons[s] - oseasons[s] + pf.plot_map_and_save(plot_name, case_nickname, data.ref_nickname, [syear_cases[case_idx],eyear_cases[case_idx]], [syear_baseline,eyear_baseline], From 3be03dbfeeba1a283f1291bd02ff5ad7395ffdf1 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 1 Dec 2023 13:37:49 -0700 Subject: [PATCH 31/46] correct zonal mean error message --- scripts/plotting/zonal_mean_B.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/plotting/zonal_mean_B.py b/scripts/plotting/zonal_mean_B.py index 1374dadd8..8580741fc 100644 --- a/scripts/plotting/zonal_mean_B.py +++ b/scripts/plotting/zonal_mean_B.py @@ -236,7 +236,7 @@ def zonal_mean_B(adfobj): print(f"Error: expecting lev for both case: {has_lev} and ref: {has_lev_ref}") continue if len(mdata['lev']) != len(odata['lev']): - print(f"Error: zonal mean contour expects `lev` dim to have same size, got {len(mdata['lev'])} and {odata['lev']}") + print(f"Error: zonal mean contour expects `lev` dim to have same size, got {len(mdata['lev'])} and {len(odata['lev'])}") continue plot_name_log = plot_loc / f"{var}_{s}_Zonal_logp_Mean.{plot_type}" if plot_name_log not in logp_zonal_skip: From 6b9f67138b5f92e78f0e00670d920aa5318178a3 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 19 Apr 2024 08:37:22 -0600 Subject: [PATCH 32/46] updated adf_dataset.py --- lib/adf_dataset.py | 49 ++++++++++++++++++++++++++++++++++------------ 1 file changed, 37 insertions(+), 12 deletions(-) diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py index 6f06bfd32..1043799c5 100644 --- a/lib/adf_dataset.py +++ b/lib/adf_dataset.py @@ -90,15 +90,14 @@ def set_reference(self): def get_reference_climo_file(self, var): """Return a list of files to be used as reference (aka baseline) for variable var.""" if self.reference_is_obs: - return [self.ref_var_loc[var]] - else: - self.ref_loc = self.adf.get_baseline_info("cam_climo_loc") - # NOTE: originally had this looking for *_baseline.nc - fils = sorted(Path(self.ref_loc).glob(f"{self.ref_case_label}_{var}_climo.nc")) - if fils: - return fils - else: - return None + fils = self.ref_var_loc.get(var, None) + return [fils] if fils is not None else None + self.ref_loc = self.adf.get_baseline_info("cam_climo_loc") + # NOTE: originally had this looking for *_baseline.nc + fils = sorted(Path(self.ref_loc).glob(f"{self.ref_case_label}_{var}_climo.nc")) + if fils: + return fils + return None def load_reference_dataset(self, var): fils = self.get_reference_climo_file(var) @@ -121,8 +120,29 @@ def load_reference_da(self, variablename): da.attrs['units'] = vres.get("new_unit", da.attrs.get('units', 'none')) return da - def get_climo_file(self, variablename): - pass + + def load_climo_da(self, case, variablename): + """Return DataArray from climo file""" + fils = self.get_climo_file(case, variablename) + return self.load_da(fils, variablename) + + + def load_climo_file(self, case, variablename): + """Return Dataset for climo of variablename""" + fils = self.get_climo_file(case, variablename) + if not fils: + print(f"ERROR: Did not find climo file for variable: {variablename}. Will try to skip.") + return None + return self.load_dataset(fils) + + + def get_climo_file(self, case, variablename): + """Retrieve the climo file path(s) for variablename for a specific case.""" + a = self.adf.get_cam_info("cam_climo_loc", required=True) # list of paths (could be multiple cases) + caseindex = (self.case_names).index(case) # the entry for specified case + # print(f"Checking if case name is in the climo loc entry: {case in a[caseindex]}") + model_cl_loc = Path(a[caseindex]) + return sorted(model_cl_loc.glob(f"{case}_{variablename}_climo.nc")) def get_timeseries_file(self, case, field): ts_locs = self.adf.get_cam_info("cam_ts_loc", required=True) @@ -171,13 +191,18 @@ def load_dataset(self, fils): ds = xr.open_mfdataset(fils, combine='by_coords') else: sfil = str(fils[0]) + assert Path(sfil).is_file(), f"Needs to be a file: {sfil}" ds = xr.open_dataset(sfil) if ds is None: warnings.warn(f"invalid data on load_dataset") return ds def load_da(self, fils, variablename): - da = (self.load_dataset(fils)[variablename]).squeeze() + ds = self.load_dataset(fils) + if ds is None: + print(f"ERROR: Load failed for {variablename}") + return None + da = (ds[variablename]).squeeze() if variablename in self.res: vres = self.res[variablename] da = da * vres.get("scale_factor",1) + vres.get("add_offset", 0) From 5d46906dd8e33b47b05a87b0fc04aa63093f5006 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Wed, 12 Jun 2024 14:29:30 -0600 Subject: [PATCH 33/46] starting to address PR comments. First round done on adf_dataset.py --- lib/adf_dataset.py | 48 ++++++++++++------------- scripts/plotting/global_latlon_map_B.py | 4 +-- scripts/plotting/zonal_mean_B.py | 4 +-- 3 files changed, 28 insertions(+), 28 deletions(-) diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py index 1043799c5..9ef774d62 100644 --- a/lib/adf_dataset.py +++ b/lib/adf_dataset.py @@ -50,8 +50,7 @@ def __init__(self, adfobj): self.model_rgrid_loc = adfobj.get_basic_info("cam_regrid_loc", required=True) # variables (and info for unit transform) - self.var_list = adfobj.diag_var_list - self.res = adfobj.variable_defaults + # use self.adf.diag_var_list and self.adf.self.adf.variable_defaults # case names and nicknames self.case_names = adfobj.get_cam_info("cam_case_name", required=True) @@ -60,27 +59,26 @@ def __init__(self, adfobj): self.ref_nickname = self.base_nickname # define reference data - self.reference_is_obs = adfobj.get_basic_info("compare_obs") self.set_reference() # specify "ref_labels" -> called "data_list" in zonal_mean (name of data source) def set_reference(self): """Set attributes for reference (aka baseline) data location, names, and variables.""" - if self.reference_is_obs: + if self.adf.compare_obs: self.ref_var_loc = {v: self.adf.var_obs_dict[v]['obs_file'] for v in self.adf.var_obs_dict} self.ref_labels = {v: self.adf.var_obs_dict[v]['obs_name'] for v in self.adf.var_obs_dict} self.ref_var_nam = {v: self.adf.var_obs_dict[v]['obs_var'] for v in self.adf.var_obs_dict} if not self.adf.var_obs_dict: - print("\t WARNING: reference is observations, but no observations found to plot against.") + warnings.warn("\t WARNING: reference is observations, but no observations found to plot against.") else: self.ref_var_loc = {} self.ref_var_nam = {} self.ref_labels = {} # when using a reference simulation, allow a "special" attribute with the case name: self.ref_case_label = self.adf.get_baseline_info("cam_case_name", required=True) - for v in self.var_list: + for v in self.adf.diag_var_list: f = self.get_reference_climo_file(v) if f is None: - print(f"\t WARNING: ADFData found no reference climo file for {v}") + warnings.warn(f"\t WARNING: ADFData found no reference climo file for {v}") continue else: self.ref_var_loc[v] = f @@ -89,12 +87,12 @@ def set_reference(self): def get_reference_climo_file(self, var): """Return a list of files to be used as reference (aka baseline) for variable var.""" - if self.reference_is_obs: + if self.adf.compare_obs: fils = self.ref_var_loc.get(var, None) return [fils] if fils is not None else None - self.ref_loc = self.adf.get_baseline_info("cam_climo_loc") + ref_loc = self.adf.get_baseline_info("cam_climo_loc") # NOTE: originally had this looking for *_baseline.nc - fils = sorted(Path(self.ref_loc).glob(f"{self.ref_case_label}_{var}_climo.nc")) + fils = sorted(Path(ref_loc).glob(f"{self.ref_case_label}_{var}_climo.nc")) if fils: return fils return None @@ -102,15 +100,15 @@ def get_reference_climo_file(self, var): def load_reference_dataset(self, var): fils = self.get_reference_climo_file(var) if not fils: - print(f"ERROR: Did not find any reference files for variable: {var}. Will try to skip.") + warnings.warn(f"ERROR: Did not find any reference files for variable: {var}. Will try to skip.") return None return self.load_dataset(fils) def load_reference_da(self, variablename): da = self.load_reference_dataset(variablename)[self.ref_var_nam[variablename]] - if variablename in self.res: - vres = self.res[variablename] - if self.reference_is_obs: + if variablename in self.adf.variable_defaults: + vres = self.adf.variable_defaults[variablename] + if self.adf.compare_obs: scale_factor = vres.get("obs_scale_factor",1) add_offset = vres.get("obs_add_offset", 0) else: @@ -145,14 +143,15 @@ def get_climo_file(self, case, variablename): return sorted(model_cl_loc.glob(f"{case}_{variablename}_climo.nc")) def get_timeseries_file(self, case, field): - ts_locs = self.adf.get_cam_info("cam_ts_loc", required=True) - ts_loc = Path(ts_locs[case]) + ts_locs = self.adf.get_cam_info("cam_ts_loc", required=True) # list of paths (could be multiple cases) + caseindex = (self.case_names).index(case) + ts_loc = Path(ts_locs[caseindex]) ts_filenames = f'{case}.*.{field}.*nc' ts_files = sorted(ts_loc.glob(ts_filenames)) return ts_files def get_ref_timeseries_file(self, field): - if self.reference_is_obs: + if self.adf.compare_obs: return None else: ts_loc = Path(self.adf.get_baseline_info("cam_ts_loc", required=True)) @@ -169,19 +168,17 @@ def get_regrid_file(self, case, field): def load_regrid_dataset(self, case, field): fils = self.get_regrid_file(case, field) if not fils: - print(f"ERROR: Did not find regrid file(s) for case: {case}, variable: {field}") + warnings.warn(f"ERROR: Did not find regrid file(s) for case: {case}, variable: {field}") return None return self.load_dataset(fils) def load_regrid_da(self, case, field): fils = self.get_regrid_file(case, field) if not fils: - print(f"ERROR: Did not find regrid file(s) for case: {case}, variable: {field}") + warnings.warn(f"ERROR: Did not find regrid file(s) for case: {case}, variable: {field}") return None return self.load_da(fils, field) - def get_file_list(): - pass def load_dataset(self, fils): if (len(fils) == 0): @@ -191,20 +188,23 @@ def load_dataset(self, fils): ds = xr.open_mfdataset(fils, combine='by_coords') else: sfil = str(fils[0]) - assert Path(sfil).is_file(), f"Needs to be a file: {sfil}" + if not Path(sfil).is_file(): + warnings.warn(f"Expecting to find file: {sfil}") + return None ds = xr.open_dataset(sfil) if ds is None: warnings.warn(f"invalid data on load_dataset") return ds + def load_da(self, fils, variablename): ds = self.load_dataset(fils) if ds is None: print(f"ERROR: Load failed for {variablename}") return None da = (ds[variablename]).squeeze() - if variablename in self.res: - vres = self.res[variablename] + if variablename in self.adf.variable_defaults: + vres = self.adf.variable_defaults[variablename] da = da * vres.get("scale_factor",1) + vres.get("add_offset", 0) da.attrs['units'] = vres.get("new_unit", da.attrs.get('units', 'none')) return da diff --git a/scripts/plotting/global_latlon_map_B.py b/scripts/plotting/global_latlon_map_B.py index fb3c4a4be..9b93ab76b 100644 --- a/scripts/plotting/global_latlon_map_B.py +++ b/scripts/plotting/global_latlon_map_B.py @@ -243,7 +243,7 @@ def global_latlon_map_B(adfobj): [syear_cases[case_idx],eyear_cases[case_idx]], [syear_baseline,eyear_baseline], mseasons[s], oseasons[s], dseasons[s], - obs=data.reference_is_obs, **vres) + obs=data.adf.compare_obs, **vres) #Add plot to website (if enabled): adfobj.add_website_data(plot_name, var, case_name, category=web_category, @@ -283,7 +283,7 @@ def global_latlon_map_B(adfobj): [syear_cases[case_idx],eyear_cases[case_idx]], [syear_baseline,eyear_baseline], mseasons[s].sel(lev=pres), oseasons[s].sel(lev=pres), dseasons[s].sel(lev=pres), - obs=data.reference_is_obs, **vres) + obs=data.adf.compare_obs, **vres) #Add plot to website (if enabled): adfobj.add_website_data(plot_name, f"{var}_{pres}hpa", case_name, category=web_category, diff --git a/scripts/plotting/zonal_mean_B.py b/scripts/plotting/zonal_mean_B.py index 8580741fc..ffed220e5 100644 --- a/scripts/plotting/zonal_mean_B.py +++ b/scripts/plotting/zonal_mean_B.py @@ -224,7 +224,7 @@ def zonal_mean_B(adfobj): pf.plot_zonal_mean_and_save(plot_name, case_nickname, data.ref_nickname, [syear_cases[case_idx],eyear_cases[case_idx]], [syear_baseline,eyear_baseline], - mseasons[s], oseasons[s], has_lev, log_p=False, obs=data.reference_is_obs, **vres) + mseasons[s], oseasons[s], has_lev, log_p=False, obs=data.adf.compare_obs, **vres) #Add plot to website (if enabled): adfobj.add_website_data(plot_name, var, case_name, season=s, plot_type="Zonal") @@ -247,7 +247,7 @@ def zonal_mean_B(adfobj): pf.plot_zonal_mean_and_save(plot_name_log, case_nickname, data.ref_nickname, [syear_cases[case_idx],eyear_cases[case_idx]], [syear_baseline,eyear_baseline], - mseasons[s], oseasons[s], has_lev, log_p=True, obs=data.reference_is_obs, **vres) + mseasons[s], oseasons[s], has_lev, log_p=True, obs=data.adf.compare_obs, **vres) #Add plot to website (if enabled): adfobj.add_website_data(plot_name_log, f"{var}_logp", case_name, season=s, plot_type="Zonal", category="Log-P") From ce50839f135b84b1876d396ed8c77ddfe66808a3 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Wed, 12 Jun 2024 14:44:37 -0600 Subject: [PATCH 34/46] merged my conflicted adf_dataset.py --- lib/adf_dataset.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py index 9ef774d62..91ec00f52 100644 --- a/lib/adf_dataset.py +++ b/lib/adf_dataset.py @@ -129,7 +129,7 @@ def load_climo_file(self, case, variablename): """Return Dataset for climo of variablename""" fils = self.get_climo_file(case, variablename) if not fils: - print(f"ERROR: Did not find climo file for variable: {variablename}. Will try to skip.") + warnings.warning(f"ERROR: Did not find climo file for variable: {variablename}. Will try to skip.") return None return self.load_dataset(fils) @@ -138,7 +138,6 @@ def get_climo_file(self, case, variablename): """Retrieve the climo file path(s) for variablename for a specific case.""" a = self.adf.get_cam_info("cam_climo_loc", required=True) # list of paths (could be multiple cases) caseindex = (self.case_names).index(case) # the entry for specified case - # print(f"Checking if case name is in the climo loc entry: {case in a[caseindex]}") model_cl_loc = Path(a[caseindex]) return sorted(model_cl_loc.glob(f"{case}_{variablename}_climo.nc")) @@ -200,7 +199,7 @@ def load_dataset(self, fils): def load_da(self, fils, variablename): ds = self.load_dataset(fils) if ds is None: - print(f"ERROR: Load failed for {variablename}") + warnings.warn(f"ERROR: Load failed for {variablename}") return None da = (ds[variablename]).squeeze() if variablename in self.adf.variable_defaults: From 142d4e5c779da4bc272fdbef94b1a9145b2a7392 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Wed, 12 Jun 2024 15:19:55 -0600 Subject: [PATCH 35/46] addressing PR comments. Instantiate AdfData from AdfDiag --- lib/adf_diag.py | 7 +++++-- lib/plotting_functions.py | 2 -- scripts/averaging/create_climo_files.py | 1 - scripts/plotting/global_latlon_map_B.py | 26 ++++++++++++------------- scripts/plotting/zonal_mean_B.py | 24 +++++++++++------------ 5 files changed, 30 insertions(+), 30 deletions(-) diff --git a/lib/adf_diag.py b/lib/adf_diag.py index 9c57c1f36..38b6cb0dd 100644 --- a/lib/adf_diag.py +++ b/lib/adf_diag.py @@ -92,9 +92,9 @@ # +++++++++++++++++++++++++++++ -# Finally, import needed ADF module: +# Finally, import needed ADF modules: from adf_web import AdfWeb - +from adf_dataset import AdfData ################# # Helper functions @@ -182,6 +182,9 @@ def __init__(self, config_file, debug=False): # Add plotting script names: self.__plotting_scripts = self.read_config_var("plotting_scripts") + # Provide convenience functions for data handling: + self.data = AdfData(self) + # Create property needed to return "plotting_scripts" variable to user: @property def plotting_scripts(self): diff --git a/lib/plotting_functions.py b/lib/plotting_functions.py index 71c3cea96..05a7d32b3 100644 --- a/lib/plotting_functions.py +++ b/lib/plotting_functions.py @@ -2317,9 +2317,7 @@ def square_contour_difference(fld1, fld2, **kwargs): mnorm = mpl.colors.Normalize(mn, mx) coord1, coord2 = fld1.coords # ASSUMES xarray WITH coords AND 2-dimensions - # print(f"{coord1}, {coord2}") xx, yy = np.meshgrid(fld1[coord2], fld1[coord1]) - # print(f"shape of meshgrid: {xx.shape}") img1 = ax1.contourf(xx, yy, fld1.transpose()) if (coord1 == 'month') and (fld1.shape[0] ==12): diff --git a/scripts/averaging/create_climo_files.py b/scripts/averaging/create_climo_files.py index 78c6d257b..1591fa93f 100644 --- a/scripts/averaging/create_climo_files.py +++ b/scripts/averaging/create_climo_files.py @@ -214,7 +214,6 @@ def process_variable(ts_files, syr, eyr, output_file): cam_ts_data.assign_coords(time=time) cam_ts_data = xr.decode_cf(cam_ts_data) #Extract data subset using provided year bounds: - #cam_ts_data = cam_ts_data.sel(time=slice(syr, eyr)) tslice = get_time_slice_by_year(cam_ts_data.time, int(syr), int(eyr)) cam_ts_data = cam_ts_data.isel(time=tslice) #Group time series values by month, and average those months together: diff --git a/scripts/plotting/global_latlon_map_B.py b/scripts/plotting/global_latlon_map_B.py index 9b93ab76b..8c68d828a 100644 --- a/scripts/plotting/global_latlon_map_B.py +++ b/scripts/plotting/global_latlon_map_B.py @@ -18,7 +18,7 @@ import warnings # use to warn user about missing files. import plotting_functions as pf -from adf_dataset import AdfData +# from adf_dataset import AdfData #Format warning messages: def my_formatwarning(msg, *args, **kwargs): @@ -83,7 +83,7 @@ def global_latlon_map_B(adfobj): # # Use ADF api to get all necessary information # - data = AdfData(adfobj) + # data = AdfData(adfobj) NO LONGER NEEDED var_list = adfobj.diag_var_list #Special ADF variable which contains the output paths for #all generated plots and tables for each case: @@ -128,7 +128,7 @@ def global_latlon_map_B(adfobj): # probably want to do this one variable at a time: for var in var_list: - if var not in data.ref_var_nam: + if var not in adfobj.data.ref_var_nam: dmsg = f"No reference data found for variable `{var}`, zonal mean plotting skipped." adfobj.debug_log(dmsg) continue @@ -156,7 +156,7 @@ def global_latlon_map_B(adfobj): vres['central_longitude'] = pf.get_central_longitude(adfobj) # load reference data (observational or baseline) - odata = data.load_reference_da(var) + odata = adfobj.data.load_reference_da(var) if odata is None: continue has_dims = pf.lat_lon_validate_dims(odata) # T iff dims are (lat,lon) -- can't plot unless we have both @@ -165,10 +165,10 @@ def global_latlon_map_B(adfobj): continue #Loop over model cases: - for case_idx, case_name in enumerate(data.case_names): + for case_idx, case_name in enumerate(adfobj.data.case_names): #Set case nickname: - case_nickname = data.test_nicknames[case_idx] + case_nickname = adfobj.data.test_nicknames[case_idx] #Set output plot location: plot_loc = Path(plot_locations[case_idx]) @@ -179,7 +179,7 @@ def global_latlon_map_B(adfobj): plot_loc.mkdir(parents=True) #Load re-gridded model files: - mdata = data.load_regrid_da(case_name, var) + mdata = adfobj.data.load_regrid_da(case_name, var) #Skip this variable/case if the regridded climo file doesn't exist: if mdata is None: @@ -239,11 +239,11 @@ def global_latlon_map_B(adfobj): # difference: each entry should be (lat, lon) dseasons[s] = mseasons[s] - oseasons[s] - pf.plot_map_and_save(plot_name, case_nickname, data.ref_nickname, + pf.plot_map_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], dseasons[s], - obs=data.adf.compare_obs, **vres) + obs=adfobj.compare_obs, **vres) #Add plot to website (if enabled): adfobj.add_website_data(plot_name, var, case_name, category=web_category, @@ -283,7 +283,7 @@ def global_latlon_map_B(adfobj): [syear_cases[case_idx],eyear_cases[case_idx]], [syear_baseline,eyear_baseline], mseasons[s].sel(lev=pres), oseasons[s].sel(lev=pres), dseasons[s].sel(lev=pres), - obs=data.adf.compare_obs, **vres) + obs=adfobj.compare_obs, **vres) #Add plot to website (if enabled): adfobj.add_website_data(plot_name, f"{var}_{pres}hpa", case_name, category=web_category, @@ -343,14 +343,14 @@ def plot_file_op(adfobj, plot_name, var, case_name, season, web_category, redo_p if plot_name.is_file(): if redo_plot: plot_name.unlink() - return 1 + return True else: #Add already-existing plot to website (if enabled): adfobj.add_website_data(plot_name, var, case_name, category=web_category, season=season, plot_type=plot_type) - return None # None tells caller that file exists and not to overwrite + return False # False tells caller that file exists and not to overwrite else: - return 1 + return True ############## #END OF SCRIPT diff --git a/scripts/plotting/zonal_mean_B.py b/scripts/plotting/zonal_mean_B.py index ffed220e5..998c91c80 100644 --- a/scripts/plotting/zonal_mean_B.py +++ b/scripts/plotting/zonal_mean_B.py @@ -2,7 +2,7 @@ import numpy as np import xarray as xr import plotting_functions as pf -from adf_dataset import AdfData +# from adf_dataset import AdfData import warnings # use to warn user about missing files. @@ -47,7 +47,7 @@ def zonal_mean_B(adfobj): # Start by instantiating the AdfData object # and Extract needed quantities from ADF object: - data = AdfData(adfobj) + # data = AdfData(adfobj) var_list = adfobj.diag_var_list @@ -93,7 +93,7 @@ def zonal_mean_B(adfobj): logp_zonal_skip = [] #Loop over model cases: - for case_idx, case_name in enumerate(data.case_names): + for case_idx, case_name in enumerate(adfobj.data.case_names): #Set output plot location: plot_loc = Path(plot_locations[case_idx]) @@ -146,7 +146,7 @@ def zonal_mean_B(adfobj): # #Loop over variables: for var in var_list: - if var not in data.ref_var_nam: + if var not in adfobj.data.ref_var_nam: dmsg = f"No obs found for variable `{var}`, zonal mean plotting skipped." adfobj.debug_log(dmsg) continue @@ -165,20 +165,20 @@ def zonal_mean_B(adfobj): #End if # load reference data (observational or baseline) - odata = data.load_reference_da(var) + odata = adfobj.data.load_reference_da(var) has_lat_ref, has_lev_ref = pf.zm_validate_dims(odata) #Loop over model cases: - for case_idx, case_name in enumerate(data.case_names): + for case_idx, case_name in enumerate(adfobj.data.case_names): #Set case nickname: - case_nickname = data.test_nicknames[case_idx] + case_nickname = adfobj.data.test_nicknames[case_idx] #Set output plot location: plot_loc = Path(plot_locations[case_idx]) # load re-gridded model files: - mdata = data.load_regrid_da(case_name, var) + mdata = adfobj.data.load_regrid_da(case_name, var) # determine whether it's 2D or 3D # 3D triggers search for surface pressure @@ -221,10 +221,10 @@ def zonal_mean_B(adfobj): # This could be re-visited for efficiency or improved code structure. #Create new plot: - pf.plot_zonal_mean_and_save(plot_name, case_nickname, data.ref_nickname, + pf.plot_zonal_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, log_p=False, obs=data.adf.compare_obs, **vres) + mseasons[s], oseasons[s], has_lev, log_p=False, obs=adfobj.compare_obs, **vres) #Add plot to website (if enabled): adfobj.add_website_data(plot_name, var, case_name, season=s, plot_type="Zonal") @@ -244,10 +244,10 @@ def zonal_mean_B(adfobj): mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) - pf.plot_zonal_mean_and_save(plot_name_log, case_nickname, data.ref_nickname, + pf.plot_zonal_mean_and_save(plot_name_log, case_nickname, adfobj.data.ref_nickname, [syear_cases[case_idx],eyear_cases[case_idx]], [syear_baseline,eyear_baseline], - mseasons[s], oseasons[s], has_lev, log_p=True, obs=data.adf.compare_obs, **vres) + mseasons[s], oseasons[s], has_lev, log_p=True, obs=adfobj.compare_obs, **vres) #Add plot to website (if enabled): adfobj.add_website_data(plot_name_log, f"{var}_logp", case_name, season=s, plot_type="Zonal", category="Log-P") From cae606fe52851afef6c5c9c381935834ce0333bf Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Wed, 12 Jun 2024 17:04:18 -0600 Subject: [PATCH 36/46] testing AdfData. Fixes for history file naming. Fix to force loading data for derived time series. --- lib/adf_diag.py | 11 +++++++---- lib/adf_info.py | 24 +++++++++++++++++++++--- 2 files changed, 28 insertions(+), 7 deletions(-) diff --git a/lib/adf_diag.py b/lib/adf_diag.py index 38b6cb0dd..983c21bf0 100644 --- a/lib/adf_diag.py +++ b/lib/adf_diag.py @@ -541,11 +541,13 @@ def call_ncrcat(cmd): # Aerosol Calcs #-------------- #Always make sure PMID is made if aerosols are desired in config file + # Since there's no requirement for `aerosol_zonal_list` to be included, allow it to be absent: + azl = res.get("aerosol_zonal_list", []) if "PMID" not in diag_var_list: - if any(item in res["aerosol_zonal_list"] for item in diag_var_list): + if any(item in azl for item in diag_var_list): diag_var_list += ["PMID"] if "T" not in diag_var_list: - if any(item in res["aerosol_zonal_list"] for item in diag_var_list): + if any(item in azl for item in diag_var_list): diag_var_list += ["T"] #End aerosol calcs @@ -1059,7 +1061,7 @@ def derive_variables(self, res=None, vars_to_derive=None, ts_dir=None, overwrite print(ermsg) else: #Open a new dataset with all the constituent files/variables - ds = xr.open_mfdataset(constit_files) + ds = xr.open_mfdataset(constit_files).compute() # create new file name for derived variable derived_file = constit_files[0].replace(constit_list[0], var) @@ -1091,7 +1093,8 @@ def derive_variables(self, res=None, vars_to_derive=None, ts_dir=None, overwrite #These will be multiplied by rho (density of dry air) ds_pmid_done = False ds_t_done = False - if var in res["aerosol_zonal_list"]: + azl = res.get("aerosol_zonal_list", []) # User-defined defaults might not include aerosol zonal list + if var in azl: #Only calculate once for all aerosol vars if not ds_pmid_done: diff --git a/lib/adf_info.py b/lib/adf_info.py index 42fe47fc3..5081ddd93 100644 --- a/lib/adf_info.py +++ b/lib/adf_info.py @@ -673,20 +673,38 @@ def get_climo_yrs_from_ts(self, input_ts_loc, case_name): #Search for first variable in var_list to get a time series file to read #NOTE: it is assumed all the variables have the same dates! - ts_files = sorted(input_location.glob(f"{case_name}*.{var_list[0]}.*nc")) + + #Read hist_str (component.hist_num) from the yaml file, or set to default + hist_str = self.get_basic_info('hist_str') + #If hist_str is not present, then default to 'cam.h0': + if not hist_str: + hist_str = 'cam.h0' + #End if + + ts_files = sorted(input_location.glob(f"{case_name}.{hist_str}.{var_list[0]}.*nc")) #Read in file(s) if len(ts_files) == 1: cam_ts_data = xr.open_dataset(ts_files[0], decode_times=True) else: - cam_ts_data = xr.open_mfdataset(ts_files, decode_times=True, combine='by_coords') + try: + cam_ts_data = xr.open_mfdataset(ts_files, decode_times=True, combine='by_coords') + except: + print(" ----------- ERROR ------------") + print(ts_files) #Average time dimension over time bounds, if bounds exist: if 'time_bnds' in cam_ts_data: + timeBoundsName = 'time_bnds' + elif 'time_bounds' in cam_ts_data: + timeBoundsName = 'time_bounds' + else: + timeBoundsName = None + if timeBoundsName: time = cam_ts_data['time'] #NOTE: force `load` here b/c if dask & time is cftime, #throws a NotImplementedError: - time = xr.DataArray(cam_ts_data['time_bnds'].load().mean(dim='nbnd').values, dims=time.dims, attrs=time.attrs) + time = xr.DataArray(cam_ts_data[timeBoundsName].load().mean(dim='nbnd').values, dims=time.dims, attrs=time.attrs) cam_ts_data['time'] = time cam_ts_data.assign_coords(time=time) cam_ts_data = xr.decode_cf(cam_ts_data) From 0183799260f627769300ce7212b3f09850f85861 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Wed, 12 Jun 2024 17:23:54 -0600 Subject: [PATCH 37/46] removed commented lines --- scripts/plotting/global_latlon_map_B.py | 1 - scripts/plotting/zonal_mean_B.py | 1 - 2 files changed, 2 deletions(-) diff --git a/scripts/plotting/global_latlon_map_B.py b/scripts/plotting/global_latlon_map_B.py index 8c68d828a..c8444c6de 100644 --- a/scripts/plotting/global_latlon_map_B.py +++ b/scripts/plotting/global_latlon_map_B.py @@ -18,7 +18,6 @@ import warnings # use to warn user about missing files. import plotting_functions as pf -# from adf_dataset import AdfData #Format warning messages: def my_formatwarning(msg, *args, **kwargs): diff --git a/scripts/plotting/zonal_mean_B.py b/scripts/plotting/zonal_mean_B.py index 998c91c80..1b78841c3 100644 --- a/scripts/plotting/zonal_mean_B.py +++ b/scripts/plotting/zonal_mean_B.py @@ -2,7 +2,6 @@ import numpy as np import xarray as xr import plotting_functions as pf -# from adf_dataset import AdfData import warnings # use to warn user about missing files. From 26a059c330149d9641f19115815deafddf89dc8d Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Tue, 18 Jun 2024 14:14:37 -0600 Subject: [PATCH 38/46] addressing Jesse's comments on PR --- lib/adf_dataset.py | 2 +- lib/adf_info.py | 6 +----- scripts/plotting/global_latlon_map_B.py | 1 - scripts/plotting/zonal_mean_B.py | 4 ---- 4 files changed, 2 insertions(+), 11 deletions(-) diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py index 91ec00f52..5fd9521bc 100644 --- a/lib/adf_dataset.py +++ b/lib/adf_dataset.py @@ -129,7 +129,7 @@ def load_climo_file(self, case, variablename): """Return Dataset for climo of variablename""" fils = self.get_climo_file(case, variablename) if not fils: - warnings.warning(f"ERROR: Did not find climo file for variable: {variablename}. Will try to skip.") + warnings.warn(f"ERROR: Did not find climo file for variable: {variablename}. Will try to skip.") return None return self.load_dataset(fils) diff --git a/lib/adf_info.py b/lib/adf_info.py index 5081ddd93..7fed436b7 100644 --- a/lib/adf_info.py +++ b/lib/adf_info.py @@ -687,11 +687,7 @@ def get_climo_yrs_from_ts(self, input_ts_loc, case_name): if len(ts_files) == 1: cam_ts_data = xr.open_dataset(ts_files[0], decode_times=True) else: - try: - cam_ts_data = xr.open_mfdataset(ts_files, decode_times=True, combine='by_coords') - except: - print(" ----------- ERROR ------------") - print(ts_files) + cam_ts_data = xr.open_mfdataset(ts_files, decode_times=True, combine='by_coords') #Average time dimension over time bounds, if bounds exist: if 'time_bnds' in cam_ts_data: diff --git a/scripts/plotting/global_latlon_map_B.py b/scripts/plotting/global_latlon_map_B.py index c8444c6de..2e13d62d5 100644 --- a/scripts/plotting/global_latlon_map_B.py +++ b/scripts/plotting/global_latlon_map_B.py @@ -82,7 +82,6 @@ def global_latlon_map_B(adfobj): # # Use ADF api to get all necessary information # - # data = AdfData(adfobj) NO LONGER NEEDED var_list = adfobj.diag_var_list #Special ADF variable which contains the output paths for #all generated plots and tables for each case: diff --git a/scripts/plotting/zonal_mean_B.py b/scripts/plotting/zonal_mean_B.py index 1b78841c3..019065d9f 100644 --- a/scripts/plotting/zonal_mean_B.py +++ b/scripts/plotting/zonal_mean_B.py @@ -44,10 +44,6 @@ def zonal_mean_B(adfobj): print("\n Generating zonal mean plots...") - # Start by instantiating the AdfData object - # and Extract needed quantities from ADF object: - # data = AdfData(adfobj) - var_list = adfobj.diag_var_list #Special ADF variable which contains the output paths for From 7d6cdc9c1d236f3abc44326e2433928194c73e54 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 21 Jun 2024 13:28:27 -0600 Subject: [PATCH 39/46] Replace and with versions using new class --- scripts/plotting/global_latlon_map.py | 516 +++++++++--------------- scripts/plotting/global_latlon_map_B.py | 354 ---------------- scripts/plotting/zonal_mean.py | 281 +++++-------- scripts/plotting/zonal_mean_B.py | 260 ------------ 4 files changed, 288 insertions(+), 1123 deletions(-) delete mode 100644 scripts/plotting/global_latlon_map_B.py delete mode 100644 scripts/plotting/zonal_mean_B.py diff --git a/scripts/plotting/global_latlon_map.py b/scripts/plotting/global_latlon_map.py index e04e1b91e..2e13d62d5 100644 --- a/scripts/plotting/global_latlon_map.py +++ b/scripts/plotting/global_latlon_map.py @@ -8,6 +8,8 @@ my_formatwarning(msg, *args, **kwargs) format warning messages (private method) +plot_file_op + Check on status of output plot file. """ #Import standard modules: from pathlib import Path @@ -15,13 +17,15 @@ import xarray as xr import warnings # use to warn user about missing files. +import plotting_functions as pf + #Format warning messages: def my_formatwarning(msg, *args, **kwargs): """Issue `msg` as warning.""" return str(msg) + '\n' warnings.formatwarning = my_formatwarning -def global_latlon_map(adfobj): +def global_latlon_map_B(adfobj): """ This script/function is designed to generate global 2-D lat/lon maps of model fields with continental overlays. @@ -37,32 +41,22 @@ def global_latlon_map(adfobj): Notes ----- - This function imports `pandas` and `plotting_functions` It uses the AdfDiag object's methods to get necessary information. - Specificially: + This version passes the AdfDiag object to instatiate an AdfData object. + Explicitly accesses: adfobj.diag_var_list List of variables - adfobj.get_basic_info - Regrid data path, checks `compare_obs`, checks `redo_plot`, checks `plot_press_levels` adfobj.plot_location output plot path - adfobj.get_cam_info - Get `cam_case_name` and `case_nickname` adfobj.climo_yrs start and end climo years of the case(s), `syears` & `eyears` start and end climo years of the reference, `syear_baseline` & `eyear_baseline` - adfobj.var_obs_dict - reference data (conditional) - adfobj.get_baseline_info - get reference case, `cam_case_name` adfobj.variable_defaults dict of variable-specific plot preferences adfobj.read_config_var dict of basic info, `diag_basic_info` Then use to check `plot_type` - adfobj.compare_obs - Used to set data path adfobj.debug_log Issues debug message adfobj.add_website_data @@ -82,20 +76,6 @@ def global_latlon_map(adfobj): Checks on pressure level dimension """ - #Import necessary modules: - #------------------------ - import pandas as pd - - #CAM diagnostic plotting functions: - import plotting_functions as pf - #------------------------- - - # Steps: - # - load regridded climo files for model and obs - # - calculate all-time and seasonal fields (from individual months) - # - Take difference, calculate statistics - # - make plot - #Notify user that script has started: print("\n Generating lat/lon maps...") @@ -103,50 +83,18 @@ def global_latlon_map(adfobj): # Use ADF api to get all necessary information # var_list = adfobj.diag_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 for each case: 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("No observations found to plot against, so no lat/lon 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 YAML. @@ -162,19 +110,11 @@ def global_latlon_map(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) - #----------------------- - #Determine if user wants to plot 3-D variables on #pressure levels: pres_levs = adfobj.get_basic_info("plot_press_levels") - #For now, let's always do seasonal weighting: - weight_season = True + weight_season = True #always do seasonal weighting #Set seasonal ranges: seasons = {"ANN": np.arange(1,13,1), @@ -183,29 +123,13 @@ def global_latlon_map(adfobj): "MAM": [3, 4, 5], "SON": [9, 10, 11] } - + # probably want to do this one variable at a time: for var in var_list: - - 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}`, lat/lon map plotting skipped." - adfobj.debug_log(dmsg) - continue - else: - #Set "data_var" for consistent use below: - data_var = var - #End if + if var not in adfobj.data.ref_var_nam: + dmsg = f"No reference data found for variable `{var}`, zonal mean plotting skipped." + adfobj.debug_log(dmsg) + continue #Notify user of variable being plotted: print("\t - lat/lon maps for {}".format(var)) @@ -229,258 +153,202 @@ def global_latlon_map(adfobj): # otherwise defaults to 180 vres['central_longitude'] = pf.get_central_longitude(adfobj) - #loop over different data sets to plot model against: - for data_src in data_list: + # load reference data (observational or baseline) + odata = adfobj.data.load_reference_da(var) + if odata is None: + continue + has_dims = pf.lat_lon_validate_dims(odata) # T iff dims are (lat,lon) -- can't plot unless we have both + if not has_dims: + print(f"\t = skipping global map for {var} as REFERENCE does not have both lat and lon") + continue - # load data (observational) commparison 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")) + #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]) - oclim_ds = pf.load_dataset(oclim_fils) - if oclim_ds is None: - print("WARNING: Did not find any oclim_fils. Will try to skip.") - print(f"INFO: Data Location, dclimo_loc is {dclimo_loc}") - print(f"INFO: The glob is: {data_src}_{var}_*.nc") + #Check if plot output directory exists, and if not, then create it: + if not plot_loc.is_dir(): + print(" {} not found, making new directory".format(plot_loc)) + plot_loc.mkdir(parents=True) + + #Load re-gridded model files: + mdata = adfobj.data.load_regrid_da(case_name, var) + + #Skip this variable/case if the regridded climo file doesn't exist: + if mdata is None: continue - #End if - #Loop over model cases: - for case_idx, case_name in enumerate(case_names): + #Determine dimensions of variable: + has_dims_cam = pf.lat_lon_validate_dims(mdata) # T iff dims are (lat,lon) -- can't plot unless we have both + _, has_lev = pf.zm_validate_dims(mdata) # has_lev T if lev in mdata + if not has_dims_cam: + print(f"\t = skipping global map for {var} for case {case_name} as it does not have both lat and lon") + continue + else: # i.e., has lat&lon + if pres_levs and (not has_lev): + print(f"\t - skipping global map for {var} as it has more than lat/lon dims, but no pressure levels were provided") + continue - #Set case nickname: - case_nickname = test_nicknames[case_idx] + # Check output file. If file does not exist, proceed. + # If file exists: + # if redo_plot is true: delete it now and make plot + # if redo_plot is false: add to website and move on + doplot = {} + if not pres_levs: + for s in seasons: + plot_name = plot_loc / f"{var}_{s}_LatLon_Mean.{plot_type}" + doplot[plot_name] = plot_file_op(adfobj, plot_name, var, case_name, s, web_category, redo_plot, "LatLon") + else: + for pres in pres_levs: + for s in seasons: + plot_name = plot_loc / f"{var}_{pres}hpa_{s}_LatLon_Mean.{plot_type}" + doplot[plot_name] = plot_file_op(adfobj, plot_name, f"{var}_{pres}hpa", case_name, s, web_category, redo_plot, "LatLon") + if all(value is None for value in doplot.values()): + print(f"All plots exist for {var}. Redo is {redo_plot}. Existing plots added to website data. Continue.") + continue - #Set output plot location: - plot_loc = Path(plot_locations[case_idx]) + #Create new dictionaries: + mseasons = {} + oseasons = {} + dseasons = {} # hold the differences - #Check if plot output directory exists, and if not, then create it: - if not plot_loc.is_dir(): - print(" {} not found, making new directory".format(plot_loc)) - plot_loc.mkdir(parents=True) + if not pres_levs: - #Load re-gridded model files: - mclim_fils = sorted(mclimo_rg_loc.glob(f"{data_src}_{case_name}_{var}_*.nc")) - mclim_ds = pf.load_dataset(mclim_fils) + #Loop over season dictionary: + for s in seasons: + plot_name = plot_loc / f"{var}_{s}_LatLon_Mean.{plot_type}" + if doplot[plot_name] is None: + continue + + if weight_season: + mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) + oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) + else: + #Just average months as-is: + mseasons[s] = mdata.sel(time=seasons[s]).mean(dim='time') + oseasons[s] = odata.sel(time=seasons[s]).mean(dim='time') + #End if + + # difference: each entry should be (lat, lon) + dseasons[s] = mseasons[s] - oseasons[s] + + pf.plot_map_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], dseasons[s], + obs=adfobj.compare_obs, **vres) + + #Add plot to website (if enabled): + adfobj.add_website_data(plot_name, var, case_name, category=web_category, + season=s, plot_type="LatLon") + + else: # => pres_levs has values, & we already checked that lev is in mdata (has_lev) + + for pres in pres_levs: + + #Check that the user-requested pressure level + #exists in the model data, which should already + #have been interpolated to the standard reference + #pressure levels: + if not (pres in mdata['lev']): + print(f"plot_press_levels value '{pres}' not present in {var}, so skipping.") + continue + + #Loop over seasons: + for s in seasons: + plot_name = plot_loc / f"{var}_{pres}hpa_{s}_LatLon_Mean.{plot_type}" + if doplot[plot_name] is None: + continue - #Skip this variable/case if the regridded climo file doesn't exist: - if mclim_ds is None: - print("WARNING: Did not find any regridded climo files. Will try to skip.") - print(f"INFO: Data Location, mclimo_rg_loc, is {mclimo_rg_loc}") - print(f"INFO: The glob is: {data_src}_{case_name}_{var}_*.nc") - 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 dimensions of variable: - has_dims = pf.lat_lon_validate_dims(odata) - if has_dims: - #If observations/baseline CAM have the correct - #dimensions, does the input CAM run have correct - #dimensions as well? - has_dims_cam = pf.lat_lon_validate_dims(mdata) - - #If both fields have the required dimensions, then - #proceed with plotting: - if has_dims_cam: - - # - # Seasonal Averages - # Note: xarray can do seasonal averaging, - # but depends on having time accessor, - # which these prototype climo files do not have. - # - - #Create new dictionaries: - mseasons = {} - oseasons = {} - dseasons = {} # hold the differences - - #Loop over season dictionary: - for s in seasons: - # time to make plot; here we'd probably loop over whatever plots we want for this variable - # I'll just call this one "LatLon_Mean" ... would this work as a pattern [operation]_[AxesDescription] ? - plot_name = plot_loc / f"{var}_{s}_LatLon_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, category=web_category, - season=s, plot_type="LatLon") - - #Continue to next iteration: - continue - elif (redo_plot) and plot_name.is_file(): - plot_name.unlink() - - - if weight_season: - mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) - oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) - else: - #Just average months as-is: - mseasons[s] = mdata.sel(time=seasons[s]).mean(dim='time') - oseasons[s] = odata.sel(time=seasons[s]).mean(dim='time') - #End if - - # difference: each entry should be (lat, lon) - dseasons[s] = mseasons[s] - oseasons[s] - - - #Create new plot: - # NOTE: send vres as kwarg dictionary. --> ONLY vres, not the full res - # This relies on `plot_map_and_save` knowing how to deal with the options - # currently knows how to handle: - # colormap, contour_levels, diff_colormap, diff_contour_levels, tiString, tiFontSize, mpl - # *Any other entries will be ignored. - # NOTE: If we were doing all the plotting here, we could use whatever we want from the provided YAML file. - - pf.plot_map_and_save(plot_name, case_nickname, base_nickname, - [syear_cases[case_idx],eyear_cases[case_idx]], - [syear_baseline,eyear_baseline], - mseasons[s], oseasons[s], dseasons[s], - obs, **vres) - - #Add plot to website (if enabled): - adfobj.add_website_data(plot_name, var, case_name, category=web_category, - season=s, plot_type="LatLon") - - else: #mdata dimensions check - print(f"\t - skipping lat/lon map for {var} as it doesn't have only lat/lon dims.") - #End if (dimensions check) - - elif pres_levs: #Is the user wanting to interpolate to a specific pressure level? - - #Check that case inputs have the correct dimensions (including "lev"): - _, has_lev = pf.zm_validate_dims(mdata) - - if has_lev: - - #Calculate monthly weights (if applicable): if weight_season: - #Add date-stamp to time dimension: - #Note: For now using made-up dates, but in the future - #it might be good to extract this info from the files - #themselves. - timefix = pd.date_range(start='1/1/1980', end='12/1/1980', freq='MS') - mdata['time']=timefix - odata['time']=timefix - - #Calculate monthly weights based on number of days: - month_length = mdata.time.dt.days_in_month - weights = (month_length.groupby("time.season") / month_length.groupby("time.season").sum()) + mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) + oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) + else: + #Just average months as-is: + mseasons[s] = mdata.sel(time=seasons[s]).mean(dim='time') + oseasons[s] = odata.sel(time=seasons[s]).mean(dim='time') #End if - #Loop over pressure levels: - for pres in pres_levs: - - #Check that the user-requested pressure level - #exists in the model data, which should already - #have been interpolated to the standard reference - #pressure levels: - if not (pres in mclim_ds['lev']): - #Move on to the next pressure level: - print(f"plot_press_levels value '{pres}' not a standard reference pressure, so skipping.") - continue - #End if - - #Create new dictionaries: - mseasons = {} - oseasons = {} - dseasons = {} - - #Loop over seasons: - for s in seasons: - plot_name = plot_loc / f"{var}_{pres}hpa_{s}_LatLon_Mean.{plot_type}" - - # Check redo_plot. If set to True: remove old plot, if it already exists: - redo_plot = adfobj.get_basic_info('redo_plot') - 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, f"{var}_{pres}hpa", case_name, category=web_category, - season=s, plot_type="LatLon") - - #Continue to next iteration: - continue - elif (redo_plot) and plot_name.is_file(): - plot_name.unlink() - - #If requested, then calculate the monthly-weighted seasonal averages: - if weight_season: - mseasons[s] = (pf.seasonal_mean(mdata, season=s, is_climo=True)).sel(lev=pres) - oseasons[s] = (pf.seasonal_mean(odata, season=s, is_climo=True)).sel(lev=pres) - else: - #Just average months as-is: - mseasons[s] = mdata.sel(time=seasons[s], lev=pres).mean(dim='time') - oseasons[s] = odata.sel(time=seasons[s], lev=pres).mean(dim='time') - #End if - - # difference: each entry should be (lat, lon) - dseasons[s] = mseasons[s] - oseasons[s] - - # time to make plot; here we'd probably loop over whatever plots we want for this variable - # I'll just call this one "LatLon_Mean" ... would this work as a pattern [operation]_[AxesDescription] ? - - #Create new plot: - # NOTE: send vres as kwarg dictionary. --> ONLY vres, not the full res - # This relies on `plot_map_and_save` knowing how to deal with the options - # currently knows how to handle: - # colormap, contour_levels, diff_colormap, diff_contour_levels, tiString, tiFontSize, mpl - # *Any other entries will be ignored. - # NOTE: If we were doing all the plotting here, we could use whatever we want from the provided YAML file. - pf.plot_map_and_save(plot_name, case_nickname, base_nickname, - [syear_cases[case_idx],eyear_cases[case_idx]], - [syear_baseline,eyear_baseline], - mseasons[s], oseasons[s], dseasons[s], - obs, **vres) - - #Add plot to website (if enabled): - adfobj.add_website_data(plot_name, f"{var}_{pres}hpa", case_name, category=web_category, - season=s, plot_type="LatLon") - - #End for (seasons) - #End for (pressure levels) - - else: - print(f"\t - variable '{var}' has no vertical dimension but is not just time/lat/lon, so skipping.") - #End if (has_lev) - else: - print(f"\t - skipping polar map for {var} as it has more than lat/lon dims, but no pressure levels were provided") - #End if (dimensions check and plotting pressure levels) - #End for (case loop) - #End for (obs/baseline loop) + # difference: each entry should be (lat, lon) + dseasons[s] = mseasons[s] - oseasons[s] + + pf.plot_map_and_save(plot_name, case_nickname, data.ref_nickname, + [syear_cases[case_idx],eyear_cases[case_idx]], + [syear_baseline,eyear_baseline], + mseasons[s].sel(lev=pres), oseasons[s].sel(lev=pres), dseasons[s].sel(lev=pres), + obs=adfobj.compare_obs, **vres) + + #Add plot to website (if enabled): + adfobj.add_website_data(plot_name, f"{var}_{pres}hpa", case_name, category=web_category, + season=s, plot_type="LatLon") + #End for (seasons) + #End for (pressure levels) + #End if (plotting pressure levels) + #End for (case loop) #End for (variable loop) #Notify user that script has ended: print(" ...lat/lon maps have been generated successfully.") -######### -# Helpers -######### +def plot_file_op(adfobj, plot_name, var, case_name, season, web_category, redo_plot, plot_type): + """Check if output plot needs to be made or remade. + + Parameters + ---------- + adfobj : AdfDiag + The diagnostics object that contains all the configuration information + + plot_name : Path + path of the output plot + + var : str + name of variable + + case_name : str + case name + + season : str + season being plotted + + web_category : str + the category for this variable + + redo_plot : bool + whether to overwrite existing plot with this file name + + plot_type : str + the file type for the output plot + + Returns + ------- + int, None + Returns 1 if existing file is removed or no existing file. + Returns None if file exists and redo_plot is False + + Notes + ----- + The long list of parameters is because add_website_data is called + when the file exists and will not be overwritten. + + """ + # Check redo_plot. If set to True: remove old plot, if it already exists: + if plot_name.is_file(): + if redo_plot: + plot_name.unlink() + return True + else: + #Add already-existing plot to website (if enabled): + adfobj.add_website_data(plot_name, var, case_name, category=web_category, + season=season, plot_type=plot_type) + return False # False tells caller that file exists and not to overwrite + else: + return True ############## -#END OF SCRIPT \ No newline at end of file +#END OF SCRIPT diff --git a/scripts/plotting/global_latlon_map_B.py b/scripts/plotting/global_latlon_map_B.py deleted file mode 100644 index 2e13d62d5..000000000 --- a/scripts/plotting/global_latlon_map_B.py +++ /dev/null @@ -1,354 +0,0 @@ -""" -Generate global maps of 2-D fields - -Functions ---------- -global_latlon_map(adfobj) - use ADF object to make maps -my_formatwarning(msg, *args, **kwargs) - format warning messages - (private method) -plot_file_op - Check on status of output plot file. -""" -#Import standard modules: -from pathlib import Path -import numpy as np -import xarray as xr -import warnings # use to warn user about missing files. - -import plotting_functions as pf - -#Format warning messages: -def my_formatwarning(msg, *args, **kwargs): - """Issue `msg` as warning.""" - return str(msg) + '\n' -warnings.formatwarning = my_formatwarning - -def global_latlon_map_B(adfobj): - """ - This script/function is designed to generate global - 2-D lat/lon maps of model fields with continental overlays. - - Parameters - ---------- - adfobj : AdfDiag - The diagnostics object that contains all the configuration information - - Returns - ------- - Does not return a value; produces plots and saves files. - - Notes - ----- - - It uses the AdfDiag object's methods to get necessary information. - This version passes the AdfDiag object to instatiate an AdfData object. - Explicitly accesses: - adfobj.diag_var_list - List of variables - adfobj.plot_location - output plot path - adfobj.climo_yrs - start and end climo years of the case(s), `syears` & `eyears` - start and end climo years of the reference, `syear_baseline` & `eyear_baseline` - adfobj.variable_defaults - dict of variable-specific plot preferences - adfobj.read_config_var - dict of basic info, `diag_basic_info` - Then use to check `plot_type` - adfobj.debug_log - Issues debug message - adfobj.add_website_data - Communicates information to the website generator - - - The `plotting_functions` module is needed for: - pf.get_central_longitude() - determine central longitude for global plots - pf.lat_lon_validate_dims() - makes sure latitude and longitude are valid - pf.seasonal_mean() - calculate seasonal mean - pf.plot_map_and_save() - send information to make the plot and save the file - pf.zm_validate_dims() - Checks on pressure level dimension - """ - - #Notify user that script has started: - print("\n Generating lat/lon maps...") - - # - # Use ADF api to get all necessary information - # - var_list = adfobj.diag_var_list - #Special ADF variable which contains the output paths for - #all generated plots and tables for each case: - plot_locations = adfobj.plot_location - - #Grab case years - syear_cases = adfobj.climo_yrs["syears"] - eyear_cases = adfobj.climo_yrs["eyears"] - - #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"] - - res = adfobj.variable_defaults # will be dict of variable-specific plot preferences - # 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 - # -- So check for it, and default to png - basic_info_dict = adfobj.read_config_var("diag_basic_info") - plot_type = basic_info_dict.get('plot_type', 'png') - print(f"\t NOTE: Plot type is set to {plot_type}") - - # check if existing plots need to be redone - redo_plot = adfobj.get_basic_info('redo_plot') - print(f"\t NOTE: redo_plot is set to {redo_plot}") - #----------------------------------------- - - #Determine if user wants to plot 3-D variables on - #pressure levels: - pres_levs = adfobj.get_basic_info("plot_press_levels") - - weight_season = True #always do seasonal weighting - - #Set seasonal ranges: - seasons = {"ANN": np.arange(1,13,1), - "DJF": [12, 1, 2], - "JJA": [6, 7, 8], - "MAM": [3, 4, 5], - "SON": [9, 10, 11] - } - - # probably want to do this one variable at a time: - for var in var_list: - if var not in adfobj.data.ref_var_nam: - dmsg = f"No reference data found for variable `{var}`, zonal mean plotting skipped." - adfobj.debug_log(dmsg) - continue - - #Notify user of variable being plotted: - print("\t - lat/lon maps for {}".format(var)) - - # 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"global_latlon_map: Found variable defaults for {var}") - - #Extract category (if available): - web_category = vres.get("category", None) - - else: - vres = {} - web_category = None - #End if - - # For global maps, also set the central longitude: - # can be specified in adfobj basic info as 'central_longitude' or supplied as a number, - # otherwise defaults to 180 - vres['central_longitude'] = pf.get_central_longitude(adfobj) - - # load reference data (observational or baseline) - odata = adfobj.data.load_reference_da(var) - if odata is None: - continue - has_dims = pf.lat_lon_validate_dims(odata) # T iff dims are (lat,lon) -- can't plot unless we have both - if not has_dims: - print(f"\t = skipping global map for {var} as REFERENCE does not have both lat and lon") - continue - - #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]) - - #Check if plot output directory exists, and if not, then create it: - if not plot_loc.is_dir(): - print(" {} not found, making new directory".format(plot_loc)) - plot_loc.mkdir(parents=True) - - #Load re-gridded model files: - mdata = adfobj.data.load_regrid_da(case_name, var) - - #Skip this variable/case if the regridded climo file doesn't exist: - if mdata is None: - continue - - #Determine dimensions of variable: - has_dims_cam = pf.lat_lon_validate_dims(mdata) # T iff dims are (lat,lon) -- can't plot unless we have both - _, has_lev = pf.zm_validate_dims(mdata) # has_lev T if lev in mdata - if not has_dims_cam: - print(f"\t = skipping global map for {var} for case {case_name} as it does not have both lat and lon") - continue - else: # i.e., has lat&lon - if pres_levs and (not has_lev): - print(f"\t - skipping global map for {var} as it has more than lat/lon dims, but no pressure levels were provided") - continue - - # Check output file. If file does not exist, proceed. - # If file exists: - # if redo_plot is true: delete it now and make plot - # if redo_plot is false: add to website and move on - doplot = {} - if not pres_levs: - for s in seasons: - plot_name = plot_loc / f"{var}_{s}_LatLon_Mean.{plot_type}" - doplot[plot_name] = plot_file_op(adfobj, plot_name, var, case_name, s, web_category, redo_plot, "LatLon") - else: - for pres in pres_levs: - for s in seasons: - plot_name = plot_loc / f"{var}_{pres}hpa_{s}_LatLon_Mean.{plot_type}" - doplot[plot_name] = plot_file_op(adfobj, plot_name, f"{var}_{pres}hpa", case_name, s, web_category, redo_plot, "LatLon") - if all(value is None for value in doplot.values()): - print(f"All plots exist for {var}. Redo is {redo_plot}. Existing plots added to website data. Continue.") - continue - - #Create new dictionaries: - mseasons = {} - oseasons = {} - dseasons = {} # hold the differences - - if not pres_levs: - - #Loop over season dictionary: - for s in seasons: - plot_name = plot_loc / f"{var}_{s}_LatLon_Mean.{plot_type}" - if doplot[plot_name] is None: - continue - - if weight_season: - mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) - oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) - else: - #Just average months as-is: - mseasons[s] = mdata.sel(time=seasons[s]).mean(dim='time') - oseasons[s] = odata.sel(time=seasons[s]).mean(dim='time') - #End if - - # difference: each entry should be (lat, lon) - dseasons[s] = mseasons[s] - oseasons[s] - - pf.plot_map_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], dseasons[s], - obs=adfobj.compare_obs, **vres) - - #Add plot to website (if enabled): - adfobj.add_website_data(plot_name, var, case_name, category=web_category, - season=s, plot_type="LatLon") - - else: # => pres_levs has values, & we already checked that lev is in mdata (has_lev) - - for pres in pres_levs: - - #Check that the user-requested pressure level - #exists in the model data, which should already - #have been interpolated to the standard reference - #pressure levels: - if not (pres in mdata['lev']): - print(f"plot_press_levels value '{pres}' not present in {var}, so skipping.") - continue - - #Loop over seasons: - for s in seasons: - plot_name = plot_loc / f"{var}_{pres}hpa_{s}_LatLon_Mean.{plot_type}" - if doplot[plot_name] is None: - continue - - if weight_season: - mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) - oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) - else: - #Just average months as-is: - mseasons[s] = mdata.sel(time=seasons[s]).mean(dim='time') - oseasons[s] = odata.sel(time=seasons[s]).mean(dim='time') - #End if - - # difference: each entry should be (lat, lon) - dseasons[s] = mseasons[s] - oseasons[s] - - pf.plot_map_and_save(plot_name, case_nickname, data.ref_nickname, - [syear_cases[case_idx],eyear_cases[case_idx]], - [syear_baseline,eyear_baseline], - mseasons[s].sel(lev=pres), oseasons[s].sel(lev=pres), dseasons[s].sel(lev=pres), - obs=adfobj.compare_obs, **vres) - - #Add plot to website (if enabled): - adfobj.add_website_data(plot_name, f"{var}_{pres}hpa", case_name, category=web_category, - season=s, plot_type="LatLon") - #End for (seasons) - #End for (pressure levels) - #End if (plotting pressure levels) - #End for (case loop) - #End for (variable loop) - - #Notify user that script has ended: - print(" ...lat/lon maps have been generated successfully.") - - -def plot_file_op(adfobj, plot_name, var, case_name, season, web_category, redo_plot, plot_type): - """Check if output plot needs to be made or remade. - - Parameters - ---------- - adfobj : AdfDiag - The diagnostics object that contains all the configuration information - - plot_name : Path - path of the output plot - - var : str - name of variable - - case_name : str - case name - - season : str - season being plotted - - web_category : str - the category for this variable - - redo_plot : bool - whether to overwrite existing plot with this file name - - plot_type : str - the file type for the output plot - - Returns - ------- - int, None - Returns 1 if existing file is removed or no existing file. - Returns None if file exists and redo_plot is False - - Notes - ----- - The long list of parameters is because add_website_data is called - when the file exists and will not be overwritten. - - """ - # Check redo_plot. If set to True: remove old plot, if it already exists: - if plot_name.is_file(): - if redo_plot: - plot_name.unlink() - return True - else: - #Add already-existing plot to website (if enabled): - adfobj.add_website_data(plot_name, var, case_name, category=web_category, - season=season, plot_type=plot_type) - return False # False tells caller that file exists and not to overwrite - else: - return True - -############## -#END OF SCRIPT diff --git a/scripts/plotting/zonal_mean.py b/scripts/plotting/zonal_mean.py index 944683752..019065d9f 100644 --- a/scripts/plotting/zonal_mean.py +++ b/scripts/plotting/zonal_mean.py @@ -2,6 +2,7 @@ import numpy as np import xarray as xr import plotting_functions as pf + import warnings # use to warn user about missing files. def my_formatwarning(msg, *args, **kwargs): @@ -10,91 +11,53 @@ def my_formatwarning(msg, *args, **kwargs): warnings.formatwarning = my_formatwarning -def zonal_mean(adfobj): +def zonal_mean_B(adfobj): """ - This script plots zonal averages. + Plots zonal average from climatological files (annual and seasonal). Compare CAM climatologies against other climatological data (observations or baseline runs). - Description of needed inputs from ADF: - case_name -> Name of CAM case provided by "cam_case_name". - model_rgrid_loc -> Location of re-gridded CAM climo files provided by "cam_regrid_loc". - data_name -> Name of data set CAM case is being compared against, - which is always either "obs" or the baseline CAM case name, - depending on whether "compare_obs" is true or false. - data_loc -> Location of comparison data, which is either "obs_climo_loc" - or "cam_baseline_climo_loc", depending on whether - "compare_obs" is true or false. - var_list -> List of CAM output variables provided by "diag_var_list" - data_list -> List of data sets CAM will be compared against, which - is simply the baseline case name in situations when - "compare_obs" is false. - climo_yrs -> Dictionary containing the start and end years of the test - and baseline model data (if applicable). - plot_location -> Location where plot files will be written to, which is - specified by "cam_diag_plot_loc". - variable_defaults -> optional, - Dict that has keys that are variable names and values that are plotting preferences/defaults. - Notes: - The script produces plots of 2-D and 3-D variables, - but needs to determine which type along the way. - For 3-D variables, the default behavior is to interpolate - climo files to pressure levels, which requires the hybrid-sigma - coefficients and surface pressure. That ASSUMES that the climo - files are using native hybrid-sigma levels rather than being - transformed to pressure levels. + + 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: + diag_var_list, climo_yrs, variable_defaults, read_config_var, + get_basic_info, add_website_data, debug_log + + Determines whether `lev` dimension is present. If not, makes + a line plot, but if so it makes a contour plot. + TODO: There's a flag to plot linear vs log pressure, but no + method to infer what the user wants. """ - #Notify user that script has started: print("\n Generating zonal mean plots...") - #Extract needed quantities from ADF object: - #----------------------------------------- var_list = adfobj.diag_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 zonal-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 YAML. @@ -110,12 +73,6 @@ def zonal_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), @@ -131,7 +88,7 @@ def zonal_mean(adfobj): logp_zonal_skip = [] #Loop over model cases: - for case_idx, case_name in enumerate(case_names): + for case_idx, case_name in enumerate(adfobj.data.case_names): #Set output plot location: plot_loc = Path(plot_locations[case_idx]) @@ -145,7 +102,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}_logp_{s}_Zonal_Mean.{plot_type}" + plot_name_log = plot_loc / f"{var}_{s}_Zonal_logp_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(): @@ -165,7 +122,6 @@ def zonal_mean(adfobj): if (not redo_plot) and plot_name.is_file(): zonal_skip.append(plot_name) #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="Zonal") @@ -185,27 +141,10 @@ def zonal_mean(adfobj): # #Loop over variables: for var in var_list: - - 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}`, zonal mean plotting skipped." - adfobj.debug_log(dmsg) - continue - #End if - else: - #Set "data_var" for consistent use below: - data_var = var - #End if + if var not in adfobj.data.ref_var_nam: + dmsg = f"No obs found for variable `{var}`, zonal mean plotting skipped." + adfobj.debug_log(dmsg) + continue #Notify user of variable being plotted: print(f"\t - zonal mean plots for {var}") @@ -220,81 +159,48 @@ def zonal_mean(adfobj): 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")) - #End if - oclim_ds = pf.load_dataset(oclim_fils) - - #Loop over model cases: - for case_idx, case_name in enumerate(case_names): - - #Set case nickname: - case_nickname = test_nicknames[case_idx] - - #Set output plot location: - plot_loc = Path(plot_locations[case_idx]) - - # load re-gridded model files: - mclim_fils = sorted(mclimo_rg_loc.glob(f"{data_src}_{case_name}_{var}_*.nc")) - mclim_ds = pf.load_dataset(mclim_fils) - - # stop if data is invalid: - if (oclim_ds is None) or (mclim_ds is None): - warnings.warn(f"invalid data, skipping zonal mean plot of {var}") - continue + # load reference data (observational or baseline) + odata = adfobj.data.load_reference_da(var) + has_lat_ref, has_lev_ref = pf.zm_validate_dims(odata) - #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. + #Loop over model cases: + for case_idx, case_name in enumerate(adfobj.data.case_names): - # determine whether it's 2D or 3D - # 3D triggers search for surface pressure - has_lat, has_lev = pf.zm_validate_dims(mdata) # assumes will work for both mdata & odata + #Set case nickname: + case_nickname = adfobj.data.test_nicknames[case_idx] - #Notify user of level dimension: - if has_lev: - print(f"\t {var} has lev dimension.") + #Set output plot location: + plot_loc = Path(plot_locations[case_idx]) - # - # Seasonal Averages - # + # load re-gridded model files: + mdata = adfobj.data.load_regrid_da(case_name, var) - #Create new dictionaries: - mseasons = {} - oseasons = {} - - #Loop over season dictionary: - for s in seasons: - - # time to make plot; here we'd probably loop over whatever plots we want for this variable - # I'll just call this one "Zonal_Mean" ... would this work as a pattern [operation]_[AxesDescription] ? - # NOTE: Up to this point, nothing really differs from global_latlon_map, - # so we could have made one script instead of two. - # Merging would make overall timing better because looping twice will double I/O steps. - # + # determine whether it's 2D or 3D + # 3D triggers search for surface pressure + has_lat, has_lev = pf.zm_validate_dims(mdata) # assumes will work for both mdata & odata + + #Notify user of level dimension: + if has_lev: + print(f"\t {var} has lev dimension.") + + # + # Seasonal Averages + # + + #Create new dictionaries: + mseasons = {} + oseasons = {} + + #Loop over season dictionary: + for s in seasons: + + # time to make plot; here we'd probably loop over whatever plots we want for this variable + # I'll just call this one "Zonal_Mean" ... would this work as a pattern [operation]_[AxesDescription] ? + # NOTE: Up to this point, nothing really differs from global_latlon_map, + # so we could have made one script instead of two. + # Merging would make overall timing better because looping twice will double I/O steps. + # + if not has_lev: plot_name = plot_loc / f"{var}_{s}_Zonal_Mean.{plot_type}" if plot_name not in zonal_skip: @@ -302,7 +208,7 @@ def zonal_mean(adfobj): #Seasonal Averages mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) - + # difference: each entry should be (lat, lon) or (plev, lat, lon) # dseasons[s] = mseasons[s] - oseasons[s] # difference will be calculated in plot_zonal_mean_and_save; @@ -310,40 +216,45 @@ def zonal_mean(adfobj): # This could be re-visited for efficiency or improved code structure. #Create new plot: - pf.plot_zonal_mean_and_save(plot_name, case_nickname, base_nickname, + pf.plot_zonal_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, log_p=False, obs=obs, **vres) + mseasons[s], oseasons[s], has_lev, log_p=False, obs=adfobj.compare_obs, **vres) #Add plot to website (if enabled): adfobj.add_website_data(plot_name, var, case_name, season=s, plot_type="Zonal") - #Create new plot with log-p: - if has_lev: - plot_name_log = plot_loc / f"{var}_logp_{s}_Zonal_Mean.{plot_type}" + #Create new plot with log-p: + # NOTE: The log-p should be an option here. + else: + if (not has_lev_ref) or (not has_lev): + print(f"Error: expecting lev for both case: {has_lev} and ref: {has_lev_ref}") + continue + if len(mdata['lev']) != len(odata['lev']): + print(f"Error: zonal mean contour expects `lev` dim to have same size, got {len(mdata['lev'])} and {len(odata['lev'])}") + continue + plot_name_log = plot_loc / f"{var}_{s}_Zonal_logp_Mean.{plot_type}" + if plot_name_log not in logp_zonal_skip: + #Seasonal Averages + mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) + oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) - if plot_name_log not in logp_zonal_skip: - pf.plot_zonal_mean_and_save(plot_name_log, case_nickname, base_nickname, - [syear_cases[case_idx],eyear_cases[case_idx]], - [syear_baseline,eyear_baseline], - mseasons[s], oseasons[s], has_lev, log_p=True, obs=obs, **vres) + pf.plot_zonal_mean_and_save(plot_name_log, case_nickname, adfobj.data.ref_nickname, + [syear_cases[case_idx],eyear_cases[case_idx]], + [syear_baseline,eyear_baseline], + mseasons[s], oseasons[s], has_lev, log_p=True, obs=adfobj.compare_obs, **vres) - #Add plot to website (if enabled): - adfobj.add_website_data(plot_name_log, f"{var}_logp", case_name, season=s, plot_type="Zonal", category="Log-P") + #Add plot to website (if enabled): + adfobj.add_website_data(plot_name_log, f"{var}_logp", case_name, season=s, plot_type="Zonal", category="Log-P") - #End for (seasons loop) - #End for (case names loop) - #End for (obs/baseline loop) + #End for (seasons loop) + #End for (case names loop) #End for (variables loop) #Notify user that script has ended: print(" ...Zonal mean plots have been generated successfully.") -######### -# Helpers -######### - - ############## -#END OF SCRIPT \ No newline at end of file +#END OF SCRIPT + diff --git a/scripts/plotting/zonal_mean_B.py b/scripts/plotting/zonal_mean_B.py deleted file mode 100644 index 019065d9f..000000000 --- a/scripts/plotting/zonal_mean_B.py +++ /dev/null @@ -1,260 +0,0 @@ -from pathlib import Path -import numpy as np -import xarray as xr -import plotting_functions as pf - -import warnings # use to warn user about missing files. - -def my_formatwarning(msg, *args, **kwargs): - # ignore everything except the message - return str(msg) + '\n' - -warnings.formatwarning = my_formatwarning - -def zonal_mean_B(adfobj): - - """ - Plots zonal average from climatological files (annual and seasonal). - Compare 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: - diag_var_list, climo_yrs, variable_defaults, read_config_var, - get_basic_info, add_website_data, debug_log - - Determines whether `lev` dimension is present. If not, makes - a line plot, but if so it makes a contour plot. - TODO: There's a flag to plot linear vs log pressure, but no - method to infer what the user wants. - """ - - print("\n Generating zonal mean plots...") - - var_list = adfobj.diag_var_list - - #Special ADF variable which contains the output paths for - #all generated plots and tables: - plot_locations = adfobj.plot_location - - #Grab case years - syear_cases = adfobj.climo_yrs["syears"] - eyear_cases = adfobj.climo_yrs["eyears"] - - #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"] - - res = adfobj.variable_defaults # will be dict of variable-specific plot preferences - # 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 - # -- So check for it, and default to png - basic_info_dict = adfobj.read_config_var("diag_basic_info") - plot_type = basic_info_dict.get('plot_type', 'png') - print(f"\t NOTE: Plot type is set to {plot_type}") - - # check if existing plots need to be redone - redo_plot = adfobj.get_basic_info('redo_plot') - print(f"\t NOTE: redo_plot is set to {redo_plot}") - #----------------------------------------- - - - #Set seasonal ranges: - seasons = {"ANN": np.arange(1,13,1), - "DJF": [12, 1, 2], - "JJA": [6, 7, 8], - "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 attempt to skip calcs to - #speed up preformance a bit if re-running the ADF - zonal_skip = [] - logp_zonal_skip = [] - - #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: - #Check zonal log-p: - plot_name_log = plot_loc / f"{var}_{s}_Zonal_logp_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(): - logp_zonal_skip.append(plot_name_log) - #Continue to next iteration: - adfobj.add_website_data(plot_name_log, f"{var}_logp", case_name, season=s, - plot_type="Zonal", category="Log-P") - pass - - elif (redo_plot) and plot_name_log.is_file(): - plot_name_log.unlink() - #End if - - #Check regular zonal - plot_name = plot_loc / f"{var}_{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.is_file(): - zonal_skip.append(plot_name) - #Add already-existing plot to website (if enabled): - adfobj.add_website_data(plot_name, var, case_name, season=s, - plot_type="Zonal") - - 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: - if var not in adfobj.data.ref_var_nam: - dmsg = f"No obs found for variable `{var}`, zonal mean plotting skipped." - adfobj.debug_log(dmsg) - continue - - #Notify user of variable being plotted: - print(f"\t - zonal mean plots for {var}") - - # 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"zonal_mean: Found variable defaults for {var}") - - else: - vres = {} - #End if - - # load reference data (observational or baseline) - odata = adfobj.data.load_reference_da(var) - has_lat_ref, has_lev_ref = pf.zm_validate_dims(odata) - - #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]) - - # load re-gridded model files: - mdata = adfobj.data.load_regrid_da(case_name, var) - - # determine whether it's 2D or 3D - # 3D triggers search for surface pressure - has_lat, has_lev = pf.zm_validate_dims(mdata) # assumes will work for both mdata & odata - - #Notify user of level dimension: - if has_lev: - print(f"\t {var} has lev dimension.") - - # - # Seasonal Averages - # - - #Create new dictionaries: - mseasons = {} - oseasons = {} - - #Loop over season dictionary: - for s in seasons: - - # time to make plot; here we'd probably loop over whatever plots we want for this variable - # I'll just call this one "Zonal_Mean" ... would this work as a pattern [operation]_[AxesDescription] ? - # NOTE: Up to this point, nothing really differs from global_latlon_map, - # so we could have made one script instead of two. - # Merging would make overall timing better because looping twice will double I/O steps. - # - if not has_lev: - plot_name = plot_loc / f"{var}_{s}_Zonal_Mean.{plot_type}" - - if plot_name not in zonal_skip: - - #Seasonal Averages - mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) - oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) - - # difference: each entry should be (lat, lon) or (plev, lat, lon) - # dseasons[s] = mseasons[s] - oseasons[s] - # difference will be calculated in plot_zonal_mean_and_save; - # because we can let any pressure-level interpolation happen there - # This could be re-visited for efficiency or improved code structure. - - #Create new plot: - pf.plot_zonal_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, log_p=False, obs=adfobj.compare_obs, **vres) - - #Add plot to website (if enabled): - adfobj.add_website_data(plot_name, var, case_name, season=s, plot_type="Zonal") - - #Create new plot with log-p: - # NOTE: The log-p should be an option here. - else: - if (not has_lev_ref) or (not has_lev): - print(f"Error: expecting lev for both case: {has_lev} and ref: {has_lev_ref}") - continue - if len(mdata['lev']) != len(odata['lev']): - print(f"Error: zonal mean contour expects `lev` dim to have same size, got {len(mdata['lev'])} and {len(odata['lev'])}") - continue - plot_name_log = plot_loc / f"{var}_{s}_Zonal_logp_Mean.{plot_type}" - if plot_name_log not in logp_zonal_skip: - #Seasonal Averages - mseasons[s] = pf.seasonal_mean(mdata, season=s, is_climo=True) - oseasons[s] = pf.seasonal_mean(odata, season=s, is_climo=True) - - pf.plot_zonal_mean_and_save(plot_name_log, case_nickname, adfobj.data.ref_nickname, - [syear_cases[case_idx],eyear_cases[case_idx]], - [syear_baseline,eyear_baseline], - mseasons[s], oseasons[s], has_lev, log_p=True, obs=adfobj.compare_obs, **vres) - - #Add plot to website (if enabled): - adfobj.add_website_data(plot_name_log, f"{var}_logp", case_name, season=s, plot_type="Zonal", category="Log-P") - - #End for (seasons loop) - #End for (case names loop) - #End for (variables loop) - - #Notify user that script has ended: - print(" ...Zonal mean plots have been generated successfully.") - - -############## -#END OF SCRIPT - From df14b1b3934598a3cd4024aa3ecd7847cd00dbd2 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 21 Jun 2024 13:35:52 -0600 Subject: [PATCH 40/46] try to merge amwg_table from ADF/main --- scripts/analysis/amwg_table.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/scripts/analysis/amwg_table.py b/scripts/analysis/amwg_table.py index 794d5a678..b2bac7332 100644 --- a/scripts/analysis/amwg_table.py +++ b/scripts/analysis/amwg_table.py @@ -166,10 +166,6 @@ def amwg_table(adf): raise AdfError(errmsg) #Write to debug log if enabled: adf.debug_log(f"DEBUG: location of files is {str(input_location)}") - #Check if analysis directory exists, and if not, then create it: - if not output_location.is_dir(): - print(f"\t {output_locs[case_idx]} not found, making new directory") - output_location.mkdir(parents=True) #Create output file name: output_csv_file = output_location / f"amwg_table_{case_name}.csv" From b966637ad69630bab5c6bcbe93e0dac56c08d563 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 21 Jun 2024 13:42:53 -0600 Subject: [PATCH 41/46] resolve upstream conflicts --- lib/adf_info.py | 39 +++++++++++++-------------- scripts/plotting/tape_recorder.py | 45 ++++++++++++++++++++----------- 2 files changed, 48 insertions(+), 36 deletions(-) diff --git a/lib/adf_info.py b/lib/adf_info.py index 7fed436b7..cdde4d6f3 100644 --- a/lib/adf_info.py +++ b/lib/adf_info.py @@ -111,6 +111,7 @@ def __init__(self, config_file, debug=False): if not hist_str: hist_str = 'cam.h0' #End if + self.__hist_str = hist_str #Initialize ADF variable list: self.__diag_var_list = self.read_config_var('diag_var_list', required=True) @@ -309,10 +310,10 @@ def __init__(self, config_file, debug=False): #Extract cam history files location: cam_hist_locs = self.get_cam_info('cam_hist_loc') - + #Check if using pre-made ts files cam_ts_done = self.get_cam_info("cam_ts_done") - + #Grab case time series file location(s) input_ts_locs = self.get_cam_info("cam_ts_loc", required=True) @@ -416,13 +417,19 @@ def __init__(self, config_file, debug=False): #Set the final directory name and save it to plot_location: direc_name = f"{case_name}_vs_{data_name}" - self.__plot_location.append(os.path.join(plot_dir, direc_name)) + plot_loc = os.path.join(plot_dir, direc_name) + self.__plot_location.append(plot_loc) #If first iteration, then save directory name for use by baseline: if case_idx == 0: first_case_dir = direc_name #End if + #Go ahead and make the diag plot location if it doesn't exist already + diag_location = Path(plot_loc) + if not diag_location.is_dir(): + print(f"\t {diag_location} not found, making new directory") + diag_location.mkdir(parents=True) #End for self.__syears = syears_fixed @@ -577,6 +584,11 @@ def case_nicknames(self): return {"test_nicknames":test_nicknames,"base_nickname":base_nickname} + @property + def hist_string(self): + """ Return the history string name to the user if requested.""" + return self.__hist_str + ######### #Utility function to access expanded 'diag_basic_info' variables: @@ -673,15 +685,7 @@ def get_climo_yrs_from_ts(self, input_ts_loc, case_name): #Search for first variable in var_list to get a time series file to read #NOTE: it is assumed all the variables have the same dates! - - #Read hist_str (component.hist_num) from the yaml file, or set to default - hist_str = self.get_basic_info('hist_str') - #If hist_str is not present, then default to 'cam.h0': - if not hist_str: - hist_str = 'cam.h0' - #End if - - ts_files = sorted(input_location.glob(f"{case_name}.{hist_str}.{var_list[0]}.*nc")) + ts_files = sorted(input_location.glob(f"{case_name}*.{var_list[0]}.*nc")) #Read in file(s) if len(ts_files) == 1: @@ -691,16 +695,11 @@ def get_climo_yrs_from_ts(self, input_ts_loc, case_name): #Average time dimension over time bounds, if bounds exist: if 'time_bnds' in cam_ts_data: - timeBoundsName = 'time_bnds' - elif 'time_bounds' in cam_ts_data: - timeBoundsName = 'time_bounds' - else: - timeBoundsName = None - if timeBoundsName: time = cam_ts_data['time'] - #NOTE: force `load` here b/c if dask & time is cftime, + #NOTE: force `load` here b/c if dask & time is cftime, #throws a NotImplementedError: - time = xr.DataArray(cam_ts_data[timeBoundsName].load().mean(dim='nbnd').values, dims=time.dims, attrs=time.attrs) + time = xr.DataArray(cam_ts_data['time_bnds'].load().mean(dim='nbnd').values, + dims=time.dims, attrs=time.attrs) cam_ts_data['time'] = time cam_ts_data.assign_coords(time=time) cam_ts_data = xr.decode_cf(cam_ts_data) diff --git a/scripts/plotting/tape_recorder.py b/scripts/plotting/tape_recorder.py index fe048dc5d..2ec3048e0 100644 --- a/scripts/plotting/tape_recorder.py +++ b/scripts/plotting/tape_recorder.py @@ -36,6 +36,9 @@ def tape_recorder(adfobj): plot_location = adfobj.plot_location plot_loc = Path(plot_location[0]) + #Grab history string: + hist_str = adfobj.hist_string + #Grab test case name(s) case_names = adfobj.get_cam_info('cam_case_name', required=True) @@ -47,9 +50,7 @@ def tape_recorder(adfobj): end_years = adfobj.climo_yrs["eyears"] #Grab test case nickname(s) - test_nicknames = adfobj.get_cam_info('case_nickname') - if test_nicknames == None: - test_nicknames = case_names + test_nicknames = adfobj.case_nicknames['test_nicknames'] # CAUTION: # "data" here refers to either obs or a baseline simulation, @@ -64,9 +65,7 @@ def tape_recorder(adfobj): data_ts_loc = adfobj.get_baseline_info("cam_ts_loc", required=True) case_ts_locs = case_ts_locs+[data_ts_loc] - base_nickname = adfobj.get_baseline_info('case_nickname') - if base_nickname == None: - base_nickname = data_name + base_nickname = adfobj.case_nicknames['base_nickname'] test_nicknames = test_nicknames+[base_nickname] data_start_year = adfobj.climo_yrs["syear_baseline"] @@ -88,17 +87,24 @@ def tape_recorder(adfobj): # check if existing plots need to be redone redo_plot = adfobj.get_basic_info('redo_plot') print(f"\t NOTE: redo_plot is set to {redo_plot}") + + #Set location for observations + obs_loc = Path(adfobj.get_basic_info("obs_data_loc")) #----------------------------------------- + #Set var to Q for now + #TODO: add option to look for H2O if Q is not available, and vice-versa + var = "Q" + #This may have to change if other variables are desired in this plot type? - plot_name = plot_loc / f"Q_TapeRecorder_ANN_Special_Mean.{plot_type}" - print(f"\t - Plotting annual tape recorder for Q") + plot_name = plot_loc / f"{var}_TapeRecorder_ANN_Special_Mean.{plot_type}" + print(f"\t - Plotting annual tape recorder for {var}") # 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, "Q_TapeRecorder", None, season="ANN", multi_case=True) + adfobj.add_website_data(plot_name, f"{var}_TapeRecorder", None, season="ANN", multi_case=True) return elif (redo_plot) and plot_name.is_file(): @@ -110,7 +116,7 @@ def tape_recorder(adfobj): runs_LT2[val] = case_ts_locs[i] # MLS data - mls = xr.open_dataset("/glade/campaign/cgd/cas/islas/CAM7validation/MLS/mls_h2o_latNpressNtime_3d_monthly_v5.nc") + mls = xr.open_dataset(obs_loc / "mls_h2o_latNpressNtime_3d_monthly_v5.nc") mls = mls.rename(x='lat', y='lev', t='time') time = pd.date_range("2004-09","2021-11",freq='M') mls['time'] = time @@ -120,17 +126,22 @@ def tape_recorder(adfobj): mls = mls*18.015280/(1e6*28.964) # ERA5 data - era5 = xr.open_dataset("/glade/campaign/cgd/cas/islas/CAM7validation/ERA5/ERA5_Q_10Sto10N_1980to2020.nc") + era5 = xr.open_dataset(obs_loc / "ERA5_Q_10Sto10N_1980to2020.nc") era5 = era5.groupby('time.month').mean('time') + era5_data = era5.Q alldat=[] runname_LT=[] for idx,key in enumerate(runs_LT2): - fils= sorted(Path(runs_LT2[key]).glob('*h0.Q.*.nc')) + fils= sorted(Path(runs_LT2[key]).glob(f'*{hist_str}.{var}.*.nc')) dat = pf.load_dataset(fils) + if not dat: + dmsg = f"No data for `{var}` found in {fils}, case will be skipped in tape recorder plot." + print(dmsg) + continue dat = fixcesmtime(dat,start_years[idx],end_years[idx]) datzm = dat.mean('lon') - dat_tropics = cosweightlat(datzm.Q, -10, 10) + dat_tropics = cosweightlat(datzm[var], -10, 10) dat_mon = dat_tropics.groupby('time.month').mean('time').load() alldat.append(dat_mon) runname_LT.append(key) @@ -149,10 +160,12 @@ def tape_recorder(adfobj): x1[0],x2[0],y1[0],y2[0],cmap=cmap, paxis='lev', taxis='month',climo_yrs="2004-2021") - ax = plot_pre_mon(fig, era5.Q, plot_step,plot_min,plot_max, + ax = plot_pre_mon(fig, era5_data, plot_step,plot_min,plot_max, 'ERA5',x1[1],x2[1],y1[1],y2[1], cmap=cmap, paxis='pre', taxis='month',climo_yrs="1980-2020") + + #Start count at 2 to account for MLS and ERA5 plots above count=2 for irun in np.arange(0,alldat_concat_LT.run.size,1): @@ -181,7 +194,7 @@ def tape_recorder(adfobj): y1_loc = y1[count]-0.03 y2_loc = y1[count]-0.02 - ax = plotcolorbar(fig, plot_step, plot_min, plot_max, 'Q (kg/kg)', + ax = plotcolorbar(fig, plot_step, plot_min, plot_max, f'{var} (kg/kg)', x1_loc, x2_loc, y1_loc, y2_loc, cmap=cmap) @@ -189,7 +202,7 @@ def tape_recorder(adfobj): fig.savefig(plot_name, bbox_inches='tight', facecolor='white') #Add plot to website (if enabled): - adfobj.add_website_data(plot_name, "Q_TapeRecorder", None, season="ANN", multi_case=True) + adfobj.add_website_data(plot_name, f"{var}_TapeRecorder", None, season="ANN", multi_case=True) #Notify user that script has ended: print(" ...Tape recorder plots have been generated successfully.") From 39da27e9ffda09c35a1000276ceeb06f854537ad Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Fri, 21 Jun 2024 15:25:05 -0600 Subject: [PATCH 42/46] add back useful changes to adf_info --- lib/adf_info.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/lib/adf_info.py b/lib/adf_info.py index cdde4d6f3..d20d9c090 100644 --- a/lib/adf_info.py +++ b/lib/adf_info.py @@ -687,6 +687,12 @@ def get_climo_yrs_from_ts(self, input_ts_loc, case_name): #NOTE: it is assumed all the variables have the same dates! ts_files = sorted(input_location.glob(f"{case_name}*.{var_list[0]}.*nc")) + #Read hist_str (component.hist_num) from the yaml file, or set to default + hist_str = self.get_basic_info('hist_str') + #If hist_str is not present, then default to 'cam.h0': + if not hist_str: + hist_str = 'cam.h0' + #Read in file(s) if len(ts_files) == 1: cam_ts_data = xr.open_dataset(ts_files[0], decode_times=True) @@ -695,10 +701,17 @@ def get_climo_yrs_from_ts(self, input_ts_loc, case_name): #Average time dimension over time bounds, if bounds exist: if 'time_bnds' in cam_ts_data: + timeBoundsName = 'time_bnds' + elif 'time_bounds' in cam_ts_data: + timeBoundsName = 'time_bounds' + else: + timeBoundsName = None + + if timeBoundsName: time = cam_ts_data['time'] #NOTE: force `load` here b/c if dask & time is cftime, #throws a NotImplementedError: - time = xr.DataArray(cam_ts_data['time_bnds'].load().mean(dim='nbnd').values, + time = xr.DataArray(cam_ts_data[timeBoundsName].load().mean(dim='nbnd').values, dims=time.dims, attrs=time.attrs) cam_ts_data['time'] = time cam_ts_data.assign_coords(time=time) From 231aaaaa1b3da69926901d9997e5d80c24240099 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Thu, 27 Jun 2024 16:02:59 -0600 Subject: [PATCH 43/46] bug fixes. --- lib/adf_dataset.py | 21 +++++++++++++++++++++ scripts/plotting/global_latlon_map.py | 26 ++++++++++++++------------ scripts/plotting/zonal_mean.py | 2 +- 3 files changed, 36 insertions(+), 13 deletions(-) diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py index 5fd9521bc..a6c0be005 100644 --- a/lib/adf_dataset.py +++ b/lib/adf_dataset.py @@ -119,6 +119,22 @@ def load_reference_da(self, variablename): return da + def load_reference_regrid_dataset(self, case, field): + fils = self.get_ref_regrid_file(case, field) + if not fils: + warnings.warn(f"ERROR: Did not find regrid file(s) for case: {case}, variable: {field}") + return None + return self.load_dataset(fils) + + + def load_reference_regrid_da(self, case, field): + fils = self.get_ref_regrid_file(case, field) + if not fils: + warnings.warn(f"ERROR: Did not find regrid file(s) for case: {case}, variable: {field}") + return None + return self.load_da(fils, field) + + def load_climo_da(self, case, variablename): """Return DataArray from climo file""" fils = self.get_climo_file(case, variablename) @@ -159,6 +175,11 @@ def get_ref_timeseries_file(self, field): return ts_files + def get_ref_regrid_file(self, case, field): + model_rg_loc = Path(self.adf.get_basic_info("cam_regrid_loc", required=True)) + return sorted(model_rg_loc.glob(f"{case}_{field}_*.nc")) + + def get_regrid_file(self, case, field): model_rg_loc = Path(self.adf.get_basic_info("cam_regrid_loc", required=True)) rlbl = self.ref_labels[field] # rlbl = "reference label" = the name of the reference data that defines target grid diff --git a/scripts/plotting/global_latlon_map.py b/scripts/plotting/global_latlon_map.py index 2e13d62d5..6939d98e7 100644 --- a/scripts/plotting/global_latlon_map.py +++ b/scripts/plotting/global_latlon_map.py @@ -43,7 +43,7 @@ def global_latlon_map_B(adfobj): ----- It uses the AdfDiag object's methods to get necessary information. - This version passes the AdfDiag object to instatiate an AdfData object. + Makes use of AdfDiag's data sub-class. Explicitly accesses: adfobj.diag_var_list List of variables @@ -61,6 +61,8 @@ def global_latlon_map_B(adfobj): Issues debug message adfobj.add_website_data Communicates information to the website generator + adfobj.compare_obs + Logical to determine if comparing to observations The `plotting_functions` module is needed for: @@ -154,11 +156,12 @@ def global_latlon_map_B(adfobj): vres['central_longitude'] = pf.get_central_longitude(adfobj) # load reference data (observational or baseline) - odata = adfobj.data.load_reference_da(var) + # odata = adfobj.data.load_reference_da(var) + odata = adfobj.data.load_reference_regrid_da(adfobj.data.ref_case_label, var) if odata is None: continue - has_dims = pf.lat_lon_validate_dims(odata) # T iff dims are (lat,lon) -- can't plot unless we have both - if not has_dims: + o_has_dims = pf.validate_dims(odata, ["lat", "lon", "lev"]) # T iff dims are (lat,lon) -- can't plot unless we have both + if (not o_has_dims['has_lat']) or (not o_has_dims['has_lon']): print(f"\t = skipping global map for {var} as REFERENCE does not have both lat and lon") continue @@ -184,13 +187,12 @@ def global_latlon_map_B(adfobj): continue #Determine dimensions of variable: - has_dims_cam = pf.lat_lon_validate_dims(mdata) # T iff dims are (lat,lon) -- can't plot unless we have both - _, has_lev = pf.zm_validate_dims(mdata) # has_lev T if lev in mdata - if not has_dims_cam: + has_dims = pf.validate_dims(mdata, ["lat", "lon", "lev"]) + if (not has_dims['has_lat']) or (not has_dims['has_lon']): print(f"\t = skipping global map for {var} for case {case_name} as it does not have both lat and lon") continue else: # i.e., has lat&lon - if pres_levs and (not has_lev): + if pres_levs and (not has_dims['has_lev']): print(f"\t - skipping global map for {var} as it has more than lat/lon dims, but no pressure levels were provided") continue @@ -217,7 +219,7 @@ def global_latlon_map_B(adfobj): oseasons = {} dseasons = {} # hold the differences - if not pres_levs: + if not has_dims['has_lev']: # strictly 2-d data #Loop over season dictionary: for s in seasons: @@ -255,8 +257,8 @@ def global_latlon_map_B(adfobj): #exists in the model data, which should already #have been interpolated to the standard reference #pressure levels: - if not (pres in mdata['lev']): - print(f"plot_press_levels value '{pres}' not present in {var}, so skipping.") + if (not (pres in mdata['lev'])) or (not (pres in odata['lev'])): + print(f"plot_press_levels value '{pres}' not present in {var} [test: {(pres in mdata['lev'])}, ref: {pres in odata['lev']}], so skipping.") continue #Loop over seasons: @@ -277,7 +279,7 @@ def global_latlon_map_B(adfobj): # difference: each entry should be (lat, lon) dseasons[s] = mseasons[s] - oseasons[s] - pf.plot_map_and_save(plot_name, case_nickname, data.ref_nickname, + pf.plot_map_and_save(plot_name, case_nickname, adfobj.data.ref_nickname, [syear_cases[case_idx],eyear_cases[case_idx]], [syear_baseline,eyear_baseline], mseasons[s].sel(lev=pres), oseasons[s].sel(lev=pres), dseasons[s].sel(lev=pres), diff --git a/scripts/plotting/zonal_mean.py b/scripts/plotting/zonal_mean.py index 019065d9f..a8d5acdc5 100644 --- a/scripts/plotting/zonal_mean.py +++ b/scripts/plotting/zonal_mean.py @@ -160,7 +160,7 @@ def zonal_mean_B(adfobj): #End if # load reference data (observational or baseline) - odata = adfobj.data.load_reference_da(var) + odata = adfobj.data.load_reference_regrid_da(var) has_lat_ref, has_lev_ref = pf.zm_validate_dims(odata) #Loop over model cases: From 6e20daeb3e1724530c3038c9b53dc57e7a011410 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Thu, 27 Jun 2024 16:05:16 -0600 Subject: [PATCH 44/46] correct arguments for load_reference_regrid_da --- scripts/plotting/zonal_mean.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/plotting/zonal_mean.py b/scripts/plotting/zonal_mean.py index a8d5acdc5..a150009f1 100644 --- a/scripts/plotting/zonal_mean.py +++ b/scripts/plotting/zonal_mean.py @@ -160,7 +160,7 @@ def zonal_mean_B(adfobj): #End if # load reference data (observational or baseline) - odata = adfobj.data.load_reference_regrid_da(var) + odata = adfobj.data.load_reference_regrid_da(adfobj.data.ref_case_label, var) has_lat_ref, has_lev_ref = pf.zm_validate_dims(odata) #Loop over model cases: From e36585f5b606aa8374e3bcf85c8ba71437e270ff Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Thu, 11 Jul 2024 15:00:41 -0600 Subject: [PATCH 45/46] trying to fix linting errors --- lib/adf_diag.py | 105 +++++++++++++++++++++++------------------------- lib/adf_info.py | 17 ++++---- 2 files changed, 59 insertions(+), 63 deletions(-) diff --git a/lib/adf_diag.py b/lib/adf_diag.py index 84efc08b2..0053e278b 100644 --- a/lib/adf_diag.py +++ b/lib/adf_diag.py @@ -1,4 +1,3 @@ - """ Location of the "AdfDiag" object, which is used to store all relevant data and @@ -441,7 +440,6 @@ def call_ncrcat(cmd): # Note: could use `open_mfdataset`, but that can become very slow; # This approach effectively assumes that all files contain the same variables. - # Check what kind of vertical coordinate (if any) is being used for this model run: # ------------------------ if "lev" in hist_file_ds: @@ -488,8 +486,8 @@ def call_ncrcat(cmd): print(wmsg) vert_coord_type = None - # End if (long name) - # End if (vert_coord) + # End if (long name) + # End if (vert_coord) else: # No level dimension found, so assume there is no vertical coordinate: vert_coord_type = None @@ -520,8 +518,8 @@ def call_ncrcat(cmd): diag_var_list = self.diag_var_list # Aerosol Calcs - #-------------- - #Always make sure PMID is made if aerosols are desired in config file + # -------------- + # Always make sure PMID is made if aerosols are desired in config file # Since there's no requirement for `aerosol_zonal_list` to be included, allow it to be absent: azl = res.get("aerosol_zonal_list", []) @@ -531,9 +529,9 @@ def call_ncrcat(cmd): if "T" not in diag_var_list: if any(item in azl for item in diag_var_list): diag_var_list += ["T"] - #End aerosol calcs + # End aerosol calcs - #Initialize dictionary for derived variable with needed list of constituents + # Initialize dictionary for derived variable with needed list of constituents constit_dict = {} for var in diag_var_list: @@ -548,27 +546,27 @@ def call_ncrcat(cmd): constit_errmsg += "\n\tPlease add list of constituents to 'derivable_from' " constit_errmsg += f"for {var} in variable defaults yaml file." - #Check if current variable is a derived quantity + # Check if current variable is a derived quantity if var not in hist_file_var_list: vres = res.get(var, {}) - #Initialiaze list for constituents - #NOTE: This is if the variable is NOT derivable but needs + # Initialiaze list for constituents + # NOTE: This is if the variable is NOT derivable but needs # an empty list as a check later constit_list = [] - #intialize boolean to check if variable is derivable + # intialize boolean to check if variable is derivable derive = False # assume it can't be derived and update if it can - #intialize boolean for regular CAM variable constituents + # intialize boolean for regular CAM variable constituents try_cam_constits = True - #Check first if variable is potentially part of a CAM-CHEM run + # Check first if variable is potentially part of a CAM-CHEM run if "derivable_from_cam_chem" in vres: constit_list = vres["derivable_from_cam_chem"] if constit_list: if all(item in hist_file_ds.data_vars for item in constit_list): - #Set check to look for regular CAM constituents in variable defaults + # Set check to look for regular CAM constituents in variable defaults try_cam_constits = False derive = True msg = f"create time series for {case_name}:" @@ -577,10 +575,10 @@ def call_ncrcat(cmd): self.debug_log(msg) else: self.debug_log(constit_errmsg) - #End if - #End if + # End if + # End if - #If not CAM-CHEM, check regular CAM runs + # If not CAM-CHEM, check regular CAM runs if try_cam_constits: if "derivable_from" in vres: derive = True @@ -595,32 +593,32 @@ def call_ncrcat(cmd): der_from_msg += "or set appropriate argument in variable " der_from_msg += "defaults yaml file." self.debug_log(der_from_msg) - #End if + # End if - #Check if this variable can be derived + # Check if this variable can be derived if (derive) and (constit_list): for constit in constit_list: if constit not in diag_var_list: diag_var_list.append(constit) - #Add variable to list to derive + # Add variable to list to derive vars_to_derive.append(var) - #Add constituent list to variable key in dictionary + # Add constituent list to variable key in dictionary constit_dict[var] = constit_list continue - #Log if this variable can be derived but is missing list of constituents + # Log if this variable can be derived but is missing list of constituents elif (derive) and (not constit_list): self.debug_log(constit_errmsg) continue - #Lastly, raise error if the variable is not a derived quanitity but is also not - #in the history file(s) + # Lastly, raise error if the variable is not a derived quanitity but is also not + # in the history file(s) else: msg = f"WARNING: {var} is not in the file {hist_files[0]} " msg += "nor can it be derived.\n" msg += "\t ** No time series will be generated." print(msg) continue - #End if - #End if (var in var_diag_list) + # End if + # End if (var in var_diag_list) # Check if variable has a "lev" dimension according to first file: has_lev = bool("lev" in hist_file_ds[var].dims) @@ -954,7 +952,7 @@ def setup_run_cvdp(self): ) # End if - #intialize objects that might not be declared later + # intialize objects that might not be declared later case_name_baseline = None baseline_ts_loc = None syears_baseline = None @@ -1099,33 +1097,33 @@ def derive_variables(self, res=None, hist_str=None, vars_to_derive=None, ts_dir= """ - #Loop through derived variables + # Loop through derived variables for var in vars_to_derive: print(f"\t - deriving time series for {var}") - #Grab list of constituents for this variable + # Grab list of constituents for this variable constit_list = constit_dict[var] - #Grab all required time series files for derived variable + # Grab all required time series files for derived variable constit_files = [] for constit in constit_list: - #Check if the constituent file is present, if so add it to list + # Check if the constituent file is present, if so add it to list if hist_str: const_glob_str = f"*{hist_str}*.{constit}.*.nc" else: const_glob_str = f"*.{constit}.*.nc" - #end if + # end if if glob.glob(os.path.join(ts_dir, const_glob_str)): constit_files.append(glob.glob(os.path.join(ts_dir, const_glob_str ))[0]) - #Check if all the necessary constituent files were found + # Check if all the necessary constituent files were found if len(constit_files) != len(constit_list): ermsg = f"\t ** Not all constituent files present; {var} cannot be calculated." ermsg += f" Please remove {var} from 'diag_var_list' or find the " ermsg += "relevant CAM files.\n" print(ermsg) if constit_files: - #Add what's missing to debug log + # Add what's missing to debug log dmsg = "create time series:" dmsg += "\n\tneeded constituents for derivation of " dmsg += f"{var}:\n\t\t- {constit_list}\n\tfound constituent file(s) in " @@ -1140,13 +1138,13 @@ def derive_variables(self, res=None, hist_str=None, vars_to_derive=None, ts_dir= self.debug_log(dmsg) else: - #Open a new dataset with all the constituent files/variables + # Open a new dataset with all the constituent files/variables ds = xr.open_mfdataset(constit_files).compute() - + # create new file name for derived variable derived_file = constit_files[0].replace(constit_list[0], var) - #Check if clobber is true for file + # Check if clobber is true for file if Path(derived_file).is_file(): if overwrite: Path(derived_file).unlink() @@ -1156,29 +1154,29 @@ def derive_variables(self, res=None, hist_str=None, vars_to_derive=None, ts_dir= print(msg) continue - #NOTE: this will need to be changed when derived equations are more complex! - JR + # NOTE: this will need to be changed when derived equations are more complex! - JR if var == "RESTOM": der_val = ds["FSNT"]-ds["FLNT"] else: - #Loop through all constituents and sum + # Loop through all constituents and sum der_val = 0 for v in constit_list: der_val += ds[v] - #Set derived variable name and add to dataset + # Set derived variable name and add to dataset der_val.name = var ds[var] = der_val - #Aerosol Calculations - #---------------------------------------------------------------------------------- - #These will be multiplied by rho (density of dry air) + # Aerosol Calculations + # ---------------------------------------------------------------------------------- + # These will be multiplied by rho (density of dry air) ds_pmid_done = False ds_t_done = False # User-defined defaults might not include aerosol zonal list azl = res.get("aerosol_zonal_list", []) if var in azl: - #Only calculate once for all aerosol vars + # Only calculate once for all aerosol vars if not ds_pmid_done: ds_pmid = _load_dataset(glob.glob(os.path.join(ts_dir, "*.PMID.*"))[0]) ds_pmid_done = True @@ -1200,16 +1198,16 @@ def derive_variables(self, res=None, hist_str=None, vars_to_derive=None, ts_dir= print(errmsg) continue - #Multiply aerosol by dry air density (rho): (P/Rd*T) + # Multiply aerosol by dry air density (rho): (P/Rd*T) ds[var] = ds[var]*(ds_pmid["PMID"]/(res["Rgas"]*ds_t["T"])) - #Sulfate conversion factor + # Sulfate conversion factor if var == "SO4": ds[var] = ds[var]*(96./115.) - #---------------------------------------------------------------------------------- + # ---------------------------------------------------------------------------------- - #Drop all constituents from final saved dataset - #These are not necessary because they have their own time series files + # Drop all constituents from final saved dataset + # These are not necessary because they have their own time series files ds_final = ds.drop_vars(constit_list) ds_final.to_netcdf(derived_file, unlimited_dims='time', mode='w') @@ -1221,7 +1219,6 @@ def setup_run_mdtf(self): """ - copy_files_only = False # True (copy files but don't run), False (copy files and run MDTF) # Note that the MDTF variable test_mode (set in the mdtf_info of the yaml file) # has a different meaning: Data is fetched but PODs are not run. @@ -1266,7 +1263,6 @@ def setup_run_mdtf(self): if mdtf_info[var] == "default": mdtf_info[var] = plot_path - # # Write the input settings json file # @@ -1435,10 +1431,9 @@ def move_tsfiles_for_mdtf(self, verbose): # end for case - ######## -#Helper Function(s) +# Helper Function(s) def _load_dataset(fils): @@ -1475,5 +1470,5 @@ def my_formatwarning(msg, *args, **kwargs): else: return xr.open_dataset(fils[0]) #End if -#End def +# End def ######## diff --git a/lib/adf_info.py b/lib/adf_info.py index d1fd5889e..3b35495f7 100644 --- a/lib/adf_info.py +++ b/lib/adf_info.py @@ -209,7 +209,8 @@ def __init__(self, config_file, debug=False): input_ts_loc = Path(input_ts_baseline) #Get years from pre-made timeseries file(s) - found_syear_baseline, found_eyear_baseline = self.get_climo_yrs_from_ts(input_ts_loc, data_name) + found_syear_baseline, found_eyear_baseline = self.get_climo_yrs_from_ts( + input_ts_loc, data_name) found_yr_range = np.arange(found_syear_baseline,found_eyear_baseline,1) #History file path isn't needed if user is running ADF directly on time series. @@ -633,7 +634,7 @@ def case_nicknames(self): @property def hist_string(self): """ Return the history string name to the user if requested.""" - return self.__hist_str + return self.get_cam_info('hist_str') ######### @@ -776,18 +777,18 @@ def get_climo_yrs_from_ts(self, input_ts_loc, case_name): #Average time dimension over time bounds, if bounds exist: if 'time_bnds' in cam_ts_data: - timeBoundsName = 'time_bnds' + time_bounds_name = 'time_bnds' elif 'time_bounds' in cam_ts_data: - timeBoundsName = 'time_bounds' + time_bounds_name = 'time_bounds' else: - timeBoundsName = None - - if timeBoundsName: + time_bounds_name = None + + if time_bounds_name: time = cam_ts_data['time'] #NOTE: force `load` here b/c if dask & time is cftime, #throws a NotImplementedError: - time = xr.DataArray(cam_ts_data[timeBoundsName].load().mean(dim='nbnd').values, + time = xr.DataArray(cam_ts_data[time_bounds_name].load().mean(dim='nbnd').values, dims=time.dims, attrs=time.attrs) cam_ts_data['time'] = time cam_ts_data.assign_coords(time=time) From febddf5e1a110e44bbb90148cc778d18db35b460 Mon Sep 17 00:00:00 2001 From: Brian Medeiros Date: Thu, 11 Jul 2024 15:43:28 -0600 Subject: [PATCH 46/46] added load method for timeseries files --- lib/adf_dataset.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/lib/adf_dataset.py b/lib/adf_dataset.py index a6c0be005..35d26a865 100644 --- a/lib/adf_dataset.py +++ b/lib/adf_dataset.py @@ -165,6 +165,7 @@ def get_timeseries_file(self, case, field): ts_files = sorted(ts_loc.glob(ts_filenames)) return ts_files + def get_ref_timeseries_file(self, field): if self.adf.compare_obs: return None @@ -175,6 +176,33 @@ def get_ref_timeseries_file(self, field): return ts_files + def load_timeseries_dataset(self, fils): + if (len(fils) == 0): + warnings.warn("Input file list is empty.") + return None + elif (len(fils) > 1): + ds = xr.open_mfdataset(fils, decode_times=False) + else: + sfil = str(fils[0]) + if not Path(sfil).is_file(): + warnings.warn(f"Expecting to find file: {sfil}") + return None + ds = xr.open_dataset(sfil, decode_times=False) + if ds is None: + warnings.warn(f"invalid data on load_dataset") + # assign time to midpoint of interval (even if it is already) + if 'time_bnds' in ds: + t = ds['time_bnds'].mean(dim='nbnd') + t.attrs = ds['time'].attrs + ds = ds.assign_coords({'time':t}) + elif 'time_bounds' in ds: + t = ds['time_bounds'].mean(dim='nbnd') + t.attrs = ds['time'].attrs + ds = ds.assign_coords({'time':t}) + else: + warnings.warn("Timeseries file does not have time bounds info.") + return xr.decode_cf(ds) + def get_ref_regrid_file(self, case, field): model_rg_loc = Path(self.adf.get_basic_info("cam_regrid_loc", required=True)) return sorted(model_rg_loc.glob(f"{case}_{field}_*.nc"))