Add V0 Ozone Diagnostics #285
Open
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Adds in V0 of the Ozone diagnostics package
|
Addresses #281 |
justin-richling
added
plotting
This was referenced May 8, 2024
justin-richling
requested changes
May 31, 2024
There was a problem hiding this comment.
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!
| @@ -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' | |||
There was a problem hiding this comment.
I would undo this change
Comment on lines
+11
to
+30
| 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) |
There was a problem hiding this comment.
It might be good to have some one line descriptions for these functions here and maybe remove the arguments for clarity?
| #----------------------------------------------------------------------------------------- | ||
| #Process ozone sonde climatology | ||
| #----------------------------------------------------------------------------------------- | ||
| def open_process_sonde_data_simone(obsdir): |
There was a problem hiding this comment.
It might be good to have a docstring describing what this function does. The same suggestion holds for all other functions below.
Comment on lines
+665
to
+667
| 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 |
There was a problem hiding this comment.
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
| 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: |
There was a problem hiding this comment.
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: |
| 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) |
There was a problem hiding this comment.
The variable ncases doesn't seem to be used. You can remove it if you want.
Comment on lines
+467
to
+497
| 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 |
There was a problem hiding this comment.
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 |
|
|
||
| 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 |
There was a problem hiding this comment.
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 |
| import os | ||
| import matplotlib.pyplot as plt | ||
| import matplotlib as mpl | ||
| import Ngl |
There was a problem hiding this comment.
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 |
| plt.clf() | ||
| plt.close(fig) | ||
| except: | ||
| print("ERROR: Could not save file ",Output_FName,flush=True) |
There was a problem hiding this comment.
Output_FName doesn't seem to exist, should it be changed to Output_IMG?
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Adds V0 ozone diagnostics that output ozone seasonal cycle and profile plots with model comparisons to ozone sonde climatology.