From b928839fbdde5c3a13b01e237eb4b213fc2457d2 Mon Sep 17 00:00:00 2001 From: justin-richling Date: Tue, 19 Mar 2024 11:16:06 -0600 Subject: [PATCH 01/10] Fix incorrect variable name --- scripts/regridding/regrid_and_vert_interp.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/regridding/regrid_and_vert_interp.py b/scripts/regridding/regrid_and_vert_interp.py index 37aa5b345..a501de16d 100644 --- a/scripts/regridding/regrid_and_vert_interp.py +++ b/scripts/regridding/regrid_and_vert_interp.py @@ -532,7 +532,7 @@ def _regrid_and_interpolate_levs(model_dataset, var_name, regrid_dataset=None, r #Regrid model data to match target grid: rgdata = regrid_data(mdata, tgrid, method=1) if mdat_ofrac: - rgofrac = regrid_data(mdat_ofrac, tgrd, method=1) + rgofrac = regrid_data(mdat_ofrac, tgrid, method=1) #Regrid surface pressure if need be: if has_lev: if not regridded_ps: From 1ec9d80034ebd85698d816771e071cb10c8c21d2 Mon Sep 17 00:00:00 2001 From: justin-richling Date: Tue, 19 Mar 2024 11:17:43 -0600 Subject: [PATCH 02/10] Update adf_variable_defaults.yaml Add more aerosol info, aerosol calc list, and constants section. Also add subsection for difference colorbar options and flag for non_linear scale plotting/color scale --- lib/adf_variable_defaults.yaml | 93 ++++++++++++++++++++++++++++++++-- 1 file changed, 90 insertions(+), 3 deletions(-) diff --git a/lib/adf_variable_defaults.yaml b/lib/adf_variable_defaults.yaml index f6a765cb7..fa3664ea3 100644 --- a/lib/adf_variable_defaults.yaml +++ b/lib/adf_variable_defaults.yaml @@ -65,10 +65,14 @@ #+++++++++++++ # Available ADF Default Plot Types #+++++++++++++ - default_ptypes: ["Tables","LatLon","LatLon_Vector","Zonal","Meridional", "NHPolar","SHPolar","Special"] +#+++++++++++++ +# Constants +#+++++++++++++ +Rgas: 287.04 #[J/K/Kg]=8.314/0.028965 + #+++++++++++++ # Category: Microphysics #+++++++++++++ @@ -125,6 +129,9 @@ FICE: # Category: Aerosols #+++++++++++ +#List of zonal areosols +aerosol_zonal_list: ["BC","POM","SO4","SOA","NH4HSO4","DUST","SeaSalt"] + AODDUST: category: "Aerosols" colormap: "Oranges" @@ -185,6 +192,87 @@ SO2: SOAG: category: "Aerosols" +BC: + colormap: "RdBu_r" + diff_colormap: "BrBG" + scale_factor: 1000000000 + add_offset: 0 + new_unit: "ug/m3" + mpl: + colorbar: + label : "ug/m3" + category: "Aerosols" + derivable_from: ["bc_a1", "bc_a4"] + +POM: + colormap: "RdBu_r" + diff_colormap: "BrBG" + scale_factor: 1000000000 + add_offset: 0 + new_unit: "ug/m3" + mpl: + colorbar: + label : "ug/m3" + category: "Aerosols" + derivable_from: ["pom_a1", "pom_a4"] + +SO4: + colormap: "RdBu_r" + diff_colormap: "BrBG" + scale_factor: 1000000000 + add_offset: 0 + new_unit: "ug/m3" + mpl: + colorbar: + label : "ug/m3" + category: "Aerosols" + derivable_from: ["so4_a1", "so4_a2", "so4_a3", "so4_a5"] + +SOA: + colormap: "RdBu_r" + diff_colormap: "BrBG" + scale_factor: 1000000000 + add_offset: 0 + new_unit: "ug/m3" + mpl: + colorbar: + label : "ug/m3" + category: "Aerosols" + derivable_from: ["soa_a1", "soa_a2"] + +DUST: + colormap: "RdBu_r" + contour_levels: [0,0.1,0.25,0.4,0.6,0.8,1.4,2,3,4,8,12,30,48,114,180] + non_linear: True + diff_colormap: "BrBG" + scale_factor: 1000000000 + add_offset: 0 + new_unit: "ug/m3" + mpl: + colorbar: + label : "ug/m3" + category: "Aerosols" + derivable_from: ["dst_a1", "dst_a2", "dst_a3"] + +SeaSalt: + colormap: "RdBu_r" + contour_levels: [0,0.05,0.075,0.2,0.3,0.4,0.7,1,1.5,2,4,6,15,24,57,90] + non_linear: True + diff_colormap: "BrBG" + scale_factor: 1000000000 + add_offset: 0 + new_unit: "ug/m3" + mpl: + colorbar: + label : "ug/m3" + ticks: [0.05,0.2,0.4,1,2,6,24,90] + diff_colorbar: + label : "ug/m3" + ticks: [-10,8,6,4,2,0,-2,-4,-6,-8,-10] + category: "Aerosols" + derivable_from: ["ncl_a1", "ncl_a2", "ncl_a3"] + + #+++++++++++++++++ # Category: Budget @@ -501,7 +589,6 @@ QFLX: # Category: Surface variables #+++++++++++++++++ - PBLH: category: "Surface variables" obs_file: "PBLH_ERA5_monthly_climo_197901-202112.nc" @@ -1097,4 +1184,4 @@ utendwtem: obs_var_name: "utendwtem" #----------- -#End of File +#End of File \ No newline at end of file From 417f79553772085c1a36f84325689675c86a8811 Mon Sep 17 00:00:00 2001 From: justin-richling Date: Tue, 19 Mar 2024 11:21:25 -0600 Subject: [PATCH 03/10] Update derived variable section This moves away from the nco way of deriving variables and moves it to xarray. This is mostly so that any derived variable that has more than two constituents can be processed. The nco way only worked for two constituents. Also add check for any aerosol variables that will require rho (density of dry air) and create PMID and T time series files (if not already requested). Additionally, add section in derived variables to calculate aerosol values. --- lib/adf_diag.py | 195 ++++++++++++++++++++++++++++++++++-------------- 1 file changed, 137 insertions(+), 58 deletions(-) diff --git a/lib/adf_diag.py b/lib/adf_diag.py index 610fecb67..a5286fd0e 100644 --- a/lib/adf_diag.py +++ b/lib/adf_diag.py @@ -416,7 +416,7 @@ def call_ncrcat(cmd): emsg = f"Provided baseline 'cam_hist_loc' directory '{starting_location}' " emsg += "not found. Script is ending here." else: - emsg = "Provided 'cam_hist_loc' directory '{starting_location}' not found." + emsg = f"Provided 'cam_hist_loc' directory '{starting_location}' not found." emsg += " Script is ending here." # End if @@ -534,6 +534,17 @@ def call_ncrcat(cmd): vars_to_derive = [] # create copy of var list that can be modified for derivable variables diag_var_list = self.diag_var_list + + # Aerosol Requirements + #--------------------- + #Always make sure PMID and T are made if aerosols are desired in config file + if any(diag_var in res["aerosol_zonal_list"] for diag_var in diag_var_list): + if "PMID" not in diag_var_list: + diag_var_list += ["PMID"] + if "T" not in diag_var_list: + diag_var_list += ["T"] + #End aerosol + for var in diag_var_list: if var not in hist_file_var_list: vres = res.get(var, {}) @@ -632,9 +643,9 @@ def call_ncrcat(cmd): with mp.Pool(processes=self.num_procs) as mpool: _ = mpool.map(call_ncrcat, list_of_commands) - if vars_to_derive: + if vars_to_derive: self.derive_variables( - vars_to_derive=vars_to_derive, ts_dir=ts_dir[case_idx] + res=res, vars_to_derive=vars_to_derive, ts_dir=ts_dir[case_idx] ) # End with @@ -1008,7 +1019,7 @@ def setup_run_cvdp(self): ######### - def derive_variables(self, vars_to_derive=None, ts_dir=None, overwrite=None): + def derive_variables(self, res=None, vars_to_derive=None, ts_dir=None, overwrite=None): """ Derive variables acccording to steps given here. Since derivations will depend on the variable, each variable to derive will need its own set of steps below. @@ -1017,66 +1028,134 @@ def derive_variables(self, vars_to_derive=None, ts_dir=None, overwrite=None): If the file for the derived variable exists, the kwarg `overwrite` determines whether to overwrite the file (true) or exit with a warning message. + """ for var in vars_to_derive: - if var == "PRECT": - # PRECT can be found by simply adding PRECL and PRECC - # grab file names for the PRECL and PRECC files from the case ts directory - if glob.glob(os.path.join(ts_dir, "*PRECC*")) and glob.glob( - os.path.join(ts_dir, "*PRECL*") - ): - constit_files = sorted(glob.glob(os.path.join(ts_dir, "*PREC*"))) - else: - ermsg = "PRECC and PRECL were not both present; PRECT cannot be calculated." - ermsg += " Please remove PRECT from diag_var_list or find the relevant CAM files." - raise FileNotFoundError(ermsg) - # create new file name for PRECT - prect_file = constit_files[0].replace("PRECC", "PRECT") - if Path(prect_file).is_file(): - if overwrite: - Path(prect_file).unlink() - else: - print( - f"[{__name__}] Warning: PRECT file was found and overwrite is False. Will use existing file." - ) - continue - # append PRECC to the file containing PRECL - os.system(f"ncks -A -v PRECC {constit_files[0]} {constit_files[1]}") - # create new file with the sum of PRECC and PRECL - os.system( - f"ncap2 -s 'PRECT=(PRECC+PRECL)' {constit_files[1]} {prect_file}" - ) - if var == "RESTOM": - # RESTOM = FSNT-FLNT - # Have to be more precise than with PRECT because FSNTOA, FSTNC, etc are valid variables - if glob.glob(os.path.join(ts_dir, "*.FSNT.*")) and glob.glob( - os.path.join(ts_dir, "*.FLNT.*") - ): - input_files = [ - sorted(glob.glob(os.path.join(ts_dir, f"*.{v}.*"))) - for v in ["FLNT", "FSNT"] - ] - constit_files = [] - for elem in input_files: - constit_files += elem - else: - ermsg = "FSNT and FLNT were not both present; RESTOM cannot be calculated." - ermsg += " Please remove RESTOM from diag_var_list or find the relevant CAM files." - raise FileNotFoundError(ermsg) - # create new file name for RESTOM - derived_file = constit_files[0].replace("FLNT", "RESTOM") + print(f"\t - deriving time series for {var}") + + #Check whether there are parts to derive from and if there is an associated equation + vres = res.get(var, {}) + if "derivable_from" in vres: + constit_list = vres['derivable_from'] + else: + print("WARNING: No constituents listed in defaults config file, moving on") + pass + + constit_files = [] + for constit in constit_list: + if glob.glob(os.path.join(ts_dir, f"*.{constit}.*.nc")): + constit_files.append(glob.glob(os.path.join(ts_dir, f"*.{constit}.*"))[0]) + + #Check if all the constituent files were found + if len(constit_files) != len(constit_list): + ermsg = f"Not all constituent files present; {var} cannot be calculated." + ermsg += f" Please remove {var} from diag_var_list or find the relevant CAM files." + print(ermsg) + else: + #Open a new dataset with all the constituent files/variables + ds = xr.open_mfdataset(constit_files) + + # create new file name for derived variable + derived_file = constit_files[0].replace(constit_list[0], var) + print(derived_file,"\n") if Path(derived_file).is_file(): if overwrite: Path(derived_file).unlink() else: print( - f"[{__name__}] Warning: RESTOM file was found and overwrite is False. Will use existing file." + f"[{__name__}] Warning: '{var}' file was found and overwrite is False. Will use existing file." ) - continue - # append FSNT to the file containing FLNT - os.system(f"ncks -A -v FLNT {constit_files[0]} {constit_files[1]}") - # create new file with the difference of FLNT and FSNT - os.system( - f"ncap2 -s 'RESTOM=(FSNT-FLNT)' {constit_files[1]} {derived_file}" - ) + return None + + variables = {} + for i in constit_list: + variables[i] = ds[constit_list[0]].dims + + #Get coordinate values from the first constituent file + coords={'lat': ds[constit_list[0]].lat.values, 'lon': ds[constit_list[0]].lon.values, + "time": ds[constit_list[0]].time.values} + + #Create data arrays from each constituent + #These variables will all be added to one file that will eventually contain the + #derived variable as well + #NOTE: This has to be done via xarray dataArrays + # - There might be a way with xarray dataSets but none that have worked thus far + data_arrays = [] + for var_const, dims in variables.items(): + values = ds[var_const].values + da = xr.DataArray(values, coords=coords, dims=dims) + data_arrays.append(da) + + #NOTE: this will need to be changed when derived equations are more than sums! - JR + ds[var] = sum(data_arrays) + + #Aerosol Calculations - used for zonal plots + #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"]: + + #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 + if not ds_pmid: + errmsg = f"Missing necessary files for rho calculation.\n" + errmsg += "Please make sure 'PMID' is in the CAM run for aerosol calculations" + print(errmsg) + continue + if not ds_t_done: + ds_t = _load_dataset(glob.glob(os.path.join(ts_dir, "*.T.*"))[0]) + ds_t_done = True + if not ds_t: + errmsg = f"Missing necessary files for rho calculation.\n" + errmsg += "Please make sure 'T' is in the CAM run for aerosol calculations" + print(errmsg) + continue + + #Multiply aerosol by rho + ds[var] = ds[var]*(ds_pmid["PMID"]/(res["Rgas"]*ds_t["T"])) + #Sulfate conversion factor + if var == "SO4": + ds[var] = ds[var]*(96/115) + + ds.to_netcdf(derived_file, unlimited_dims='time', mode='w') +######## + +#Helper Function(s) +def _load_dataset(fils): + """ + This method exists to get an xarray Dataset from input file information that can be passed into the plotting methods. + + Parameters + ---------- + fils : list + strings or paths to input file(s) + + Returns + ------- + xr.Dataset + + Notes + ----- + When just one entry is provided, use `open_dataset`, otherwise `open_mfdatset` + """ + import warnings # use to warn user about missing files. + + #Format warning messages: + def my_formatwarning(msg, *args, **kwargs): + """Issue `msg` as warning.""" + return str(msg) + '\n' + warnings.formatwarning = my_formatwarning + + if len(fils) == 0: + warnings.warn(f"Input file list is empty.") + return None + elif len(fils) > 1: + return xr.open_mfdataset(fils, combine='by_coords') + else: + return xr.open_dataset(fils[0]) + #End if +#End def +######## \ No newline at end of file From 6ea5aa1455fba8ea0c21a4c3bb8f9975da983184 Mon Sep 17 00:00:00 2001 From: justin-richling Date: Tue, 19 Mar 2024 11:22:22 -0600 Subject: [PATCH 04/10] Add non linear colorbar option and difference cmap opts --- lib/plotting_functions.py | 27 ++++++++++++++++++++++----- 1 file changed, 22 insertions(+), 5 deletions(-) diff --git a/lib/plotting_functions.py b/lib/plotting_functions.py index 5e470123e..5eb4c7660 100644 --- a/lib/plotting_functions.py +++ b/lib/plotting_functions.py @@ -98,6 +98,7 @@ import geocat.comp as gcomp from mpl_toolkits.axes_grid1.inset_locator import inset_axes from matplotlib.lines import Line2D +import matplotlib.cm as cm from adf_diag import AdfDiag from adf_base import AdfError @@ -1753,6 +1754,7 @@ def prep_contour_plot(adata, bdata, diffdata, **kwargs): - 'subplots_opt': mpl kwargs for subplots - 'contourf_opt': mpl kwargs for contourf - 'colorbar_opt': mpl kwargs for colorbar + - 'diff_colorbar_opt' : mpl kwargs for difference colorbar - 'normdiff': color normalization for difference panel - 'cmapdiff': colormap for difference panel - 'levelsdiff': contour levels for difference panel @@ -1776,16 +1778,28 @@ def prep_contour_plot(adata, bdata, diffdata, **kwargs): if 'contour_levels' in kwargs: levels1 = kwargs['contour_levels'] - norm1 = mpl.colors.Normalize(vmin=min(levels1), vmax=max(levels1)) + if ('non_linear' in kwargs) and (kwargs['non_linear']): + cmap_obj = cm.get_cmap(cmap1) + norm1 = mpl.colors.BoundaryNorm(levels1, cmap_obj.N) + else: + norm1 = mpl.colors.Normalize(vmin=min(levels1), vmax=max(levels1)) elif 'contour_levels_range' in kwargs: assert len(kwargs['contour_levels_range']) == 3, \ "contour_levels_range must have exactly three entries: min, max, step" levels1 = np.arange(*kwargs['contour_levels_range']) - norm1 = mpl.colors.Normalize(vmin=min(levels1), vmax=max(levels1)) + if ('non_linear' in kwargs) and (kwargs['non_linear']): + cmap_obj = cm.get_cmap(cmap1) + norm1 = mpl.colors.BoundaryNorm(levels1, cmap_obj.N) + else: + norm1 = mpl.colors.Normalize(vmin=min(levels1), vmax=max(levels1)) else: levels1 = np.linspace(minval, maxval, 12) - norm1 = mpl.colors.Normalize(vmin=minval, vmax=maxval) + if ('non_linear' in kwargs) and (kwargs['non_linear']): + cmap_obj = cm.get_cmap(cmap1) + norm1 = mpl.colors.BoundaryNorm(levels1, cmap_obj.N) + else: + norm1 = mpl.colors.Normalize(vmin=minval, vmax=maxval) #End if #Check if the minval and maxval are actually different. If not, @@ -1840,16 +1854,19 @@ def prep_contour_plot(adata, bdata, diffdata, **kwargs): subplots_opt = {} contourf_opt = {} colorbar_opt = {} + diff_colorbar_opt = {} # extract any MPL kwargs that should be passed on: if 'mpl' in kwargs: subplots_opt.update(kwargs['mpl'].get('subplots',{})) contourf_opt.update(kwargs['mpl'].get('contourf',{})) colorbar_opt.update(kwargs['mpl'].get('colorbar',{})) + diff_colorbar_opt.update(kwargs['mpl'].get('diff_colorbar',{})) #End if return {'subplots_opt': subplots_opt, 'contourf_opt': contourf_opt, 'colorbar_opt': colorbar_opt, + 'diff_colorbar_opt': diff_colorbar_opt, 'normdiff': normdiff, 'cmapdiff': cmapdiff, 'levelsdiff': levelsdiff, @@ -1899,6 +1916,7 @@ def plot_zonal_mean_and_save(wks, case_nickname, base_nickname, shrink: 0.4 ``` """ + from matplotlib.colors import LogNorm # style the plot: # We should think about how to do plot customization and defaults. @@ -1936,7 +1954,6 @@ def plot_zonal_mean_and_save(wks, case_nickname, base_nickname, levs_diff = np.unique(np.array(cp_info['levelsdiff'])) - if len(levs) < 2: img0, ax[0] = zonal_plot(adata['lat'], azm, ax=ax[0]) ax[0].text(0.4, 0.4, empty_message, transform=ax[0].transAxes, bbox=props) @@ -1954,7 +1971,7 @@ def plot_zonal_mean_and_save(wks, case_nickname, base_nickname, ax[2].text(0.4, 0.4, empty_message, transform=ax[2].transAxes, bbox=props) else: img2, ax[2] = zonal_plot(adata['lat'], diff, ax=ax[2], norm=cp_info['normdiff'],cmap=cp_info['cmapdiff'],levels=cp_info['levelsdiff'],**cp_info['contourf_opt']) - fig.colorbar(img2, ax=ax[2], location='right',**cp_info['colorbar_opt']) + fig.colorbar(img2, ax=ax[2], location='right',**cp_info['diff_colorbar_opt']) ax[0].set_title(case_title, loc='left', fontsize=tiFontSize) ax[1].set_title(base_title, loc='left', fontsize=tiFontSize) From a56ebf7d91bd701eb584a196b8cc92af43f9e50e Mon Sep 17 00:00:00 2001 From: justin-richling Date: Tue, 19 Mar 2024 11:33:03 -0600 Subject: [PATCH 05/10] Update adf_diag.py --- lib/adf_diag.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/adf_diag.py b/lib/adf_diag.py index a5286fd0e..908e1b8f3 100644 --- a/lib/adf_diag.py +++ b/lib/adf_diag.py @@ -643,7 +643,7 @@ def call_ncrcat(cmd): with mp.Pool(processes=self.num_procs) as mpool: _ = mpool.map(call_ncrcat, list_of_commands) - if vars_to_derive: + if vars_to_derive: self.derive_variables( res=res, vars_to_derive=vars_to_derive, ts_dir=ts_dir[case_idx] ) From 82c1317bbcb42ff87b9b05c29f0969fc56a97465 Mon Sep 17 00:00:00 2001 From: justin-richling Date: Tue, 19 Mar 2024 11:38:48 -0600 Subject: [PATCH 06/10] Fix github errors --- lib/adf_diag.py | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/lib/adf_diag.py b/lib/adf_diag.py index 908e1b8f3..a149934b2 100644 --- a/lib/adf_diag.py +++ b/lib/adf_diag.py @@ -1040,7 +1040,7 @@ def derive_variables(self, res=None, vars_to_derive=None, ts_dir=None, overwrite constit_list = vres['derivable_from'] else: print("WARNING: No constituents listed in defaults config file, moving on") - pass + continue constit_files = [] for constit in constit_list: @@ -1055,7 +1055,7 @@ def derive_variables(self, res=None, vars_to_derive=None, ts_dir=None, overwrite else: #Open a new dataset with all the constituent files/variables ds = xr.open_mfdataset(constit_files) - + # create new file name for derived variable derived_file = constit_files[0].replace(constit_list[0], var) print(derived_file,"\n") @@ -1067,26 +1067,27 @@ def derive_variables(self, res=None, vars_to_derive=None, ts_dir=None, overwrite f"[{__name__}] Warning: '{var}' file was found and overwrite is False. Will use existing file." ) return None - + variables = {} for i in constit_list: variables[i] = ds[constit_list[0]].dims - + #Get coordinate values from the first constituent file - coords={'lat': ds[constit_list[0]].lat.values, 'lon': ds[constit_list[0]].lon.values, - "time": ds[constit_list[0]].time.values} + coords={'lat': ds[constit_list[0]].lat.values, + 'lon': ds[constit_list[0]].lon.values, + "time": ds[constit_list[0]].time.values} #Create data arrays from each constituent #These variables will all be added to one file that will eventually contain the #derived variable as well #NOTE: This has to be done via xarray dataArrays - # - There might be a way with xarray dataSets but none that have worked thus far + # - There might be a way with xarray dataSets but none that have worked thus far data_arrays = [] for var_const, dims in variables.items(): values = ds[var_const].values da = xr.DataArray(values, coords=coords, dims=dims) data_arrays.append(da) - + #NOTE: this will need to be changed when derived equations are more than sums! - JR ds[var] = sum(data_arrays) @@ -1095,13 +1096,13 @@ def derive_variables(self, res=None, vars_to_derive=None, ts_dir=None, overwrite ds_pmid_done = False ds_t_done = False if var in res["aerosol_zonal_list"]: - + #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 if not ds_pmid: - errmsg = f"Missing necessary files for rho calculation.\n" + errmsg = "Missing necessary files for rho calculation.\n" errmsg += "Please make sure 'PMID' is in the CAM run for aerosol calculations" print(errmsg) continue @@ -1109,7 +1110,7 @@ def derive_variables(self, res=None, vars_to_derive=None, ts_dir=None, overwrite ds_t = _load_dataset(glob.glob(os.path.join(ts_dir, "*.T.*"))[0]) ds_t_done = True if not ds_t: - errmsg = f"Missing necessary files for rho calculation.\n" + errmsg = "Missing necessary files for rho calculation.\n" errmsg += "Please make sure 'T' is in the CAM run for aerosol calculations" print(errmsg) continue @@ -1119,7 +1120,7 @@ def derive_variables(self, res=None, vars_to_derive=None, ts_dir=None, overwrite #Sulfate conversion factor if var == "SO4": ds[var] = ds[var]*(96/115) - + ds.to_netcdf(derived_file, unlimited_dims='time', mode='w') ######## @@ -1150,7 +1151,7 @@ def my_formatwarning(msg, *args, **kwargs): warnings.formatwarning = my_formatwarning if len(fils) == 0: - warnings.warn(f"Input file list is empty.") + warnings.warn("Input file list is empty.") return None elif len(fils) > 1: return xr.open_mfdataset(fils, combine='by_coords') From 8ea38385fca88e04156b0268ea8a79d83535c246 Mon Sep 17 00:00:00 2001 From: justin-richling Date: Tue, 19 Mar 2024 11:39:48 -0600 Subject: [PATCH 07/10] Update plotting_functions.py --- lib/plotting_functions.py | 1 - 1 file changed, 1 deletion(-) diff --git a/lib/plotting_functions.py b/lib/plotting_functions.py index 5eb4c7660..ae4a61049 100644 --- a/lib/plotting_functions.py +++ b/lib/plotting_functions.py @@ -1916,7 +1916,6 @@ def plot_zonal_mean_and_save(wks, case_nickname, base_nickname, shrink: 0.4 ``` """ - from matplotlib.colors import LogNorm # style the plot: # We should think about how to do plot customization and defaults. From c896ad0d088e5c9d3fc492263b2ad1114c3ebd84 Mon Sep 17 00:00:00 2001 From: justin-richling Date: Tue, 16 Apr 2024 14:03:44 -0600 Subject: [PATCH 08/10] Change greek symbols for units --- lib/adf_variable_defaults.yaml | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/lib/adf_variable_defaults.yaml b/lib/adf_variable_defaults.yaml index fa3664ea3..3eca5028c 100644 --- a/lib/adf_variable_defaults.yaml +++ b/lib/adf_variable_defaults.yaml @@ -71,6 +71,8 @@ default_ptypes: ["Tables","LatLon","LatLon_Vector","Zonal","Meridional", #+++++++++++++ # Constants #+++++++++++++ + +#Dry Air Gas Constant: Rgas: 287.04 #[J/K/Kg]=8.314/0.028965 #+++++++++++++ @@ -197,10 +199,10 @@ BC: diff_colormap: "BrBG" scale_factor: 1000000000 add_offset: 0 - new_unit: "ug/m3" + new_unit: '$\mu$g/m3' mpl: colorbar: - label : "ug/m3" + label : '$\mu$g/m3' category: "Aerosols" derivable_from: ["bc_a1", "bc_a4"] @@ -209,10 +211,10 @@ POM: diff_colormap: "BrBG" scale_factor: 1000000000 add_offset: 0 - new_unit: "ug/m3" + new_unit: '$\mu$g/m3' mpl: colorbar: - label : "ug/m3" + label : '$\mu$g/m3' category: "Aerosols" derivable_from: ["pom_a1", "pom_a4"] @@ -221,10 +223,10 @@ SO4: diff_colormap: "BrBG" scale_factor: 1000000000 add_offset: 0 - new_unit: "ug/m3" + new_unit: '$\mu$g/m3' mpl: colorbar: - label : "ug/m3" + label : '$\mu$g/m3' category: "Aerosols" derivable_from: ["so4_a1", "so4_a2", "so4_a3", "so4_a5"] @@ -233,10 +235,10 @@ SOA: diff_colormap: "BrBG" scale_factor: 1000000000 add_offset: 0 - new_unit: "ug/m3" + new_unit: '$\mu$g/m3' mpl: colorbar: - label : "ug/m3" + label : '$\mu$g/m3' category: "Aerosols" derivable_from: ["soa_a1", "soa_a2"] @@ -247,10 +249,10 @@ DUST: diff_colormap: "BrBG" scale_factor: 1000000000 add_offset: 0 - new_unit: "ug/m3" + new_unit: '$\mu$g/m3' mpl: colorbar: - label : "ug/m3" + label : '$\mu$g/m3' category: "Aerosols" derivable_from: ["dst_a1", "dst_a2", "dst_a3"] @@ -261,13 +263,13 @@ SeaSalt: diff_colormap: "BrBG" scale_factor: 1000000000 add_offset: 0 - new_unit: "ug/m3" + new_unit: '$\mu$g/m3' mpl: colorbar: - label : "ug/m3" + label : '$\mu$g/m3' ticks: [0.05,0.2,0.4,1,2,6,24,90] diff_colorbar: - label : "ug/m3" + label : '$\mu$g/m3' ticks: [-10,8,6,4,2,0,-2,-4,-6,-8,-10] category: "Aerosols" derivable_from: ["ncl_a1", "ncl_a2", "ncl_a3"] From e440ecbef8270e93af3a035c80010b3dd49028d8 Mon Sep 17 00:00:00 2001 From: justin-richling Date: Tue, 16 Apr 2024 14:08:15 -0600 Subject: [PATCH 09/10] Update adf_diag.py Clean up derived variable function to grab coordinates and dimensions for new variable netcdf file. Also, remove constituents from derived variable file since they will have their own. --- lib/adf_diag.py | 76 ++++++++++++++++++++++++------------------------- 1 file changed, 37 insertions(+), 39 deletions(-) diff --git a/lib/adf_diag.py b/lib/adf_diag.py index a149934b2..9c57c1f36 100644 --- a/lib/adf_diag.py +++ b/lib/adf_diag.py @@ -535,15 +535,16 @@ def call_ncrcat(cmd): # create copy of var list that can be modified for derivable variables diag_var_list = self.diag_var_list - # Aerosol Requirements - #--------------------- - #Always make sure PMID and T are made if aerosols are desired in config file - if any(diag_var in res["aerosol_zonal_list"] for diag_var in diag_var_list): - if "PMID" not in diag_var_list: + # Aerosol Calcs + #-------------- + #Always make sure PMID is made if aerosols are desired in config file + if "PMID" not in diag_var_list: + if any(item in res["aerosol_zonal_list"] for item in diag_var_list): diag_var_list += ["PMID"] - if "T" not in diag_var_list: + if "T" not in diag_var_list: + if any(item in res["aerosol_zonal_list"] for item in diag_var_list): diag_var_list += ["T"] - #End aerosol + #End aerosol calcs for var in diag_var_list: if var not in hist_file_var_list: @@ -1042,6 +1043,7 @@ def derive_variables(self, res=None, vars_to_derive=None, ts_dir=None, overwrite print("WARNING: No constituents listed in defaults config file, moving on") continue + #Grab all required time series files for derived var constit_files = [] for constit in constit_list: if glob.glob(os.path.join(ts_dir, f"*.{constit}.*.nc")): @@ -1055,10 +1057,11 @@ def derive_variables(self, res=None, vars_to_derive=None, ts_dir=None, overwrite else: #Open a new dataset with all the constituent files/variables ds = xr.open_mfdataset(constit_files) - + # create new file name for derived variable derived_file = constit_files[0].replace(constit_list[0], var) - print(derived_file,"\n") + + #Check if clobber is true for file if Path(derived_file).is_file(): if overwrite: Path(derived_file).unlink() @@ -1066,43 +1069,33 @@ def derive_variables(self, res=None, vars_to_derive=None, ts_dir=None, overwrite print( f"[{__name__}] Warning: '{var}' file was found and overwrite is False. Will use existing file." ) - return None - - variables = {} - for i in constit_list: - variables[i] = ds[constit_list[0]].dims - - #Get coordinate values from the first constituent file - coords={'lat': ds[constit_list[0]].lat.values, - 'lon': ds[constit_list[0]].lon.values, - "time": ds[constit_list[0]].time.values} - - #Create data arrays from each constituent - #These variables will all be added to one file that will eventually contain the - #derived variable as well - #NOTE: This has to be done via xarray dataArrays - # - There might be a way with xarray dataSets but none that have worked thus far - data_arrays = [] - for var_const, dims in variables.items(): - values = ds[var_const].values - da = xr.DataArray(values, coords=coords, dims=dims) - data_arrays.append(da) - - #NOTE: this will need to be changed when derived equations are more than sums! - JR - ds[var] = sum(data_arrays) + continue + + #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 + der_val = 0 + for v in constit_list: + der_val += ds[v] + + #Set derived variable name and add to dataset + der_val.name = var + ds[var] = der_val #Aerosol Calculations - used for zonal plots #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"]: - + #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 if not ds_pmid: - errmsg = "Missing necessary files for rho calculation.\n" + errmsg = f"Missing necessary files for dry air density (rho) calculation.\n" errmsg += "Please make sure 'PMID' is in the CAM run for aerosol calculations" print(errmsg) continue @@ -1110,18 +1103,23 @@ def derive_variables(self, res=None, vars_to_derive=None, ts_dir=None, overwrite ds_t = _load_dataset(glob.glob(os.path.join(ts_dir, "*.T.*"))[0]) ds_t_done = True if not ds_t: - errmsg = "Missing necessary files for rho calculation.\n" + errmsg = f"Missing necessary files for dry air density (rho) calculation.\n" errmsg += "Please make sure 'T' is in the CAM run for aerosol calculations" print(errmsg) continue - #Multiply aerosol by rho + #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 if var == "SO4": - ds[var] = ds[var]*(96/115) + 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 + ds_final = ds.drop_vars(constit_list) + ds_final.to_netcdf(derived_file, unlimited_dims='time', mode='w') - ds.to_netcdf(derived_file, unlimited_dims='time', mode='w') ######## #Helper Function(s) From 5626dd5d7740bdfecc55adad638b794bbb40a2fc Mon Sep 17 00:00:00 2001 From: justin-richling Date: Tue, 16 Apr 2024 14:36:15 -0600 Subject: [PATCH 10/10] Add statement about derived variables in docstring --- lib/adf_variable_defaults.yaml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/lib/adf_variable_defaults.yaml b/lib/adf_variable_defaults.yaml index 3eca5028c..b8ecaa1ee 100644 --- a/lib/adf_variable_defaults.yaml +++ b/lib/adf_variable_defaults.yaml @@ -54,6 +54,9 @@ # derivable_from -> If not present in the available output files, the variable can be derived from # other variables that are present (e.g. PRECT can be derived from PRECC and PRECL), # which are specified in this list +# NOTE: this is not very flexible at the moment! It can only handle variables that +# are sums of the constituents. Futher flexibility is being explored. +# # # Final Note: Please do not modify this file unless you plan to push your changes back to the ADF repo. # If you would like to modify this file for your personal ADF runs then it is recommended