From 3fbe5d00777213b6ac7a6b2e1fffe96297bca042 Mon Sep 17 00:00:00 2001 From: Shawn Honomichl Date: Thu, 18 Jan 2024 11:25:16 -0700 Subject: [PATCH 1/5] V0 --- config_cam_baseline_example.yaml | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/config_cam_baseline_example.yaml b/config_cam_baseline_example.yaml index 339de56fb..9a569c502 100644 --- a/config_cam_baseline_example.yaml +++ b/config_cam_baseline_example.yaml @@ -58,7 +58,7 @@ # Note that the string 'USER-NAME-NOT-SET' is used in the jupyter script # to check for a failure to customize # -user: 'USER-NAME-NOT-SET' +user: 'shawnh' #This first set of variables specify basic info used by all diagnostic runs: @@ -87,7 +87,7 @@ diag_basic_info: obs_data_loc: /glade/work/nusbaume/SE_projects/model_diagnostics/ADF_obs #Location where re-gridded and interpolated CAM climatology files are stored: - cam_regrid_loc: /glade/scratch/${user}/ADF/regrid + cam_regrid_loc: /glade/derecho/scratch/${user}/ADF/regrid #Overwrite CAM re-gridded files? #If false, or missing, then regridding will be skipped for regridded variables @@ -95,7 +95,7 @@ diag_basic_info: cam_overwrite_regrid: false #Location where diagnostic plots are stored: - cam_diag_plot_loc: /glade/scratch/${user}/ADF/plots + cam_diag_plot_loc: /glade/derecho/scratch/${user}/ADF/plots #Location of ADF variable plotting defaults YAML file: #If left blank or missing, ADF/lib/adf_variable_defaults.yaml will be used @@ -152,7 +152,7 @@ diag_cam_climo: cam_hist_loc: /glade/p/cesm/ADF/${diag_cam_climo.cam_case_name} #Location of CAM climatologies (to be created and then used by this script) - cam_climo_loc: /glade/scratch/${user}/ADF/${diag_cam_climo.cam_case_name}/climo + cam_climo_loc: /glade/derecho/scratch/${user}/ADF/${diag_cam_climo.cam_case_name}/climo #model year when time series files should start: #Note: Leaving this entry blank will make time series @@ -180,7 +180,7 @@ diag_cam_climo: cam_overwrite_ts: false #Location where time series files are (or will be) stored: - cam_ts_loc: /glade/scratch/${user}/ADF/${diag_cam_climo.cam_case_name}/ts + cam_ts_loc: /glade/derecho/scratch/${user}/ADF/${diag_cam_climo.cam_case_name}/ts #---------------------- @@ -274,7 +274,7 @@ diag_cam_baseline_climo: cam_hist_loc: /glade/p/cesm/ADF/${diag_cam_baseline_climo.cam_case_name} #Location of baseline CAM climatologies: - cam_climo_loc: /glade/scratch/${user}/ADF/${diag_cam_baseline_climo.cam_case_name}/climo + cam_climo_loc: /glade/derecho/scratch/${user}/ADF/${diag_cam_baseline_climo.cam_case_name}/climo #model year when time series files should start: #Note: Leaving this entry blank will make time series @@ -301,7 +301,7 @@ diag_cam_baseline_climo: cam_overwrite_ts: false #Location where time series files are (or will be) stored: - cam_ts_loc: /glade/scratch/${user}/ADF/${diag_cam_baseline_climo.cam_case_name}/ts + cam_ts_loc: /glade/derecho/scratch/${user}/ADF/${diag_cam_baseline_climo.cam_case_name}/ts #This fourth set of variables provides settings for calling the Climate Variability # Diagnostics Package (CVDP). If cvdp_run is set to true the CVDP will be set up and From 85eadad2f90a084b69a80ecdc58e7626090e5e81 Mon Sep 17 00:00:00 2001 From: shawnusaf <84995386+shawnusaf@users.noreply.github.com> Date: Thu, 22 Feb 2024 12:51:41 -0500 Subject: [PATCH 2/5] O3 DIAG V0 Adds in V0 of the Ozone diagnostics package --- config_amwg_default_plots.yaml | 2 + config_cam_baseline_example.yaml | 2 + scripts/plotting/ozone_diagnostics.py | 875 ++++++++++++++++++++++++++ 3 files changed, 879 insertions(+) create mode 100755 scripts/plotting/ozone_diagnostics.py diff --git a/config_amwg_default_plots.yaml b/config_amwg_default_plots.yaml index b432bc7e9..51f5101be 100644 --- a/config_amwg_default_plots.yaml +++ b/config_amwg_default_plots.yaml @@ -297,6 +297,7 @@ plotting_scripts: - zonal_mean - polar_map - cam_taylor_diagram + - ozone_diagnostics #- tape_recorder #- tem #To plot TEM, please un-comment fill-out #the "tem_info" section below @@ -346,6 +347,7 @@ diag_var_list: - ICEFRAC - OCNFRAC - LANDFRAC + - O3 # diff --git a/config_cam_baseline_example.yaml b/config_cam_baseline_example.yaml index 9a569c502..391b61f57 100644 --- a/config_cam_baseline_example.yaml +++ b/config_cam_baseline_example.yaml @@ -363,6 +363,7 @@ plotting_scripts: - polar_map - cam_taylor_diagram - qbo + - ozone_diagnostics #- tape_recorder #- tem #To plot TEM, please un-comment fill-out #the "tem_info" section below @@ -388,6 +389,7 @@ diag_var_list: - FSNT - FLNT - LANDFRAC + - O3 # diff --git a/scripts/plotting/ozone_diagnostics.py b/scripts/plotting/ozone_diagnostics.py new file mode 100755 index 000000000..75200bba7 --- /dev/null +++ b/scripts/plotting/ozone_diagnostics.py @@ -0,0 +1,875 @@ +""" +Generates an ozone-specific diagnostics based on the previous ncl-based +diagnostics package + +Created 19 January, 2024 +Shawn Honomichl +Associate Scientist III +NCAR/ACOM +shawnh@ucar.edu + +Functions +--------- +ozone_diagnostics(adfobj) + use ADF object to make maps + +define_regions(InputRegion) + +define_stations(InputStation) + +open_process_sonde_data_simone(obsdir) + +Subplot_O3(ax,X0,Y0,X1,Model0,Model1,SondeMean,SondeLow,SondeHigh,PLevel,plottype,case_base,case_test,ylim,Quadrant,Compare_Obs) + +Plot_Seasonal_Cycle_Profile(X0,X01,Y0,X1,X11,Y1,X2,X22,Y2,X3,X33,Y3,SondeMean0,SondeMean1,SondeMean2,SondeMean3,SondeLow0,SondeLow1,SondeLow2,SondeLow3,SondeHigh0,SondeHigh1,SondeHigh2,SondeHigh3,Model0I,Model01I,Model1I,Model11I,Model2I,Model21I,Model3I,Model31I,ylim0,ylim1,ylim2,ylim3,PLevel0,PLevel1,PLevel2,PLevel3,case_base,case_test,Station,oFile,plottype,Compare_Obs) + +get_model_data(ClimoFile) + +process_model_seasonal_cycle(MinLon,MaxLon,MinLat,MaxLat,Model_Dat,pnew,intyp,kxtrp,Station_Lons,Station_Lats) + +process_model_profiles(Model_Dat,O3_0,PS_0,pnew,intyp,kxtrp,ILAT,ILON,lat_0,lon_0) + + 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 + +""" + +#Import modules: +import numpy as np +import xarray as xr +import warnings # use to warn user about missing files. +import os +import matplotlib.pyplot as plt +import matplotlib as mpl +import Ngl +from scipy.interpolate import RegularGridInterpolator +from scipy import stats +from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) +import os.path + +def my_formatwarning(msg, *args, **kwargs): + # ignore everything except the message + return str(msg) + '\n' + +warnings.formatwarning = my_formatwarning + +#----------------------------------------------------------------------------------------- +#Lookup pertinent region info +#----------------------------------------------------------------------------------------- + +#Region_Number : ['region short name', 'region long name',region_min_lat, region_max_lat,region_min_lon,region_max_lon] +def define_regions(InputRegion): + + Region_Lookup={0 : ['not_assigned','not assigned',0,0,0,0], + 1 : ['nh_polar_west','NH Polar West',70.0,90.0,-135.0+360,-45.0+360], + 2 : ['nh_polar_east','NH Polar East',58.0,90.0,-45.0,45.0], + 3 : ['canada','Canada',48.0,62.0,-135.0+360,-45.0+360], + 4 : ['west_europe','Western Europe',43.0,57.5,0.0,25.0], + 5 : ['eastern_us','Eastern US',34.0,40.0,-95.0+360,-75.0+360], + 6 : ['japan','Japan',30.0,40.0,120.0,150.0], + 7 : ['nh_tropic','NH Sub Tropics',15.0,30.0,90.0,225.0], + 8 : ['tropics1','W-Pacific/E-Indian Ocean',-20.0,0.0,90.0,225.0], + 9 : ['tropics2','Equatorial Americas',-15.0,15.0,225.0,315.0], + 10 : ['tropics3','Atlantic/Africa',-15.0,15.0,-45.0,45.0], + 11 : ['sh_midlat','SH Midlatitudes',-57.5,-40.0,135.0,180.0], + 12 : ['sh_polar','SH Polar',-90.0,-58.0,0,360], + 13 : ['Boulder','Boulder',37.0,42.0,-110.0+360.0,-100.0+360.0], + } + + #-------------------------------------------------------------------------------------- + #try to locate the region + #-------------------------------------------------------------------------------------- + try: + RegionInfo=Region_Lookup[InputRegion] + except: + print("Input region not defined") #if region not in the lookup then auto-assign to 0 + RegionInfo=Region_Lookup[0] + + #-------------------------------------------------------------------------------------- + #return the region info and the total number of regions in the table. + #-------------------------------------------------------------------------------------- + return RegionInfo,len(Region_Lookup)-1 + +#----------------------------------------------------------------------------------------- +#Lookup pertinent ozonesonde station info +#----------------------------------------------------------------------------------------- +def define_stations(InputStation): + + #Station : ['country',lat N,lon E] + Station_Lookup={0 : ['Alert','Canada',82.5,-62.33], + 1 : ['Eureka','Canada',80.05,-86.42], + 2 : ['Ny_Alesund','Norway',78.92,11.93], + 3 : ['Resolute','Canada',74.7,-95.0], + 4 : ['Scoresbysund','Greenland',70.48,-21.97], + 5 : ['Lerwick','UK',60.13,-1.1+360.0], + 6 : ['Churchill','Canada',58.77,-94.17], + 7 : ['Edmonton','Canada',53.53,-113.5], + 8 : ['Goose_Bay','Canada',53.29,-60.39], + 9 : ['Legionowo','Poland',52.4,20.97], + 10 : ['Lindenberg','Germany',52.21,14.12], + 11 : ['DeBilt','Netherlands',52.10,5.18], + 12 : ['Uccle','Belgium',50.8,4.35], + 13 : ['Praha','Czech Republic',50.01,14.45], + 14 : ['Hohenpeissenberg','Germany',47.80,11.02], + 15 : ['Payerne','Switzerland',46.82,6.95], + 16 : ['Madrid','Spain',40.45,-3.72], + 17 : ['Boulder','Colorado',39.99,-105.26], + 18 : ['Wallops_Island','Maryland',37.94,-75.46], + 19 : ['Trinidad Head','California',41.05,-124.15], + 20 : ['Huntsville','Alabama',34.73,-86.64], + 21 : ['Sapporo','Japan',43.06,141.35], + 22 : ['Tateno','Japan',36.05,140.13], + 23 : ['Kagoshima','Japan',31.6,130.6], + 24 : ['Naha','Japan',26.21,127.68], + 25 : ['Hongkong','China',22.32,114.17], + 26 : ['Paramaribo','Suriname',5.75,55.2], + 27 : ['Hilo','USA',19.72,155.07], + 28 : ['San Cristobal','Galapagos',-0.87,-89.44], + 29 : ['Nairobi','Kenya',-1.29,36.82], + 30 : ['Natal','Brazil',-5.83,-35.2], + 31 : ['Ascension','Atlantic',-7.95,-14.36], + 32 : ['Watukosek','Indonesia',-7.57,112.66], + 33 : ['Samoa','Pacific',-14.25,-170.56], + 34 : ['Fiji','Fiji',-17.71,178.07], + 35 : ['Reunion','Africa',-21.1,55.4], + 36 : ['Broadmeadows','Australia',-37.68,144.92], + 37 : ['Lauder','New Zealand',-45.04,169.68], + 38 : ['Macquarie','Australia',-54.50,158.95], + 39 : ['Marambio','Antarctica',-64.0,-56.0], + 40 : ['Neumayer','Antarctica',-70.62,8.37], + 41 : ['Syowa','Antarctica',-69.01,39.59]} + + + #-------------------------------------------------------------------------------------- + #try to locate the station + #-------------------------------------------------------------------------------------- + try : + StationInfo=Station_Lookup[InputStation] + except: + print("Input station not defined") + + return StationInfo,len(Station_Lookup)-1 + +#----------------------------------------------------------------------------------------- +#Process ozone sonde climatology +#----------------------------------------------------------------------------------------- +def open_process_sonde_data_simone(obsdir): + + #Grab and package all of the data that will be used for plotting + NRegions=define_regions(0)[1] + + O3_MeanC=[] + O3_WidthC=[] + O3_StdDevC=[] + Region=[] + PressureC=[] + + for i in range(1,NRegions+1): + Region_Info=define_regions(i)[0] + ifileID=Region_Info[0] + + Region.append(ifileID) + + Check_File=obsdir+'ozonesondes_'+ifileID+'1995_2011.nc' + + if (os.path.isfile(Check_File)): + print("Using Ozone Climatology File: ",Check_File) + O3_Mean = xr.open_dataset(Check_File).o3_mean + O3_Width = xr.open_dataset(Check_File).o3_width + O3_StdDev = xr.open_dataset(Check_File).o3_std + if (ifileID == 'Boulder'): + Pressure = xr.open_dataset(Check_File).press #hPa + else: + Pressure = xr.open_dataset(Check_File).levels #hPa + + if (i == 1): + O3_MeanC=O3_Mean + O3_WidthC=O3_Width + O3_StdDevC=O3_StdDev + PressureC=Pressure + else: + O3_MeanC = np.dstack( (O3_MeanC,O3_Mean)) + O3_WidthC = np.dstack( (O3_WidthC,O3_Width)) + O3_StdDevC = np.dstack( (O3_StdDevC,O3_Width)) + PressureC=np.dstack( (PressureC,Pressure)) + else: + print("Error opening ozonesonde file ",Check_File) + return -1 + + return O3_MeanC,O3_WidthC,O3_StdDevC,np.squeeze(PressureC),Region + +#----------------------------------------------------------------------------------------- +#Subroutine - plots the ozone data subplots +#----------------------------------------------------------------------------------------- +def Subplot_O3(ax,X0,Y0,X1,Model0,Model1,SondeMean,SondeLow,SondeHigh,PLevel,plottype,case_base,case_test,ylim,Quadrant,Compare_Obs): + + #plottype == 0 ; Seasonal Cycle + #plottype == 1 ; Profiles + + #Blue is the base case - if comparing obs will not be plotted + #Red is the test case - should always be plotted + + Model_linewidth=0.75 + + if plottype == 0: + ax.plot(Y0,X1,'tab:red',linestyle='dashed',linewidth=Model_linewidth) + ax.plot(Y0,Model1,'tab:red',linestyle='solid',linewidth=Model_linewidth) + if Compare_Obs <= 0: + ax.plot(Y0,X0,'tab:blue',linestyle='dashed',linewidth=Model_linewidth) + ax.plot(Y0,Model0,'tab:blue',linestyle='solid',linewidth=Model_linewidth) + + ax.set_title(PLevel,loc='left') + + props = dict(boxstyle='round',facecolor='wheat',alpha=0.5) + if Compare_Obs <= 0: + ax.text(0.05,0.95," Mean: \nAbs. Diff:\n r:",transform=ax.transAxes,fontsize=5,verticalalignment='top',bbox=props,color='black') + else: + ax.text(0.05,0.95," Mean: \nAbs. Diff:\n r:",transform=ax.transAxes,fontsize=5,verticalalignment='top',bbox=props,color='black') + + meanS0=np.mean(SondeMean) + meanS0S='%.1f' % meanS0 + + meanX01=np.mean( [X1,Model1]) + meanX01S='%.1f' % meanX01 + AbsDif1=np.absolute(meanX01-meanS0) + AbsDif1S='%.1f' % AbsDif1 + PCorr01=stats.pearsonr(X1,SondeMean) + PCorr01S='%.2f' % PCorr01[0] + + if Compare_Obs <= 0: + ax.text(0.70,0.92,meanX01S,fontsize=5,color='tab:red',transform=ax.transAxes) + ax.text(0.70,0.86,AbsDif1S,fontsize=5,color='tab:red',transform=ax.transAxes) + ax.text(0.70,0.80,PCorr01S,fontsize=5,color='tab:red',transform=ax.transAxes) + else: + ax.text(0.50,0.92,meanX01S,fontsize=5,color='tab:red',transform=ax.transAxes) + ax.text(0.50,0.86,AbsDif1S,fontsize=5,color='tab:red',transform=ax.transAxes) + ax.text(0.50,0.80,PCorr01S,fontsize=5,color='tab:red',transform=ax.transAxes) + + ax.text(0.30,0.92,meanS0S,fontsize=5,color='black',transform=ax.transAxes) + + if Compare_Obs <= 0: + + meanX0=np.mean([X0,Model0]) + meanX0S='%.1f' % meanX0 + AbsDif0=np.absolute(meanX0-meanS0) + AbsDif0S='%.1f' % AbsDif0 + PCorr0=stats.pearsonr(X0,SondeMean) + PCorr0S='%.2f' % PCorr0[0] + + ax.text(0.50,0.92,meanX0S,fontsize=5,color='tab:blue',transform=ax.transAxes) + ax.text(0.50,0.86,AbsDif0S,fontsize=5,color='tab:blue',transform=ax.transAxes) + ax.text(0.50,0.80,PCorr0S,fontsize=5,color='tab:blue',transform=ax.transAxes) + + #legend (bottom left) side of bottom left plot + if Quadrant == 2: + props1 = dict(boxstyle='square',facecolor='white',alpha=0.5) + #ax.text(0.03,0.03," \n \n \n \n ",transform=ax.transAxes,fontsize=5.25,verticalalignment='bottom',bbox=props1,color='black') + + if Compare_Obs >0: + ax.text(0.12,0.15,"Region Avg.",fontsize=5.0,color='tab:red',transform=ax.transAxes) + ax.text(0.12,0.09,case_test,fontsize=5.0,color='tab:red',transform=ax.transAxes) + ax.plot([0.03,0.11],[0.16,0.16],'tab:red',linestyle='dashed',transform=ax.transAxes,linewidth=0.75) + ax.plot([0.03,0.11],[0.10,0.10],'tab:red',linestyle='solid',transform=ax.transAxes,linewidth=0.75) + else: + ax.text(0.12,0.27,"Region Avg.",fontsize=5.0,color='tab:blue',transform=ax.transAxes) + ax.text(0.12,0.21,case_base,fontsize=5.0,color='tab:blue',transform=ax.transAxes) + ax.plot([0.03,0.11],[0.28,0.28],'tab:blue',linestyle='dashed',transform=ax.transAxes,linewidth=0.75) + ax.plot([0.03,0.11],[0.22,0.22],'tab:blue',linestyle='solid',transform=ax.transAxes,linewidth=0.75) + ax.text(0.12,0.15,"Region Avg.",fontsize=5.0,color='tab:red',transform=ax.transAxes) + ax.text(0.12,0.09,case_test,fontsize=5.0,color='tab:red',transform=ax.transAxes) + ax.plot([0.03,0.11],[0.16,0.16],'tab:red',linestyle='dashed',transform=ax.transAxes,linewidth=0.75) + ax.plot([0.03,0.11],[0.10,0.10],'tab:red',linestyle='solid',transform=ax.transAxes,linewidth=0.75) + + ax.plot([0.06],[0.04],'ko',markersize=1.25,transform=ax.transAxes) + ax.text(0.12,0.03,"Ozonesondes",fontsize=5.0,color='black',transform=ax.transAxes) + + #plot the ozonesonde data for a type 0 plot + ax.plot(Y0,SondeMean,'ko',markersize=1.25) + for i in range(1,13): + ax.plot([i,i],[SondeLow[i-1],SondeHigh[i-1]],color='black',linestyle='solid',linewidth=1) + + #set the appropriate x/y axis information + ax.set(xlabel="Month", ylabel="Ozone (ppb)", xlim=(1,12),ylim=ylim,xscale='linear',yscale='linear') + ax.xaxis.set_major_locator(MultipleLocator(2)) + + else: + + ax.plot(X1,Y0,'tab:red',linestyle='dashed',linewidth=Model_linewidth) + ax.plot(Model1,Y0,'tab:red',linestyle='solid',linewidth=Model_linewidth) + ax.set_title(PLevel,loc='left') + + if Compare_Obs <= 0: + ax.plot(X0,Y0,'tab:blue',linestyle='dashed',linewidth=Model_linewidth) + ax.plot(Model0,Y0,'tab:blue',linestyle='solid',linewidth=Model_linewidth) + + if Quadrant == 2: + offset=0.67 + props1 = dict(boxstyle='square',facecolor='white',alpha=0.5) + #ax.text(0.03,0.03+offset," \n \n \n \n ",transform=ax.transAxes,fontsize=5.25,verticalalignment='bottom',bbox=props1,color='black') + + if Compare_Obs > 0: + ax.text(0.12,0.16+offset,"Ozonesondes",fontsize=5.0,color='black',transform=ax.transAxes) + ax.plot([0.06],[0.17+offset],'ko',markersize=1.25,transform=ax.transAxes) + else: + ax.text(0.12,0.03+offset,"Ozonesondes",fontsize=5.0,color='black',transform=ax.transAxes) + ax.plot([0.06],[0.04+offset],'ko',markersize=1.25,transform=ax.transAxes) + + ax.text(0.12,0.27+offset,"Region Avg.",fontsize=5.0,color='tab:red',transform=ax.transAxes) + ax.text(0.12,0.21+offset,case_test,fontsize=5.0,color='tab:red',transform=ax.transAxes) + ax.plot([0.03,0.11],[0.28+offset,0.28+offset],'tab:red',linestyle='dashed',transform=ax.transAxes,linewidth=0.75) + ax.plot([0.03,0.11],[0.22+offset,0.22+offset],'tab:red',linestyle='solid',transform=ax.transAxes,linewidth=0.75) + + + if Compare_Obs <= 0: + ax.text(0.12,0.15+offset,"Region Avg.",fontsize=5.0,color='tab:blue',transform=ax.transAxes) + ax.text(0.12,0.09+offset,case_base,fontsize=5.0,color='tab:blue',transform=ax.transAxes) + ax.plot([0.03,0.11],[0.16+offset,0.16+offset],'tab:blue',linestyle='dashed',transform=ax.transAxes,linewidth=0.75) + ax.plot([0.03,0.11],[0.10+offset,0.10+offset],'tab:blue',linestyle='solid',transform=ax.transAxes,linewidth=0.75) + + #plot the ozonesonde data for a type 1 plot + Step=1 + SondeMean=SondeMean[1::Step] + Y00=Y0[1::Step] + SondeLow=SondeLow[1::Step] + SondeHigh=SondeHigh[1::Step] + ax.plot(SondeMean,Y00,'ko',markersize=1.25) + for i in range(0,len(Y00)): + ax.plot([SondeLow[i],SondeHigh[i]],[Y00[i],Y00[i]],color='black',linestyle='solid',linewidth=1) + + #set the appropriate x/y axis information + ax.set(xlabel="Ozone (ppb)", ylabel="Pressure (hPa)", xlim=ylim,ylim=[1000.0,0.0],xscale='linear',yscale='linear') + + #Set tick parameters + ax.xaxis.set_minor_locator(AutoMinorLocator()) + ax.yaxis.set_minor_locator(AutoMinorLocator()) + ax.tick_params(top=True, bottom=True, left=True, right=True) + + return True + +#----------------------------------------------------------------------------------------- +#Subroutine - plot the seasonal cycle +#----------------------------------------------------------------------------------------- +def Plot_Seasonal_Cycle_Profile(X0,X01,Y0,X1,X11,Y1,X2,X22,Y2,X3,X33,Y3,SondeMean0,SondeMean1,SondeMean2,SondeMean3,SondeLow0,SondeLow1,SondeLow2,SondeLow3,SondeHigh0,SondeHigh1,SondeHigh2,SondeHigh3,Model0I,Model01I,Model1I,Model11I,Model2I,Model21I,Model3I,Model31I,ylim0,ylim1,ylim2,ylim3,PLevel0,PLevel1,PLevel2,PLevel3,case_base,case_test,Station,oFile,plottype,Compare_Obs): + + #-------------------------------------------------------------------------------------- + #set the default font size for the plots + #-------------------------------------------------------------------------------------- + mpl.rcParams['font.size'] = 8 + plt.rcParams['figure.figsize'] = [6.4,4.8] + + #-------------------------------------------------------------------------------------- + #set the subplots + #-------------------------------------------------------------------------------------- + fig, axs = plt.subplots(nrows=2,ncols=2,sharex=False,sharey=False) + + #-------------------------------------------------------------------------------------- + #Primary subplots + #-------------------------------------------------------------------------------------- + PTL = Subplot_O3(axs[0,0],X0,Y0,X01,Model0I,Model01I,SondeMean0,SondeLow0,SondeHigh0,PLevel0,plottype,case_base,case_test,ylim0,0,Compare_Obs)#top left quadrant 0 + PTR = Subplot_O3(axs[0,1],X1,Y1,X11,Model1I,Model11I,SondeMean1,SondeLow1,SondeHigh1,PLevel1,plottype,case_base,case_test,ylim1,1,Compare_Obs)#top right quadrant 1 + PBL = Subplot_O3(axs[1,0],X2,Y2,X22,Model2I,Model21I,SondeMean2,SondeLow2,SondeHigh2,PLevel2,plottype,case_base,case_test,ylim2,2,Compare_Obs)#bottom left quadrant 2 + PBR = Subplot_O3(axs[1,1],X3,Y3,X33,Model3I,Model31I,SondeMean3,SondeLow3,SondeHigh3,PLevel3,plottype,case_base,case_test,ylim3,3,Compare_Obs)#bottom right quadrant 3 + + #-------------------------------------------------------------------------------------- + #A couple of other plotting adjustments + #-------------------------------------------------------------------------------------- + fig.suptitle(Station) + plt.subplots_adjust(wspace=0.4,hspace=0.4,left=0.3) #Adjusts the spacing between the plots + plt.rcParams['font.size'] = '7' + + Output_IMG = oFile #set the output filename + + #-------------------------------------------------------------------------------------- + #Try to save the plot as an image file + #-------------------------------------------------------------------------------------- + try: + plt.savefig(Output_IMG, bbox_inches = 'tight', dpi = 200, format = 'png') #there will be more to add here... + plt.clf() + plt.close(fig) + except: + print("ERROR: Could not save file ",Output_FName,flush=True) + + return True + +#----------------------------------------------------------------------------------------- +#Subroutine - retrieve model data from a netCDF file +#----------------------------------------------------------------------------------------- +def get_model_data(ClimoFile): + + class Model_Data: + #Get the Reference Case data + if (os.path.isfile(ClimoFile)): + print("ADF Climatology File Located: ",ClimoFile) + lon = xr.open_dataset(ClimoFile).lon + lat = xr.open_dataset(ClimoFile).lat + #lev = xr.open_dataset(ClimoFile).lev + hyam = xr.open_dataset(ClimoFile).hyam + hybm = xr.open_dataset(ClimoFile).hybm + o3=xr.open_dataset(ClimoFile).O3 + ps=xr.open_dataset(ClimoFile).PS + #time=xr.open_dataset(ClimoFile).time + + #these typically won't change with time so only select the time 0 + hyam=hyam[0,:] + hybm=hybm[0,:] + + MDAT=Model_Data() + return MDAT + +#----------------------------------------------------------------------------------------- +#Subroutine - Process Model Seasonal Cycle +#----------------------------------------------------------------------------------------- +def process_model_seasonal_cycle(MinLon,MaxLon,MinLat,MaxLat,Model_Dat,pnew,intyp,kxtrp,Station_Lons,Station_Lats): + + class model_dat_proc: + + #get the model data from the base and test cases for the region + if (MinLon < 0 and MaxLon > 0): #if the region crosses the date line - do different processing + O3_00=Model_Dat.o3.sel(lon=slice(MinLon+360.0,360.0),lat=slice(MinLat,MaxLat)) + O3_01=Model_Dat.o3.sel(lon=slice(0,MaxLon),lat=slice(MinLat,MaxLat)) + O3_0 = np.concatenate( (O3_00,O3_01),axis=3) + PS_00=Model_Dat.ps.sel(lon=slice(MinLon+360.0,360.0),lat=slice(MinLat,MaxLat)) + PS_01=Model_Dat.ps.sel(lon=slice(0,MaxLon),lat=slice(MinLat,MaxLat)) + PS_0 = np.concatenate( (PS_00,PS_01),axis=2) + lon_00=Model_Dat.lon.sel(lon=slice(MinLon+360.0,360.0)) + lon_01=Model_Dat.lon.sel(lon=slice(0,MaxLon)) + lon_0 = np.concatenate( (lon_00,lon_01)) + + #resort the arrays as needed + lon_sort=lon_0.argsort() + O3_0 = O3_0[:,:,:,lon_sort] + PS_0 = PS_0[:,:,lon_sort] + lon_0 = lon_0[lon_sort] + + O3_sfc=np.squeeze(O3_0[:,-1,:,:])*1.0e9 #get the lowest model surface level data + + else: #if the region does not cross the date line + + O3_0=Model_Dat.o3.sel(lon=slice(MinLon,MaxLon),lat=slice(MinLat,MaxLat)) + PS_0=Model_Dat.ps.sel(lon=slice(MinLon,MaxLon),lat=slice(MinLat,MaxLat)) + lon_0=Model_Dat.lon.sel(lon=slice(MinLon,MaxLon)) + O3_sfc=np.squeeze(O3_0.values[:,-1,:,:])*1.0e9 + + + lat_0=Model_Dat.lat.sel(lat=slice(MinLat,MaxLat)) + + O3_sfc_04=np.mean(O3_sfc,axis=(1,2)) + + O3_0I = Ngl.vinth2p(O3_0,Model_Dat.hyam,Model_Dat.hybm,pnew,PS_0,intyp,1000.0,1,kxtrp)*1.0e9 + + #Get the seasonal cycle of the base case at each needed pressure level + #and average over the region. + O3_01=np.mean(O3_0I[:,0,:,:],axis=(1,2)) + O3_02=np.mean(O3_0I[:,1,:,:],axis=(1,2)) + O3_03=np.mean(O3_0I[:,2,:,:],axis=(1,2)) + O3_04=np.mean(O3_0I[:,3,:,:],axis=(1,2)) + + #Interpolate the model at each ozonesonde location and extract the needed information for the plot + Station_Find = np.where( (Station_Lons >= MinLon) & (Station_Lons <= MaxLon) & (Station_Lats >= MinLat) & (Station_Lats <= MaxLat)) + ILAT=Station_Lats[Station_Find] + ILON=Station_Lons[Station_Find] + + #Set the points to interpolate to + for i in range(1,13): + for j in range(0,len(ILAT)): + if (i == 1 and j == 0): + O3_Pt=[i,float(ILAT[j]),float(ILON[j])] + else: + O3_Pt=np.vstack( (O3_Pt,[i,float(ILAT[j]),float(ILON[j])] )) + + months=[1,2,3,4,5,6,7,8,9,10,11,12] + + #set up the regular grid interpolator for each case and level + interp_0 = RegularGridInterpolator((months,lat_0,lon_0), np.squeeze(O3_0I[:,0,:,:])) + interp_1 = RegularGridInterpolator((months,lat_0,lon_0), np.squeeze(O3_0I[:,1,:,:])) + interp_2 = RegularGridInterpolator((months,lat_0,lon_0), np.squeeze(O3_0I[:,2,:,:])) + interp_3 = RegularGridInterpolator((months,lat_0,lon_0), O3_sfc) + + #interpolate the model data at each case and pressure level + O3_station_0 = interp_0(O3_Pt) + O3_station_1 = interp_1(O3_Pt) + O3_station_2 = interp_2(O3_Pt) + O3_station_3 = interp_3(O3_Pt) + + #loop through each station interpolation and average like months + O3M=np.squeeze(O3_Pt[:,0]) + + O3_station_00=[] + O3_station_10=[] + O3_station_20=[] + O3_station_30=[] + + for i in range(1,13): + Elms = np.where(O3M == i) + O3_station_00.append(np.mean(O3_station_0[Elms])) + O3_station_10.append(np.mean(O3_station_1[Elms])) + O3_station_20.append(np.mean(O3_station_2[Elms])) + O3_station_30.append(np.mean(O3_station_3[Elms])) + + #For consistency convert the O3_station variables into numpy arrays + O3_station_00=np.array(O3_station_00) + O3_station_10=np.array(O3_station_10) + O3_station_20=np.array(O3_station_20) + O3_station_30=np.array(O3_station_30) + + Processed_Model_Dat=model_dat_proc() + + return Processed_Model_Dat + +#----------------------------------------------------------------------------------------- +#Subroutine - Process Model Profiles +#----------------------------------------------------------------------------------------- +def process_model_profiles(Model_Dat,O3_0,PS_0,pnew,intyp,kxtrp,ILAT,ILON,lat_0,lon_0): + + class model_dat_proc: + + O3_0I1 = Ngl.vinth2p(O3_0,Model_Dat.hyam,Model_Dat.hybm,pnew,PS_0,intyp,1000.0,1,kxtrp)*1.0e9 + + Locate_Bad=np.where(O3_0I1 > 10000.0) + if len(Locate_Bad) > 0: + O3_0I1[Locate_Bad]=np.nan + + #Get the monthly of the case at each needed pressure level + #and average over the region. + O3_011=np.nanmean(O3_0I1[0,:,:,:],axis=(1,2)) + O3_021=np.nanmean(O3_0I1[3,:,:,:],axis=(1,2)) + O3_031=np.nanmean(O3_0I1[6,:,:,:],axis=(1,2)) + O3_041=np.nanmean(O3_0I1[9,:,:,:],axis=(1,2)) + + #Set the points to interpolate to + for i in range(0,len(pnew)): + for j in range(0,len(ILAT)): + if (i == 0 and j == 0): + O3_Pt=[pnew[i],float(ILAT[j]),float(ILON[j])] + else: + O3_Pt=np.vstack( (O3_Pt,[pnew[i],float(ILAT[j]),float(ILON[j])] )) + + #set up the regular grid interpolator for each case and level + interp_01 = RegularGridInterpolator((pnew,lat_0,lon_0), np.squeeze(O3_0I1[0,:,:,:])) + interp_11 = RegularGridInterpolator((pnew,lat_0,lon_0), np.squeeze(O3_0I1[3,:,:,:])) + interp_21 = RegularGridInterpolator((pnew,lat_0,lon_0), np.squeeze(O3_0I1[6,:,:,:])) + interp_31 = RegularGridInterpolator((pnew,lat_0,lon_0), np.squeeze(O3_0I1[9,:,:,:])) + + #interpolate the model data at each case and pressure level + O3_station_01 = interp_01(O3_Pt) + O3_station_11 = interp_11(O3_Pt) + O3_station_21 = interp_21(O3_Pt) + O3_station_31 = interp_31(O3_Pt) + + #loop through each station interpolation and average like plevs + O3M=np.squeeze(O3_Pt[:,0]) + + O3_station_001=[] + O3_station_101=[] + O3_station_201=[] + O3_station_301=[] + + for i in range(0,len(pnew)): + Elms = np.where(O3M == pnew[i]) + O3_station_001.append(np.nanmean(O3_station_01[Elms])) + O3_station_101.append(np.nanmean(O3_station_11[Elms])) + O3_station_201.append(np.nanmean(O3_station_21[Elms])) + O3_station_301.append(np.nanmean(O3_station_31[Elms])) + + #For consistency convert the O3_station variables into numpy arrays + O3_station_001=np.array(O3_station_001) + O3_station_101=np.array(O3_station_101) + O3_station_201=np.array(O3_station_201) + O3_station_301=np.array(O3_station_301) + + Processed_Model_Dat=model_dat_proc() + + return Processed_Model_Dat + +#----------------------------------------------------------------------------------------- +#Primary Ozone Diagnostics Routine +#----------------------------------------------------------------------------------------- +def ozone_diagnostics (adfobj): + + #---------------------------------------------------------------------------------- + #Extract relevant info from the ADF + #---------------------------------------------------------------------------------- + obsdir='/glade/campaign/acom/acom-climate/tilmes/obs_data/amwg/amwg_data/obs_data/cam-chem/' #location of the ozonesonde data files + cam_climo_loc = adfobj.get_cam_info('cam_climo_loc',required=True)[0] + baseline_climo_loc = adfobj.get_baseline_info('cam_climo_loc',required=True) + ncases = len(cam_climo_loc) + + #check the number of cases. For the O3 diagnostics there needs to be 2. + #If there is more than 2, then only the first two will be used. + case_test = adfobj.get_cam_info('cam_case_name', required=True)[0] + case_base = adfobj.get_baseline_info('cam_case_name',required=True) + + #Grab all case nickname(s) + test_nicknames = adfobj.case_nicknames["test_nicknames"] + base_nickname = adfobj.case_nicknames["base_nickname"] + + plot_locations = adfobj.plot_location[0] + + climo_base=baseline_climo_loc+'/'+case_base+"_O3_climo.nc" + if os.path.isfile(climo_base): + print(climo_base," located") + else: + print(climo_base," could not be located. Exiting O3 diagnostics.") + return + + climo_test=cam_climo_loc+'/'+case_test+"_O3_climo.nc" + if os.path.isfile(climo_test): + print(climo_test," located") + else: + print(climo_test," could not be located. Exiting O3 diagnostics.") + return + + #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}") + + Compare_Obs=0 #initially assumes that user is comparing two models + if adfobj.compare_obs: + Compare_Obs=1 #if comparing a model with observations then don't need to include the base case in the plots + + print("-------------------------------") + print("Processing Ozone Diagnostics...") + print("-------------------------------") + + #-------------------------------------------------------------------------------------- + #Check and make sure that if the O3 Diagnostics are selected in the .yaml file that + #ozone (O3) is listed in the diag_var_list. + #-------------------------------------------------------------------------------------- + if not ('O3' in adfobj.diag_var_list): + msg = "No ozone ('O3') variable present" + msg += " in 'diag_var_list', so O3 diagnostic plots will" + msg += " be skipped." + print(msg) + return + + #-------------------------------------------------------------------------------------- + #Get the ozone sonde data from all of the stations + #-------------------------------------------------------------------------------------- + Data = open_process_sonde_data_simone(obsdir) + Obs_Mean=Data[0] #mean values + Obs_Width=Data[1] #width of the data + Obs_StdDev=Data[2] #Standard deviation of the data + Obs_Pressure=Data[3] #Pressure of the obs data + Regions=Data[4] #regions + + #-------------------------------------------------------------------------------------- + #Get model data from get_model_data() + #-------------------------------------------------------------------------------------- + Test_Data=get_model_data(climo_test) #Grab test model data + if Compare_Obs <= 0: + Base_Data=get_model_data(climo_base) #Grab base model data + + #-------------------------------------------------------------------------------------- + #Loop through each station and get the station coordinates + #-------------------------------------------------------------------------------------- + Station_Lats=[] + Station_Lons=[] + Station_SNames=[] + + NStations=define_stations(0)[1] + for i in range(0,NStations+1): + Station_Info=define_stations(i)[0] + if (Station_Info[3] < 0.0): + Station_Lons=np.append(Station_Lons,float(Station_Info[3])+360.0) + else: + Station_Lons=np.append(Station_Lons,float(Station_Info[3])) + Station_Lats=np.append(Station_Lats,float(Station_Info[2])) + Station_SNames=np.append(Station_SNames,Station_Info[0]) + + #-------------------------------------------------------------------------------------- + #set the parameters for interpolation to pressure levels + #-------------------------------------------------------------------------------------- + pnew = [50,250,500,900] + intyp = 2 # 1=linear, 2=log, 3=log-log + kxtrp = False # True=extrapolate (when the output pressure level is outside of the range of psrf) + months=[1,2,3,4,5,6,7,8,9,10,11,12] + + #-------------------------------------------------------------------------------------- + #loop through each region and plot the Data + #-------------------------------------------------------------------------------------- + for i in range(1,len(Regions)+1): + + #----------------------------------------------------------------------------------- + #Grab the info for the current region of interest + #----------------------------------------------------------------------------------- + Region_Info=define_regions(i)[0] + SName=Region_Info[0] + LName=Region_Info[1] + MinLat=Region_Info[2] + MaxLat=Region_Info[3] + MinLon=Region_Info[4] + MaxLon=Region_Info[5] + oFile_Seasonal = plot_locations+'/O3SeasonalCycle_'+SName+'_Special.png' + oFile_Profile = plot_locations+'/O3Profile_'+SName+'_Special.png' + + #----------------------------------------------------------------------------------- + #Check if redo_plot set and if not and plots exist already then + #add them to the website (if enabled) and then exit the routine + #----------------------------------------------------------------------------------- + if (not(redo_plot)) and (os.path.isfile(oFile_Seasonal)) and (os.path.isfile(oFile_Profile)): + print(SName,' region plots exist and redo_plot is false. Adding to website and Skipping plot.') + adfobj.add_website_data(oFile_Seasonal,SName.replace("_","")+"_SeasonalCycle", None, season="ANN",multi_case=True,category="O3_DIAGNOSTICS") + adfobj.add_website_data(oFile_Profile,SName.replace("_","")+"_Profile", None, season="ANN", multi_case=True,category="O3_DIAGNOSTICS") + continue + else: + print("Plotting Region ",LName) + #----------------------------------------------------------------------------------- + #Process the model data + #----------------------------------------------------------------------------------- + Processed_Seasonal_Cycle_Test_Data = process_model_seasonal_cycle(MinLon,MaxLon,MinLat,MaxLat,Test_Data,pnew,intyp,kxtrp,Station_Lons,Station_Lats) + if Compare_Obs <= 0: + Processed_Seasonal_Cycle_Base_Data = process_model_seasonal_cycle(MinLon,MaxLon,MinLat,MaxLat,Base_Data,pnew,intyp,kxtrp,Station_Lons,Station_Lats) + else: + Processed_Seasonal_Cycle_Base_Data=Processed_Seasonal_Cycle_Test_Data #allows the plotting routine to run without having to call it different ways. + + #----------------------------------------------------------------------------------- + #Filter to the appropriate ozone sonde regional data. + #----------------------------------------------------------------------------------- + Obs_Mean_0=Obs_Mean[:,:,i-1] + Obs_Width_0=Obs_Width[:,:,i-1] + Obs_StdDev_0=Obs_StdDev[:,:,i-1] + Obs_Pressure_0=Obs_Pressure[:,i-1] + Regions_0=Regions[i-1] + + #----------------------------------------------------------------------------------- + #Get the pressure level index for each of the 4 levels to + #filter the ozone sonde data to + #----------------------------------------------------------------------------------- + Locate_P0=np.where(Obs_Pressure_0 == pnew[0])[0] + Locate_P1=np.where(Obs_Pressure_0 == pnew[1])[0] + Locate_P2=np.where(Obs_Pressure_0 == pnew[2])[0] + Locate_P3=np.where(Obs_Pressure_0 == np.max(Obs_Pressure_0))[0] + + #----------------------------------------------------------------------------------- + #For Boulder specifically, the Index should be 2. + #----------------------------------------------------------------------------------- + if LName == 'Boulder': + Locate_P3=2 + + #----------------------------------------------------------------------------------- + #Trim the ozonesonde data to the appropriate pressure levels + #----------------------------------------------------------------------------------- + Obs_Mean_00=np.squeeze(Obs_Mean_0[Locate_P0,:]) + Obs_Mean_10=np.squeeze(Obs_Mean_0[Locate_P1,:]) + Obs_Mean_20=np.squeeze(Obs_Mean_0[Locate_P2,:]) + Obs_Mean_30=np.squeeze(Obs_Mean_0[Locate_P3,:]) + + Obs_Width_00=np.squeeze(Obs_Width_0[Locate_P0,:]) + Obs_Width_10=np.squeeze(Obs_Width_0[Locate_P1,:]) + Obs_Width_20=np.squeeze(Obs_Width_0[Locate_P2,:]) + Obs_Width_30=np.squeeze(Obs_Width_0[Locate_P3,:]) + + Obs_StdDev_00=np.squeeze(Obs_StdDev_0[Locate_P0,:]) + Obs_StdDev_10=np.squeeze(Obs_StdDev_0[Locate_P1,:]) + Obs_StdDev_20=np.squeeze(Obs_StdDev_0[Locate_P2,:]) + Obs_StdDev_30=np.squeeze(Obs_StdDev_0[Locate_P3,:]) + + Obs_Low_00=Obs_Mean_00-Obs_StdDev_00 + Obs_Low_10=Obs_Mean_10-Obs_StdDev_10 + Obs_Low_20=Obs_Mean_20-Obs_StdDev_20 + Obs_Low_30=Obs_Mean_30-Obs_StdDev_30 + + Obs_High_00=Obs_Mean_00+Obs_StdDev_00 + Obs_High_10=Obs_Mean_10+Obs_StdDev_10 + Obs_High_20=Obs_Mean_20+Obs_StdDev_20 + Obs_High_30=Obs_Mean_30+Obs_StdDev_30 + + #----------------------------------------------------------------------------------- + #Assign the output seasonal cycle filename and Plot the Seasonal Cycle Plot + #----------------------------------------------------------------------------------- + Plot_Seasonal_Cycle_Profile(Processed_Seasonal_Cycle_Base_Data.O3_01,Processed_Seasonal_Cycle_Test_Data.O3_01,months,Processed_Seasonal_Cycle_Base_Data.O3_02,Processed_Seasonal_Cycle_Test_Data.O3_02,months,Processed_Seasonal_Cycle_Base_Data.O3_03,Processed_Seasonal_Cycle_Test_Data.O3_03,months,Processed_Seasonal_Cycle_Base_Data.O3_sfc_04,Processed_Seasonal_Cycle_Test_Data.O3_sfc_04,months,Obs_Mean_00,Obs_Mean_10,Obs_Mean_20,Obs_Mean_30,Obs_Low_00,Obs_Low_10,Obs_Low_20,Obs_Low_30,Obs_High_00,Obs_High_10,Obs_High_20,Obs_High_30,Processed_Seasonal_Cycle_Base_Data.O3_station_00,Processed_Seasonal_Cycle_Test_Data.O3_station_00,Processed_Seasonal_Cycle_Base_Data.O3_station_10,Processed_Seasonal_Cycle_Test_Data.O3_station_10,Processed_Seasonal_Cycle_Base_Data.O3_station_20,Processed_Seasonal_Cycle_Test_Data.O3_station_20,Processed_Seasonal_Cycle_Base_Data.O3_station_30,Processed_Seasonal_Cycle_Test_Data.O3_station_30,[0,6000],[0,750],[0,120],[0,75],'50 hPa','250 hPa','500 hPa','Sfc',base_nickname,test_nicknames,LName,oFile_Seasonal,0,Compare_Obs) + + #----------------------------------------------------------------------------------- + #Set up for the monthly profile plots - I am calculating my own pressure level data here to ensure a tight match with the ozone sonde climatology + #----------------------------------------------------------------------------------- + pnew1=[0.01,0.05,0.1,0.15,0.25,0.5,0.75,1.0,5.0,10.0,20.0,30.0,40.0,50.0,60.0,70.0,80.0,90.0,100.0,125.0,150.0,175.0,200.0,225.0,250.0,275.0,300.0,350.0,400.0,450.0,500.0,550.0,600.0,650.0,700.0,750.0,800.0,850.0,900.0,950.0,1000.0] + pnew1=np.flip(Obs_Pressure_0) + + #----------------------------------------------------------------------------------- + #Process the model data + #----------------------------------------------------------------------------------- + Processed_Profiles_Test = process_model_profiles(Test_Data,Processed_Seasonal_Cycle_Test_Data.O3_0,Processed_Seasonal_Cycle_Test_Data.PS_0,pnew1,intyp,kxtrp,Processed_Seasonal_Cycle_Test_Data.ILAT,Processed_Seasonal_Cycle_Test_Data.ILON,Processed_Seasonal_Cycle_Test_Data.lat_0,Processed_Seasonal_Cycle_Test_Data.lon_0) + if Compare_Obs<=0: + Processed_Profiles_Base = process_model_profiles(Base_Data,Processed_Seasonal_Cycle_Base_Data.O3_0,Processed_Seasonal_Cycle_Base_Data.PS_0,pnew1,intyp,kxtrp,Processed_Seasonal_Cycle_Base_Data.ILAT,Processed_Seasonal_Cycle_Base_Data.ILON,Processed_Seasonal_Cycle_Base_Data.lat_0,Processed_Seasonal_Cycle_Base_Data.lon_0) + else: + Processed_Profiles_Base=Processed_Profiles_Test + + #----------------------------------------------------------------------------------- + #Trim the ozonesonde data to the appropriate pressure levels + #----------------------------------------------------------------------------------- + Obs_Mean_001=np.squeeze(Obs_Mean_0[:,0]) + Obs_Mean_101=np.squeeze(Obs_Mean_0[:,3]) + Obs_Mean_201=np.squeeze(Obs_Mean_0[:,6]) + Obs_Mean_301=np.squeeze(Obs_Mean_0[:,9]) + + Obs_Width_001=np.squeeze(Obs_Width_0[:,0]) + Obs_Width_101=np.squeeze(Obs_Width_0[:,3]) + Obs_Width_201=np.squeeze(Obs_Width_0[:,6]) + Obs_Width_301=np.squeeze(Obs_Width_0[:,9]) + + Obs_StdDev_001=np.squeeze(Obs_StdDev_0[:,0]) + Obs_StdDev_101=np.squeeze(Obs_StdDev_0[:,3]) + Obs_StdDev_201=np.squeeze(Obs_StdDev_0[:,6]) + Obs_StdDev_301=np.squeeze(Obs_StdDev_0[:,9]) + + Obs_Low_001=Obs_Mean_001-Obs_StdDev_001 + Obs_Low_101=Obs_Mean_101-Obs_StdDev_101 + Obs_Low_201=Obs_Mean_201-Obs_StdDev_201 + Obs_Low_301=Obs_Mean_301-Obs_StdDev_301 + + Obs_High_001=Obs_Mean_001+Obs_StdDev_001 + Obs_High_101=Obs_Mean_101+Obs_StdDev_101 + Obs_High_201=Obs_Mean_201+Obs_StdDev_201 + Obs_High_301=Obs_Mean_301+Obs_StdDev_301 + + #----------------------------------------------------------------------------------- + #Assign the output monthly profile filename and Plot the Seasonal Cycle Plot + #----------------------------------------------------------------------------------- + Plot_Seasonal_Cycle_Profile(Processed_Profiles_Base.O3_011,Processed_Profiles_Test.O3_011,pnew1,Processed_Profiles_Base.O3_021,Processed_Profiles_Test.O3_021,pnew1,Processed_Profiles_Base.O3_031,Processed_Profiles_Test.O3_031,pnew1,Processed_Profiles_Base.O3_041,Processed_Profiles_Test.O3_041,pnew1,np.flip(Obs_Mean_001),np.flip(Obs_Mean_101),np.flip(Obs_Mean_201),np.flip(Obs_Mean_301),np.flip(Obs_Low_001),np.flip(Obs_Low_101),np.flip(Obs_Low_201),np.flip(Obs_Low_301),np.flip(Obs_High_001),np.flip(Obs_High_101),np.flip(Obs_High_201),np.flip(Obs_High_301),Processed_Profiles_Base.O3_station_001,Processed_Profiles_Test.O3_station_001,Processed_Profiles_Base.O3_station_101,Processed_Profiles_Test.O3_station_101,Processed_Profiles_Base.O3_station_201,Processed_Profiles_Test.O3_station_201,Processed_Profiles_Base.O3_station_301,Processed_Profiles_Test.O3_station_301,[0,125],[0,125],[0,125],[0,125],'Jan', 'Apr','Jul','Oct',base_nickname,test_nicknames,LName,oFile_Profile,1,Compare_Obs) + + #----------------------------------------------------------------------------------- + #Once the plots have successfully run, add the web page entries (if enabled). + #----------------------------------------------------------------------------------- + adfobj.add_website_data(oFile_Seasonal,SName.replace("_","")+"_SeasonalCycle", None, season="ANN",multi_case=True,category="O3_DIAGNOSTICS") + adfobj.add_website_data(oFile_Profile,SName.replace("_","")+"_Profile", None, season="ANN", multi_case=True,category="O3_DIAGNOSTICS") + + print("Ozone Diagnostics Generated Successfully!") From df801b734a4511f0be129948a4e8470e58155ce8 Mon Sep 17 00:00:00 2001 From: shawnusaf <84995386+shawnusaf@users.noreply.github.com> Date: Fri, 11 Oct 2024 10:09:36 -0400 Subject: [PATCH 3/5] Add files via upload --- scripts/plotting/ozone_diagnostics.py | 543 ++++++++++++++++++++++---- 1 file changed, 470 insertions(+), 73 deletions(-) diff --git a/scripts/plotting/ozone_diagnostics.py b/scripts/plotting/ozone_diagnostics.py index 75200bba7..5b40590b8 100755 --- a/scripts/plotting/ozone_diagnostics.py +++ b/scripts/plotting/ozone_diagnostics.py @@ -1,6 +1,7 @@ """ Generates an ozone-specific diagnostics based on the previous ncl-based -diagnostics package +diagnostics package. Current package plots ozone sesasonal cycle, profiles, and +taylor diagrams. Created 19 January, 2024 Shawn Honomichl @@ -19,6 +20,8 @@ open_process_sonde_data_simone(obsdir) +open_process_sonde_data_simone_v2(obsdir) + Subplot_O3(ax,X0,Y0,X1,Model0,Model1,SondeMean,SondeLow,SondeHigh,PLevel,plottype,case_base,case_test,ylim,Quadrant,Compare_Obs) Plot_Seasonal_Cycle_Profile(X0,X01,Y0,X1,X11,Y1,X2,X22,Y2,X3,X33,Y3,SondeMean0,SondeMean1,SondeMean2,SondeMean3,SondeLow0,SondeLow1,SondeLow2,SondeLow3,SondeHigh0,SondeHigh1,SondeHigh2,SondeHigh3,Model0I,Model01I,Model1I,Model11I,Model2I,Model21I,Model3I,Model31I,ylim0,ylim1,ylim2,ylim3,PLevel0,PLevel1,PLevel2,PLevel3,case_base,case_test,Station,oFile,plottype,Compare_Obs) @@ -28,6 +31,14 @@ process_model_seasonal_cycle(MinLon,MaxLon,MinLat,MaxLat,Model_Dat,pnew,intyp,kxtrp,Station_Lons,Station_Lats) process_model_profiles(Model_Dat,O3_0,PS_0,pnew,intyp,kxtrp,ILAT,ILON,lat_0,lon_0) + +taylor_plot_setup(ax,baseline) + +taylor_stats_single(casedata, refdata) + +plot_taylor_data(wks, df,k, **kwargs) + +taylor_plot_finalize(wks, casenames, casecolors, syear_cases, eyear_cases, SubZones, needs_bias_labels=False,needs_case_labels=False,needs_subzone_labels=False) Parameters ---------- @@ -80,11 +91,16 @@ import os import matplotlib.pyplot as plt import matplotlib as mpl -import Ngl +import geocat.comp as gcomp from scipy.interpolate import RegularGridInterpolator from scipy import stats from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) import os.path +import pandas as pd +from matplotlib.lines import Line2D +from matplotlib.legend_handler import HandlerTuple +from scipy.stats import pearsonr +import sys def my_formatwarning(msg, *args, **kwargs): # ignore everything except the message @@ -97,22 +113,24 @@ def my_formatwarning(msg, *args, **kwargs): #----------------------------------------------------------------------------------------- #Region_Number : ['region short name', 'region long name',region_min_lat, region_max_lat,region_min_lon,region_max_lon] +#NOTE that the sub-region is added for V1 + def define_regions(InputRegion): - Region_Lookup={0 : ['not_assigned','not assigned',0,0,0,0], - 1 : ['nh_polar_west','NH Polar West',70.0,90.0,-135.0+360,-45.0+360], - 2 : ['nh_polar_east','NH Polar East',58.0,90.0,-45.0,45.0], - 3 : ['canada','Canada',48.0,62.0,-135.0+360,-45.0+360], - 4 : ['west_europe','Western Europe',43.0,57.5,0.0,25.0], - 5 : ['eastern_us','Eastern US',34.0,40.0,-95.0+360,-75.0+360], - 6 : ['japan','Japan',30.0,40.0,120.0,150.0], - 7 : ['nh_tropic','NH Sub Tropics',15.0,30.0,90.0,225.0], - 8 : ['tropics1','W-Pacific/E-Indian Ocean',-20.0,0.0,90.0,225.0], - 9 : ['tropics2','Equatorial Americas',-15.0,15.0,225.0,315.0], - 10 : ['tropics3','Atlantic/Africa',-15.0,15.0,-45.0,45.0], - 11 : ['sh_midlat','SH Midlatitudes',-57.5,-40.0,135.0,180.0], - 12 : ['sh_polar','SH Polar',-90.0,-58.0,0,360], - 13 : ['Boulder','Boulder',37.0,42.0,-110.0+360.0,-100.0+360.0], + Region_Lookup={0 : ['not_assigned','not assigned','not assigned',0,0,0,0], + 1 : ['nh_polar_west','NH Polar West','High Latitudes',70.0,90.0,-135.0+360,-45.0+360], + 2 : ['nh_polar_east','NH Polar East','High Latitudes',58.0,90.0,-45.0,45.0], + 3 : ['canada','Canada','High Latitudes',48.0,62.0,-135.0+360,-45.0+360], + 4 : ['west_europe','Western Europe','Mid Latitudes',43.0,57.5,0.0,25.0], + 5 : ['eastern_us','Eastern US','Mid Latitudes',34.0,40.0,-95.0+360,-75.0+360], + 6 : ['japan','Japan','Mid Latitudes',30.0,40.0,120.0,150.0], + 7 : ['nh_tropic','NH Sub Tropics','Tropics',15.0,30.0,90.0,225.0], + 8 : ['tropics1','W-Pacific/E-Indian Ocean','Tropics',-20.0,0.0,90.0,225.0], + 9 : ['tropics2','Equatorial Americas','Tropics',-15.0,15.0,225.0,315.0], + 10 : ['tropics3','Atlantic/Africa','Tropics',-15.0,15.0,-45.0,45.0], + 11 : ['sh_midlat','SH Midlatitudes','Mid Latitudes',-57.5,-40.0,135.0,180.0], + 12 : ['sh_polar','SH Polar','High Latitudes',-90.0,-58.0,0,360], + 13 : ['Boulder','Boulder','not assigned',37.0,42.0,-110.0+360.0,-100.0+360.0], } #-------------------------------------------------------------------------------------- @@ -213,14 +231,17 @@ def open_process_sonde_data_simone(obsdir): if (os.path.isfile(Check_File)): print("Using Ozone Climatology File: ",Check_File) - O3_Mean = xr.open_dataset(Check_File).o3_mean - O3_Width = xr.open_dataset(Check_File).o3_width - O3_StdDev = xr.open_dataset(Check_File).o3_std + + ds = xr.open_dataset(Check_File) + O3_Mean = ds.o3_mean + O3_Width = ds.o3_width + O3_StdDev = ds.o3_std + if (ifileID == 'Boulder'): - Pressure = xr.open_dataset(Check_File).press #hPa + Pressure = ds.press else: - Pressure = xr.open_dataset(Check_File).levels #hPa - + Pressure = ds.levels + if (i == 1): O3_MeanC=O3_Mean O3_WidthC=O3_Width @@ -237,6 +258,62 @@ def open_process_sonde_data_simone(obsdir): return O3_MeanC,O3_WidthC,O3_StdDevC,np.squeeze(PressureC),Region +#----------------------------------------------------------------------------------------- +#Process ozone sonde climatology V2 +#----------------------------------------------------------------------------------------- +def open_process_sonde_data_simone_v2(obsdir): + + #Grab and package all of the data that will be used for plotting + NRegions=define_regions(0)[1] + + O3_MeanC=[] + O3_WidthC=[] + O3_StdDevC=[] + Region=[] + PressureC=[] + + for i in range(1,NRegions+1): + Region_Info=define_regions(i)[0] + ifileID=Region_Info[0] + + Region.append(ifileID) + + Check_File=obsdir+'ozonesondes_'+ifileID+'1995_2011.nc' + + if (os.path.isfile(Check_File)): + print("Using Ozone Climatology File: ",Check_File) + + ds = xr.open_dataset(Check_File) + O3_Mean = ds.o3_mean + O3_Width = ds.o3_width + O3_StdDev = ds.o3_std + + if (ifileID == 'Boulder'): + Pressure = ds.press + O3_Mean = O3_Mean.rename({'level':'press'}) + O3_Width=O3_Width.rename({'level':'press'}) + O3_StdDev=O3_StdDev.rename({'level':'press'}) + Pressure = Pressure.rename({'level':'press'}) + else: + Pressure = ds.levels + + if (i == 1): + O3_MeanC=O3_Mean + O3_WidthC=O3_Width + O3_StdDevC=O3_StdDev + PressureC=Pressure + else: + O3_MeanC = xr.concat([O3_MeanC,O3_Mean],'z') + O3_WidthC = xr.concat([O3_WidthC,O3_Width],'z') + O3_StdDevC = xr.concat([O3_StdDevC,O3_StdDev],'z') + PressureC=xr.concat([PressureC,Pressure],'z') + + else: + print("Error opening ozonesonde file ",Check_File) + return -1 + + return O3_MeanC,O3_WidthC,O3_StdDevC,PressureC,Region + #----------------------------------------------------------------------------------------- #Subroutine - plots the ozone data subplots #----------------------------------------------------------------------------------------- @@ -301,9 +378,6 @@ def Subplot_O3(ax,X0,Y0,X1,Model0,Model1,SondeMean,SondeLow,SondeHigh,PLevel,plo #legend (bottom left) side of bottom left plot if Quadrant == 2: - props1 = dict(boxstyle='square',facecolor='white',alpha=0.5) - #ax.text(0.03,0.03," \n \n \n \n ",transform=ax.transAxes,fontsize=5.25,verticalalignment='bottom',bbox=props1,color='black') - if Compare_Obs >0: ax.text(0.12,0.15,"Region Avg.",fontsize=5.0,color='tab:red',transform=ax.transAxes) ax.text(0.12,0.09,case_test,fontsize=5.0,color='tab:red',transform=ax.transAxes) @@ -343,9 +417,6 @@ def Subplot_O3(ax,X0,Y0,X1,Model0,Model1,SondeMean,SondeLow,SondeHigh,PLevel,plo if Quadrant == 2: offset=0.67 - props1 = dict(boxstyle='square',facecolor='white',alpha=0.5) - #ax.text(0.03,0.03+offset," \n \n \n \n ",transform=ax.transAxes,fontsize=5.25,verticalalignment='bottom',bbox=props1,color='black') - if Compare_Obs > 0: ax.text(0.12,0.16+offset,"Ozonesondes",fontsize=5.0,color='black',transform=ax.transAxes) ax.plot([0.06],[0.17+offset],'ko',markersize=1.25,transform=ax.transAxes) @@ -426,7 +497,7 @@ def Plot_Seasonal_Cycle_Profile(X0,X01,Y0,X1,X11,Y1,X2,X22,Y2,X3,X33,Y3,SondeMea plt.clf() plt.close(fig) except: - print("ERROR: Could not save file ",Output_FName,flush=True) + print("ERROR: Could not save file ",Output_IMG,flush=True) return True @@ -439,14 +510,14 @@ class Model_Data: #Get the Reference Case data if (os.path.isfile(ClimoFile)): print("ADF Climatology File Located: ",ClimoFile) - lon = xr.open_dataset(ClimoFile).lon - lat = xr.open_dataset(ClimoFile).lat - #lev = xr.open_dataset(ClimoFile).lev - hyam = xr.open_dataset(ClimoFile).hyam - hybm = xr.open_dataset(ClimoFile).hybm - o3=xr.open_dataset(ClimoFile).O3 - ps=xr.open_dataset(ClimoFile).PS - #time=xr.open_dataset(ClimoFile).time + + ds = xr.open_dataset(ClimoFile) + lon = ds.lon + lat = ds.lat + hyam = ds.hyam + hybm = ds.hybm + o3 = ds.O3 + ps = ds.PS #these typically won't change with time so only select the time 0 hyam=hyam[0,:] @@ -466,35 +537,34 @@ class model_dat_proc: if (MinLon < 0 and MaxLon > 0): #if the region crosses the date line - do different processing O3_00=Model_Dat.o3.sel(lon=slice(MinLon+360.0,360.0),lat=slice(MinLat,MaxLat)) O3_01=Model_Dat.o3.sel(lon=slice(0,MaxLon),lat=slice(MinLat,MaxLat)) - O3_0 = np.concatenate( (O3_00,O3_01),axis=3) + O3_0 = xr.concat([O3_00, O3_01], dim='lon') #Changed for using geocat.comp.interpolation function PS_00=Model_Dat.ps.sel(lon=slice(MinLon+360.0,360.0),lat=slice(MinLat,MaxLat)) PS_01=Model_Dat.ps.sel(lon=slice(0,MaxLon),lat=slice(MinLat,MaxLat)) - PS_0 = np.concatenate( (PS_00,PS_01),axis=2) + PS_0 = xr.concat([PS_00, PS_01], dim='lon') #Changed for using geocat.comp.interpolation function lon_00=Model_Dat.lon.sel(lon=slice(MinLon+360.0,360.0)) lon_01=Model_Dat.lon.sel(lon=slice(0,MaxLon)) - lon_0 = np.concatenate( (lon_00,lon_01)) + lon_0 = xr.concat([lon_00, lon_01], dim='lon') #Changed for using geocat.comp.interpolation function #resort the arrays as needed - lon_sort=lon_0.argsort() - O3_0 = O3_0[:,:,:,lon_sort] - PS_0 = PS_0[:,:,lon_sort] - lon_0 = lon_0[lon_sort] - - O3_sfc=np.squeeze(O3_0[:,-1,:,:])*1.0e9 #get the lowest model surface level data + lon_sort=np.concatenate((lon_00,lon_01)).argsort() #Changed for using geocat.comp.interpolation function + O3_0 = O3_0.isel(lon=lon_sort) + PS_0 = PS_0.isel(lon=lon_sort) + lon_0 = lon_0.isel(lon=lon_sort) + + O3_sfc=O3_0.isel(lev=[-1]).squeeze()*1.0e9 #get the lowest model surface level data else: #if the region does not cross the date line O3_0=Model_Dat.o3.sel(lon=slice(MinLon,MaxLon),lat=slice(MinLat,MaxLat)) PS_0=Model_Dat.ps.sel(lon=slice(MinLon,MaxLon),lat=slice(MinLat,MaxLat)) lon_0=Model_Dat.lon.sel(lon=slice(MinLon,MaxLon)) - O3_sfc=np.squeeze(O3_0.values[:,-1,:,:])*1.0e9 - + O3_sfc=O3_0.isel(lev=[-1]).squeeze()*1.0e9 lat_0=Model_Dat.lat.sel(lat=slice(MinLat,MaxLat)) O3_sfc_04=np.mean(O3_sfc,axis=(1,2)) - - O3_0I = Ngl.vinth2p(O3_0,Model_Dat.hyam,Model_Dat.hybm,pnew,PS_0,intyp,1000.0,1,kxtrp)*1.0e9 + + O3_0I = gcomp.interpolation.interp_hybrid_to_pressure(O3_0,PS_0,Model_Dat.hyam,Model_Dat.hybm,new_levels=np.array(pnew),method='linear',p0=100000.0,extrapolate=False)*1.0e9 #Get the seasonal cycle of the base case at each needed pressure level #and average over the region. @@ -519,6 +589,12 @@ class model_dat_proc: months=[1,2,3,4,5,6,7,8,9,10,11,12] #set up the regular grid interpolator for each case and level + + O3_0I=O3_0I.values + lat_0=lat_0.values + lon_0=lon_0.values + O3_sfc=O3_sfc.values + interp_0 = RegularGridInterpolator((months,lat_0,lon_0), np.squeeze(O3_0I[:,0,:,:])) interp_1 = RegularGridInterpolator((months,lat_0,lon_0), np.squeeze(O3_0I[:,1,:,:])) interp_2 = RegularGridInterpolator((months,lat_0,lon_0), np.squeeze(O3_0I[:,2,:,:])) @@ -561,8 +637,9 @@ class model_dat_proc: def process_model_profiles(Model_Dat,O3_0,PS_0,pnew,intyp,kxtrp,ILAT,ILON,lat_0,lon_0): class model_dat_proc: - - O3_0I1 = Ngl.vinth2p(O3_0,Model_Dat.hyam,Model_Dat.hybm,pnew,PS_0,intyp,1000.0,1,kxtrp)*1.0e9 + + O3_0I1 = gcomp.interpolation.interp_hybrid_to_pressure(data=O3_0,hyam=Model_Dat.hyam,hybm=Model_Dat.hybm,new_levels=np.array(pnew*100.0),ps=PS_0,method='linear',p0=100000.0,extrapolate=False)*1.0e9 + O3_0I1=O3_0I1.values Locate_Bad=np.where(O3_0I1 > 10000.0) if len(Locate_Bad) > 0: @@ -619,6 +696,154 @@ class model_dat_proc: Processed_Model_Dat=model_dat_proc() return Processed_Model_Dat + +#----------------------------------------------------------------------------------------- +#Sets up the Taylor Diagram Plot +#----------------------------------------------------------------------------------------- +def taylor_plot_setup(ax,baseline): + """Constructs Figure and Axes objects for basic Taylor Diagram.""" + corr_labels = np.array([0.0, .1, .2, .3, .4, .5, .6, .7, .8, .9, .95, .99, 1.]) + corr_locations = np.pi/2 - np.arccos((corr_labels)) # azim. ticks in radians. + ax.set_thetamin(0) + ax.set_thetamax(90) + ax.set_ylim([0, 1.6]) # Works better than set_rmin / set_rmax + ax.set_theta_zero_location("N") # zero at top, + ax.set_theta_direction(-1) # angle increases clockwise + thetalines, thetalabels = ax.set_thetagrids(np.degrees(corr_locations), corr_labels) + ax.grid(axis='x', linewidth=0) # turn off radial grid + ax.set_rgrids(np.arange(-0.5, 2.75, .50)) + ax.set_ylabel("Difference from Mean [Normalized]") + # Add tick marks along azimuth + tick = [ax.get_rmax(),ax.get_rmax()*0.97] + for t in corr_locations: + ax.plot([t,t], tick, lw=0.72, color="k") + ax.text(np.radians(50), ax.get_rmax()*1.1, "Correlation", ha='center', rotation=-50, fontsize=8) + ax.text(np.radians(95), 1.0, "REF", ha='center') + ax.set_title(baseline, fontsize=8,pad=5) + + return True + +#----------------------------------------------------------------------------------------- +#Calculate the correlations, sigmas, and bias +#----------------------------------------------------------------------------------------- +def taylor_stats_single(casedata, refdata): + """This replicates the basic functionality of 'taylor_stats' from NCL. + input: + casedata : input data, DataArray + refdata : reference case data, DataArray + w : if true use cos(latitude) as spatial weight, if false assume uniform weight + returns: + pattern_correlation, ratio of standard deviation (case/ref), bias + """ + + correlation = pearsonr(casedata.values,refdata.values).statistic**2 #this pearson correlation seems to work the best + a_sigma = casedata.std() + b_sigma = refdata.std() + + mean_case=casedata.mean() + mean_ref=refdata.mean() + + mean_diff = 1.0 + (mean_case - mean_ref) / mean_ref + + bias = (100*((mean_case - mean_ref)/mean_ref)).item() + + ratio = a_sigma/b_sigma + if ratio > 2.5: + ratio=2.5 + + return correlation, ratio, bias, mean_diff + +#----------------------------------------------------------------------------------------- +#Plot the Taylor Data +#----------------------------------------------------------------------------------------- +def plot_taylor_data(wks, df,k, **kwargs): + """Apply data on top of the Taylor Diagram Axes. + wks -> Axes object, probably from taylor_plot_setup + df -> DataFrame holding the Taylor stats. + kwargs -> optional arguments + look for 'use_bias' + look for 'case_color' + """ + # option is whether to stylize the markers by the bias: + use_bias = False + if 'use_bias' in kwargs: + if kwargs['use_bias']: + use_bias = True + df['bias_digi'] = np.digitize(df['bias'].values, [-20, -10, -5, -1, 1, 5, 10, 20]) + marker_list = ["v", "v", "v", "v", "o", "^", "^", "^", "^"] + marker_size = [8, 4, 2, 1, 1, 1, 2, 4, 8] + # option: has color been specified as case_color? + # --> expect the case labeling to be done external to this function + if 'case_color' in kwargs: + color = kwargs['case_color'] + if isinstance(color, int): + # assume we should use this as an index + color = mpl.cm.tab20(color) # convert to RGBA + # TODO: allow colormap to be specified. + + #base case + annos = [] # list will hold strings for legend + #k = 1 + for ndx, row in df.iterrows(): + # NOTE: ndx will be the DataFrame index, and we expect that to be the variable name + theta = np.pi/2 - np.arccos(row['corr']) # Transform DATA + if use_bias: + mk = marker_list[row['bias_digi']] + mksz = marker_size[row['bias_digi']] + wks.plot(theta, row['mean_diff'], marker=mk, markersize=mksz, color=color) + else: + wks.plot(theta, row['mean_diff'], marker='o', markersize=6, color=color) + annos.append(f"{k} - {ndx.replace('_','')}") + wks.annotate(str(k), (theta, row['mean_diff']), ha='center', va='bottom',xytext=(0,5), textcoords='offset points', fontsize='x-large',color=color) + k += 1 # increment the annotation number (THIS REQUIRES CASES TO HAVE SAME ORDER IN DataFrame) + + return wks + +#----------------------------------------------------------------------------------------- +#Plot out the final components of the Taylor Diagrams +#----------------------------------------------------------------------------------------- +def taylor_plot_finalize(wks, casenames, casecolors, syear_cases, eyear_cases, SubZones, needs_bias_labels=False,needs_case_labels=False,needs_subzone_labels=False): + """Apply final formatting to a Taylor diagram. + wks -> Axes object that has passed through taylor_plot_setup and plot_taylor_data + casenames -> list of case names for the legend + casecolors -> list of colors for the cases + needs_bias_labels -> Bool, if T make the legend for the bias-sized markers. + """ + # CASE LEGEND -- Color-coded + bottom_of_text = 0.05 + height_of_lines = 0.08 + + n = 0 + + #CASE LEGEND + if needs_case_labels: + for case_idx, (s, c) in enumerate(zip(casenames, casecolors)): + text = wks.text(0.052, bottom_of_text + n*height_of_lines, s, + color=c, ha='left', va='bottom', transform=wks.transAxes, fontsize=8) + n += 1 + + # BIAS LEGEND + if needs_bias_labels: + # produce an info-box showing the markers/sizes based on bias + bias_legend_elements = [(Line2D([0], [0], marker="v", color='k', label="> 20%", markersize=8, fillstyle='none', linewidth=0), Line2D([0], [0], marker="^", color='k', label="> 20%", markersize=8, fillstyle='none', linewidth=0)), + (Line2D([0], [0], marker="v", color='k', label="10-20%", markersize=4, linewidth=0), Line2D([0], [0], marker="^", color='k', label="10-20%", markersize=4, linewidth=0)), + (Line2D([0], [0], marker="v", color='k', label="5-10%", markersize=2, linewidth=0), Line2D([0], [0], marker="^", color='k', label="5-10%", markersize=2, linewidth=0)), + (Line2D([0], [0], marker="v", color='k', label=">1-5%", markersize=1, linewidth=0), Line2D([0], [0], marker="^", color='k', label=">1-5%", markersize=1, linewidth=0)), + Line2D([0], [0], marker="o", color='k', label="< 1%", markersize=1, linewidth=0), + ] + bias_legend_labels = ["> 20%", "10-20%", "5-10%", "1-5%", "< 1%"] + wks.legend(handles=bias_legend_elements, labels=bias_legend_labels, loc='lower left', handler_map={tuple: HandlerTuple(ndivide=None, pad=1.)}, labelspacing=1, handletextpad=1, frameon=False, title="-/+ Bias", + title_fontsize=10) + + #SubZone Legend + n = 0 + if needs_subzone_labels: + for k in range(0,len(SubZones)): + text=wks.text(0.052,bottom_of_text + n*height_of_lines,str(k+1)+' '+SubZones[k], + ha='left', va='bottom', transform=wks.transAxes, fontsize=8) + n += 1 + + return wks #----------------------------------------------------------------------------------------- #Primary Ozone Diagnostics Routine @@ -631,7 +856,6 @@ def ozone_diagnostics (adfobj): obsdir='/glade/campaign/acom/acom-climate/tilmes/obs_data/amwg/amwg_data/obs_data/cam-chem/' #location of the ozonesonde data files cam_climo_loc = adfobj.get_cam_info('cam_climo_loc',required=True)[0] baseline_climo_loc = adfobj.get_baseline_info('cam_climo_loc',required=True) - ncases = len(cam_climo_loc) #check the number of cases. For the O3 diagnostics there needs to be 2. #If there is more than 2, then only the first two will be used. @@ -662,9 +886,11 @@ def ozone_diagnostics (adfobj): redo_plot = adfobj.get_basic_info('redo_plot') print(f"\t NOTE: redo_plot is set to {redo_plot}") + #Compare_Obs=0 #initially assumes that user is comparing two models Compare_Obs=0 #initially assumes that user is comparing two models if adfobj.compare_obs: Compare_Obs=1 #if comparing a model with observations then don't need to include the base case in the plots + Compare_Obs=adfobj.compare_obs print("-------------------------------") print("Processing Ozone Diagnostics...") @@ -719,6 +945,7 @@ def ozone_diagnostics (adfobj): #set the parameters for interpolation to pressure levels #-------------------------------------------------------------------------------------- pnew = [50,250,500,900] + pnew_Pa = np.array(pnew)*100.0 intyp = 2 # 1=linear, 2=log, 3=log-log kxtrp = False # True=extrapolate (when the output pressure level is outside of the range of psrf) months=[1,2,3,4,5,6,7,8,9,10,11,12] @@ -731,15 +958,17 @@ def ozone_diagnostics (adfobj): #----------------------------------------------------------------------------------- #Grab the info for the current region of interest #----------------------------------------------------------------------------------- - Region_Info=define_regions(i)[0] + Region_Info=define_regions(i)[0] #Get information on the current region of interest SName=Region_Info[0] LName=Region_Info[1] - MinLat=Region_Info[2] - MaxLat=Region_Info[3] - MinLon=Region_Info[4] - MaxLon=Region_Info[5] - oFile_Seasonal = plot_locations+'/O3SeasonalCycle_'+SName+'_Special.png' - oFile_Profile = plot_locations+'/O3Profile_'+SName+'_Special.png' + SubRegion=Region_Info[2] + MinLat=Region_Info[3] + MaxLat=Region_Info[4] + MinLon=Region_Info[5] + MaxLon=Region_Info[6] + + oFile_Seasonal = plot_locations+'/'+SName.replace("_","")+'_SeasonalCycle_ANN_Special_Mean.png' + oFile_Profile = plot_locations+'/'+SName.replace("_","")+'_Profile_ANN_Special_Mean.png' #----------------------------------------------------------------------------------- #Check if redo_plot set and if not and plots exist already then @@ -755,9 +984,9 @@ def ozone_diagnostics (adfobj): #----------------------------------------------------------------------------------- #Process the model data #----------------------------------------------------------------------------------- - Processed_Seasonal_Cycle_Test_Data = process_model_seasonal_cycle(MinLon,MaxLon,MinLat,MaxLat,Test_Data,pnew,intyp,kxtrp,Station_Lons,Station_Lats) + Processed_Seasonal_Cycle_Test_Data = process_model_seasonal_cycle(MinLon,MaxLon,MinLat,MaxLat,Test_Data,pnew_Pa,intyp,kxtrp,Station_Lons,Station_Lats) if Compare_Obs <= 0: - Processed_Seasonal_Cycle_Base_Data = process_model_seasonal_cycle(MinLon,MaxLon,MinLat,MaxLat,Base_Data,pnew,intyp,kxtrp,Station_Lons,Station_Lats) + Processed_Seasonal_Cycle_Base_Data = process_model_seasonal_cycle(MinLon,MaxLon,MinLat,MaxLat,Base_Data,pnew_Pa,intyp,kxtrp,Station_Lons,Station_Lats) else: Processed_Seasonal_Cycle_Base_Data=Processed_Seasonal_Cycle_Test_Data #allows the plotting routine to run without having to call it different ways. @@ -768,7 +997,6 @@ def ozone_diagnostics (adfobj): Obs_Width_0=Obs_Width[:,:,i-1] Obs_StdDev_0=Obs_StdDev[:,:,i-1] Obs_Pressure_0=Obs_Pressure[:,i-1] - Regions_0=Regions[i-1] #----------------------------------------------------------------------------------- #Get the pressure level index for each of the 4 levels to @@ -793,10 +1021,11 @@ def ozone_diagnostics (adfobj): Obs_Mean_20=np.squeeze(Obs_Mean_0[Locate_P2,:]) Obs_Mean_30=np.squeeze(Obs_Mean_0[Locate_P3,:]) - Obs_Width_00=np.squeeze(Obs_Width_0[Locate_P0,:]) - Obs_Width_10=np.squeeze(Obs_Width_0[Locate_P1,:]) - Obs_Width_20=np.squeeze(Obs_Width_0[Locate_P2,:]) - Obs_Width_30=np.squeeze(Obs_Width_0[Locate_P3,:]) + #Currently not being utilized + #Obs_Width_00=np.squeeze(Obs_Width_0[Locate_P0,:]) + #Obs_Width_10=np.squeeze(Obs_Width_0[Locate_P1,:]) + #Obs_Width_20=np.squeeze(Obs_Width_0[Locate_P2,:]) + #Obs_Width_30=np.squeeze(Obs_Width_0[Locate_P3,:]) Obs_StdDev_00=np.squeeze(Obs_StdDev_0[Locate_P0,:]) Obs_StdDev_10=np.squeeze(Obs_StdDev_0[Locate_P1,:]) @@ -841,10 +1070,11 @@ def ozone_diagnostics (adfobj): Obs_Mean_201=np.squeeze(Obs_Mean_0[:,6]) Obs_Mean_301=np.squeeze(Obs_Mean_0[:,9]) - Obs_Width_001=np.squeeze(Obs_Width_0[:,0]) - Obs_Width_101=np.squeeze(Obs_Width_0[:,3]) - Obs_Width_201=np.squeeze(Obs_Width_0[:,6]) - Obs_Width_301=np.squeeze(Obs_Width_0[:,9]) + #Currently not being utilized + #Obs_Width_001=np.squeeze(Obs_Width_0[:,0]) + #Obs_Width_101=np.squeeze(Obs_Width_0[:,3]) + #Obs_Width_201=np.squeeze(Obs_Width_0[:,6]) + #Obs_Width_301=np.squeeze(Obs_Width_0[:,9]) Obs_StdDev_001=np.squeeze(Obs_StdDev_0[:,0]) Obs_StdDev_101=np.squeeze(Obs_StdDev_0[:,3]) @@ -869,7 +1099,174 @@ def ozone_diagnostics (adfobj): #----------------------------------------------------------------------------------- #Once the plots have successfully run, add the web page entries (if enabled). #----------------------------------------------------------------------------------- - adfobj.add_website_data(oFile_Seasonal,SName.replace("_","")+"_SeasonalCycle", None, season="ANN",multi_case=True,category="O3_DIAGNOSTICS") - adfobj.add_website_data(oFile_Profile,SName.replace("_","")+"_Profile", None, season="ANN", multi_case=True,category="O3_DIAGNOSTICS") + adfobj.add_website_data(oFile_Seasonal,SName.replace("_","")+"_SeasonalCycle", None, season="ANN",multi_case=True,category="Ozone Diagnostics - Seasonal Cycle") + adfobj.add_website_data(oFile_Profile,SName.replace("_","")+"_Profile", None, season="ANN", multi_case=True,category="Ozone Diagnostics - Profiles") + + #-------------------------------------------------------------------------------------- + #-------------------------------------------------------------------------------------- + #Plot out Taylor Diagrams + # + # Plot a taylor diagram at each pressure level and zone. + # + # For this taylor diagram we are trying to compare ozonesonde data to model outputs. + # + # The ozonesonde files located at /glade/campaign/acom/acom-climate/tilmes/obs_data/amwg/amwg_data/obs_data/cam-chem/ + # condaint the mean sonde profiles (press,months). These need to get compared with + # model data which is (months,press,lat,lon). Specific pressure levels are selected. + #-------------------------------------------------------------------------------------- + #-------------------------------------------------------------------------------------- + + #-------------------------------------------------------------------------------------- + #Set plotting level and zone information + #-------------------------------------------------------------------------------------- + PLEVS=[900,500,250,50] #pressure levels to plot + + Zones=['Tropics','Mid Latitudes','High Latitudes'] #zones to plot + + SubZones= {'Tropics': ['NH SubTropic','W-Pacific/E-Indian Ocean','equatorial Americas','Atlantic/Africa'], + 'Mid Latitudes': ['Western Europe','Eastern US','Japan','SH Midlatitudes'], + 'High Latitudes': ['NH Polar West','NH Polar East','Canada','SH Polar']} + + SubZoneIndexes= {'Tropics': [6,7,8,9], + 'Mid Latitudes': [3,4,5,10], + 'High Latitudes': [0,1,2,11]} + + #Define the output filename and location of the image file + Output_IMG=plot_locations+"/OzoneSonde_TaylorDiagram_ANN_Special_Mean.png" + + #----------------------------------------------------------------------------------- + #Check if redo_plot set and if not and plots exist already then + #add them to the website (if enabled) and then exit the routine + #----------------------------------------------------------------------------------- + if (not(redo_plot)) and (os.path.isfile(Output_IMG)): + print(SName,' Ozone Taylor plot exists and redo_plot is false. Adding to website and Skipping plot.') + adfobj.add_website_data(Output_IMG, "OzoneSonde_TaylorDiagram", None,season="ANN", multi_case=True) + sys.exit(0) + else: + print("Plotting Ozone Taylor Diagrams") + + #-------------------------------------------------------------------------------------- + #Get the ozone sonde and related data + #-------------------------------------------------------------------------------------- + SData,SData_Width,SData_Dev,SDataPrs,SData_Region = open_process_sonde_data_simone_v2(obsdir) + SDataPrs=SDataPrs[0,:] + + #initialize subplots + fig, axs = plt.subplots(nrows=len(PLEVS),ncols=len(Zones),sharex=False,sharey=False,subplot_kw={'projection':'polar'},figsize=(12,12)) + + #loop through each pressure level (900, 500, 250, 50) + for i in range(0,len(PLEVS)): + + CLEV=float(PLEVS[i]) #set the current pressure level + #loop through each ZONES (i.e. Tropics, Mid Latitudes, and High Latitudes) + for j in range(0,len(Zones)): + + #set up the basic plot + TayX=taylor_plot_setup(axs[i,j],Zones[j]+" "+str(PLEVS[i])+" hPa") + + #loop through each subzone and average the ozonesonde data + ListSubZones=SubZones[Zones[j]] + ListSubZonel=SubZoneIndexes[Zones[j]] + + base_x = SData[ListSubZonel,:,:] #limit to only the subregions within the main regions + + Find_Index=np.where(SDataPrs == CLEV) #find the right pressure index + + base_x=base_x[:,Find_Index[0],:] #limit to the required pressure level - AT THIS POINT BASE_X HAS SIZE (NSONDE,1,12) + + SubZoneIndex0=SubZoneIndexes[Zones[j]] + + #loop through each sub zone and plot on the relavent pressure level and primary zone plot + for k in range(0,len(SubZoneIndex0)): + + #---------------------------------------------------------------------------- + #determine the latitude and longitude ranges to average model data for + #---------------------------------------------------------------------------- + Region_Info=define_regions(SubZoneIndex0[k]+1)[0] + SName=Region_Info[0] + LName=Region_Info[1] + SubRegion=Region_Info[2] + MinLat=float(Region_Info[3]) + MaxLat=float(Region_Info[4]) + MinLon=float(Region_Info[5]) + MaxLon=float(Region_Info[6]) + + if (MinLon < 0 and MaxLon > 0): #if the region crosses the date line - do different processing + O3_00=Base_Data.o3.sel(lon=slice(MinLon+360.0,360.0),lat=slice(MinLat,MaxLat)) + O3_01=Base_Data.o3.sel(lon=slice(0,MaxLon),lat=slice(MinLat,MaxLat)) + O3_0 = xr.concat([O3_00,O3_01],'lon') + O3_0 = O3_0.sel(lev=CLEV,method='nearest') + O3_10=Test_Data.o3.sel(lon=slice(MinLon+360.0,360.0),lat=slice(MinLat,MaxLat)) + O3_11=Test_Data.o3.sel(lon=slice(0,MaxLon),lat=slice(MinLat,MaxLat)) + O3_1 = xr.concat([O3_10,O3_11],'lon') + O3_1 = O3_1.sel(lev=CLEV,method='nearest') + + else: #if the region does not cross the date line + O3_0=Base_Data.o3.sel(lon=slice(MinLon,MaxLon),lat=slice(MinLat,MaxLat)) + O3_0=O3_0.sel(lev=CLEV,method='nearest') + O3_1=Test_Data.o3.sel(lon=slice(MinLon,MaxLon),lat=slice(MinLat,MaxLat)) + O3_1=O3_1.sel(lev=CLEV,method='nearest') + + lat_0=Base_Data.lat.sel(lat=slice(MinLat,MaxLat)) + lat_0=lat_0.mean(dim='lat') + lat_1=Test_Data.lat.sel(lat=slice(MinLat,MaxLat)) + lat_1=lat_1.mean(dim='lat') + + base_x0=base_x[k,:,:] + base_x0=base_x0.squeeze('press') + + test_x_base=O3_0*1.0e9 #convert mol/mol to ppbv + test_x_test=O3_1*1.0e9 #convert mol/mol to ppbv + + #average the regions down to monthly average values that can be compared with the monthly average sonde values. + test_x_base=test_x_base.mean(dim='lat') + test_x_test=test_x_test.mean(dim='lat') + test_x_base=test_x_base.mean(dim='lon') + test_x_test=test_x_test.mean(dim='lon') + + #set up the data templates to pass to the taylor stats calculator + df_template = pd.DataFrame(index=['O3'], columns=['corr', 'ratio', 'bias','mean_diff']) + result_by_case = {cname: df_template.copy() for cname in ['base','test']} + + #calculate the taylor stats for the base and test case + result_by_case['base'].loc['O3'] = taylor_stats_single(test_x_base, base_x0) + result_by_case['test'].loc['O3'] = taylor_stats_single(test_x_test, base_x0) + + case_colors = [mpl.cm.tab20(0),mpl.cm.tab20(6)] #hard code to blue and red to match the original ncl version + + #Set plotting flags that determine which panels to plot things like legends. + if (i == 0 and j == 0 and k == 0): + Plot_Bias_Legend=False + else: + Plot_Bias_Legend=False + + if (i == 0 and j == 1 and k == 0): + Plot_Case_Legend=True + else: + Plot_Case_Legend=False + + if (i == 3): + Plot_Zone_Legend=True + else: + Plot_Zone_Legend=False + + #plot the taylor data for each case + axs0 = plot_taylor_data(axs[i,j],result_by_case['base'],k+1, use_bias=False,case_color=case_colors[0]) + axs1 = plot_taylor_data(axs[i,j],result_by_case['test'],k+1, use_bias=False,case_color=case_colors[1]) + + #finalize the data plots with legends, labels, etc. + syear_cases = adfobj.climo_yrs["syears"] + eyear_cases = adfobj.climo_yrs["eyears"] + taylor_plot_finalize(axs[i,j],[case_base,case_test],case_colors, syear_cases, eyear_cases, ListSubZones, needs_bias_labels=Plot_Bias_Legend,needs_case_labels=Plot_Case_Legend,needs_subzone_labels=Plot_Zone_Legend) + + plt.subplots_adjust(wspace=-0.3,hspace=0.3) #Adjusts the spacing between the plots + st = fig.suptitle("Comparison to Ozonesondes", fontsize=16) + st.set_y(0.95) + + plt.savefig(Output_IMG, bbox_inches = 'tight', dpi = 150, format = 'png') #there will be more to add here... + + #Add plot to website (if enabled): + adfobj.add_website_data(Output_IMG, "OzoneSonde_TaylorDiagram",None, season="ANN",multi_case=True,category="Ozone Diagnostics - Taylor Diagrams") + print("Ozone Diagnostics Generated Successfully!") From 7f74b281932d1574cb4840d378f6cb6303dd587a Mon Sep 17 00:00:00 2001 From: shawnusaf <84995386+shawnusaf@users.noreply.github.com> Date: Wed, 23 Oct 2024 15:05:54 -0400 Subject: [PATCH 4/5] Update config_cam_baseline_example.yaml --- config_cam_baseline_example.yaml | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/config_cam_baseline_example.yaml b/config_cam_baseline_example.yaml index 391b61f57..c4f680d4c 100644 --- a/config_cam_baseline_example.yaml +++ b/config_cam_baseline_example.yaml @@ -182,6 +182,20 @@ diag_cam_climo: #Location where time series files are (or will be) stored: cam_ts_loc: /glade/derecho/scratch/${user}/ADF/${diag_cam_climo.cam_case_name}/ts + #TEM diagnostics + #--------------- + #TEM history file number + #If missing or blank, ADF will default to h4 + tem_hist_str: cam.h4 + + #Location where TEM files are stored: + #NOTE: If path not specified or commented out, TEM calculation/plots will be skipped! + cam_tem_loc: /glade/derecho/scratch/${user}/${diag_cam_climo.cam_case_name}/tem/ + + #Overwrite TEM files, if found? + #If set to false, then TEM creation will be skipped if files are found: + overwrite_tem: false + #---------------------- #You can alternatively provide a list of cases, which will make the ADF @@ -303,6 +317,20 @@ diag_cam_baseline_climo: #Location where time series files are (or will be) stored: cam_ts_loc: /glade/derecho/scratch/${user}/ADF/${diag_cam_baseline_climo.cam_case_name}/ts + #TEM diagnostics + #--------------- + #TEM history file number + #If missing or blank, ADF will default to h4 + tem_hist_str: cam.h4 + + #Location where TEM files are stored: + #NOTE: If path not specified or commented out, TEM calculation/plots will be skipped! + cam_tem_loc: /glade/derecho/scratch/${user}/${diag_cam_baseline_climo.cam_case_name}/tem/ + + #Overwrite TEM files, if found? + #If set to false, then TEM creation will be skipped if files are found: + overwrite_tem: false + #This fourth set of variables provides settings for calling the Climate Variability # Diagnostics Package (CVDP). If cvdp_run is set to true the CVDP will be set up and # run in background mode, likely completing after the ADF has completed. From a912838373b77e10730e9f404df9029b764256ef Mon Sep 17 00:00:00 2001 From: shawnusaf <84995386+shawnusaf@users.noreply.github.com> Date: Wed, 30 Oct 2024 15:21:03 -0400 Subject: [PATCH 5/5] Update config_cam_baseline_example.yaml --- config_cam_baseline_example.yaml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/config_cam_baseline_example.yaml b/config_cam_baseline_example.yaml index 97f2e861c..cc9bce430 100644 --- a/config_cam_baseline_example.yaml +++ b/config_cam_baseline_example.yaml @@ -58,7 +58,7 @@ # Note that the string 'USER-NAME-NOT-SET' is used in the jupyter script # to check for a failure to customize # -user: 'shawnh' +user: 'USER-NAME-NOT-SET' #This first set of variables specify basic info used by all diagnostic runs: @@ -530,4 +530,4 @@ diag_var_list: # region_season: # region_variables: -#END OF FILE \ No newline at end of file +#END OF FILE