Add V0 Ozone Diagnostics #285

shawnusaf · 2024-02-22T18:00:11Z

Adds V0 ozone diagnostics that output ozone seasonal cycle and profile plots with model comparisons to ozone sonde climatology.

Adds in V0 of the Ozone diagnostics package

shawnusaf · 2024-02-22T18:00:51Z

Addresses #281

justin-richling

Thanks for the work @shawnusaf, this is great to have. I just have some small changes for you to fix after which should be good to go.

The biggest change requested is the move to using geocat for the hybrid to pressure interpolation. I have the updated code in the review, but let me know if you have any questions about it.

Let me know when you get the changes in and well move forward, thanks!

config_cam_baseline_example.yaml

justin-richling · 2024-05-23T19:03:09Z

scripts/plotting/ozone_diagnostics.py

+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)


It might be good to have some one line descriptions for these functions here and maybe remove the arguments for clarity?

scripts/plotting/ozone_diagnostics.py

justin-richling · 2024-05-23T19:50:41Z

scripts/plotting/ozone_diagnostics.py

+   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


I might change this to a logical boolean:

Suggested change

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

Then you can change where ever the Compare_Obs check happens

justin-richling · 2024-05-23T19:52:24Z

scripts/plotting/ozone_diagnostics.py

+   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:


If you change Compare_Obs to a boolean (see comment in ozone_diagnostics section), then all the checks equal to 0 can be turned into the following:

Suggested change

if Compare_Obs <= 0:

if not Compare_Obs:

justin-richling · 2024-05-24T16:12:32Z

scripts/plotting/ozone_diagnostics.py

+   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)


The variable ncases doesn't seem to be used. You can remove it if you want.

justin-richling · 2024-05-28T15:00:50Z

scripts/plotting/ozone_diagnostics.py

+            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


We can use the geocat.comp.interpolation interp_hybrid_to_pressure function to replace the deprecated Ngl.vinth2p function. You'll just need to make some small changes to keep everything in xarray data arrays:

Suggested change

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

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 = xr.concat([O3_00, O3_01], dim='lon')

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 = xr.concat([PS_00, PS_01], dim='lon')

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 = xr.concat([lon_00, lon_01], dim='lon')

#resort the arrays as needed

lon_sort=np.concatenate((lon_00,lon_01)).argsort()

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=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 = gcomp.interpolation.interp_hybrid_to_pressure(data=O3_0,hyam=Model_Dat.hyam,hybm=Model_Dat.hybm,

new_levels=np.array(pnew),ps=PS_0,

method='linear',p0=1000.0)*1.0e9

justin-richling · 2024-05-28T15:03:05Z

scripts/plotting/ozone_diagnostics.py

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


Just like above, we can change the Ngl.vinth2p to geocat:

Suggested change

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),ps=PS_0,

method='linear',p0=1000.0)*1.0e9

justin-richling · 2024-05-29T17:12:17Z

scripts/plotting/ozone_diagnostics.py

+import os
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+import Ngl


Instead of using the Ngl package, which is deprecated, the ADF can bring in geocat.comp which provides the equivalent functionality. See the comments further down, they will show how to make the switch:

Suggested change

import Ngl

import geocat.comp as gcomp

justin-richling · 2024-05-31T20:13:34Z

scripts/plotting/ozone_diagnostics.py

+        plt.clf() 
+        plt.close(fig)
+   except:
+        print("ERROR: Could not save file ",Output_FName,flush=True)


Output_FName doesn't seem to exist, should it be changed to Output_IMG?

shawnusaf · 2024-10-11T14:13:08Z

Changes have been addressed in the newly updated ozone_diagnostics.py file. I have also added in taylor diagrams to the code.

justin-richling · 2024-10-30T19:18:28Z

config_cam_baseline_example.yaml

-user: 'USER-NAME-NOT-SET'
+user: 'shawnh'


This just needs to be changed back

justin-richling · 2024-10-30T19:18:42Z

scripts/plotting/ozone_diagnostics.py

+   case_base = adfobj.get_baseline_info('cam_case_name',required=True)
+
+   #Grab all case nickname(s)
+   test_nicknames = adfobj.case_nicknames["test_nicknames"]


@shawnusaf I would change this line so the plot legend isn't a list for the test case

justin-richling · 2024-10-30T19:23:52Z

@shawnusaf I left just a couple of suggestions, after that it should be good to merge in.

Shawn Honomichl and others added 3 commits January 18, 2024 11:26

V0

3fbe5d0

Merge branch 'NCAR:main' into main

c11ba78

O3 DIAG V0

85eadad

Adds in V0 of the Ozone diagnostics package

justin-richling added plotting Related to plot generation chemistry Related to chemistry diagnostics labels Mar 5, 2024

This was referenced May 8, 2024

Add climatological comparison of modeled ozone to ozonesonde observations #281

Open

Wishlist of chemistry obs datasets #237

Open

justin-richling requested changes May 31, 2024

View reviewed changes

Add files via upload

df801b7

shawnusaf added 2 commits October 23, 2024 15:05

Update config_cam_baseline_example.yaml

7f74b28

Merge branch 'main' into main

013b388

justin-richling reviewed Oct 30, 2024

View reviewed changes

Update config_cam_baseline_example.yaml

a912838

justin-richling merged commit b04efe2 into NCAR:main Oct 31, 2024

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add V0 Ozone Diagnostics #285

Add V0 Ozone Diagnostics #285

shawnusaf commented Feb 22, 2024

shawnusaf commented Feb 22, 2024

justin-richling left a comment

justin-richling May 23, 2024

justin-richling May 23, 2024

justin-richling May 23, 2024

justin-richling May 24, 2024

justin-richling May 28, 2024

justin-richling May 28, 2024

justin-richling May 29, 2024

justin-richling May 31, 2024

shawnusaf commented Oct 11, 2024

justin-richling Oct 30, 2024

justin-richling Oct 30, 2024

justin-richling commented Oct 30, 2024


		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 = 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),ps=PS_0,
+                                                              method='linear',p0=1000.0)*1.0e9

		user: 'USER-NAME-NOT-SET'
		user: 'shawnh'

Add V0 Ozone Diagnostics #285

Add V0 Ozone Diagnostics #285

Conversation

shawnusaf commented Feb 22, 2024

shawnusaf commented Feb 22, 2024

justin-richling left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

shawnusaf commented Oct 11, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

justin-richling commented Oct 30, 2024