Skip to content

Commit

Permalink
add plot sondes obs
Browse files Browse the repository at this point in the history
  • Loading branch information
weihuang-jedi committed Nov 30, 2023
1 parent 1cf3166 commit cc9af07
Show file tree
Hide file tree
Showing 4 changed files with 634 additions and 30 deletions.
85 changes: 85 additions & 0 deletions panel2.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,85 @@
import os, sys
import matplotlib.pyplot as plt
from matplotlib import image
import numpy as np
import getopt

import matplotlib as mpl
mpl.rcParams['figure.dpi']= 300

#=========================================================================
class MakePanelPlot():
def __init__(self, debug=0, output=0, img1='sfcship_ps',
img2='surface_pressure'):
self.debug = debug
self.output = output
self.img1 = img1
self.img2 = img2

if(self.debug):
print('debug = ', debug)

if(self.debug > 10):
print('self.img1 = ', self.img1)
print('self.img2 = ', self.img2)

self.genit()

def genit(self):
#initialise grid of axes
#fig, axes = plt.subplots(2, 2, constrained_layout=True)
fig, axes = plt.subplots(1, 2)
axes = axes.ravel()

#create fake data
img = [
'FAKENAME_30.png',
'FAKENAME_40.png',
]

#iterate over axes
for i, ax in enumerate(axes):
#imageprefix = '%s_%s_increment_lev' %(img1, img2)
imageprefix = '%s_lev' %(img1)
imgname = img[i].replace('FAKENAME', imageprefix)
im = image.imread(imgname)
ax.axis('off')
ax.imshow(im)

plt.tight_layout()

if(self.output):
imgname = 'panelplot'
#fig.savefig(imgname, dpi=300, bbox_inches="tight")
fig.savefig(imgname, bbox_inches="tight", quality=100)
else:
plt.show()

#------------------------------------------------------------------------------
if __name__ == '__main__':
debug = 1
output = 0
img1 = 'sfcship_ps'
img2 = 'surface_pressure'

opts, args = getopt.getopt(sys.argv[1:], '', ['debug=', 'output=',
'img1=', 'img2='])

for o, a in opts:
if o in ('--debug'):
debug = int(a)
elif o in ('--output'):
output = int(a)
elif o in ('--img1'):
img1 = a
elif o in ('--img2'):
img2 = a
#else:
# assert False, 'unhandled option'

print('debug = ', debug)
print('output = ', output)

mpp = MakePanelPlot(debug=debug, output=output,
img1=img1, img2=img2)

64 changes: 34 additions & 30 deletions plot-compare-incr.py
Original file line number Diff line number Diff line change
Expand Up @@ -138,30 +138,33 @@ def set_cmapname(self, cmapname):
if __name__== '__main__':
debug = 1
output = 0
gsifile = '/work2/noaa/gsienkf/weihuang/jedi/per_core_timing/run/solver_halo_aircraft/run_80.40t1n_36p/analysis/increment/xainc.20200110_030000z.nc4'
#jedifile = '/work2/noaa/gsienkf/weihuang/jedi/per_core_timing/run/solver_RoundRobin_aircraft/run_80.40t1n_36p/analysis/increment/xainc.20200110_030000z.nc4'
jedifile = '/work2/noaa/gsienkf/weihuang/jedi/case_study/aircraft/run_80.40t1n_36p/analysis/increment/xainc.20200110_030000z.nc4'
frtdir = '/work2/noaa/gsienkf/weihuang/jedi/develop/build/fv3-jedi/test/Data/analysis/letkf/gfs/sts'
snddir = '/work2/noaa/gsienkf/weihuang/jedi/develop/build/fv3-jedi/test/Data/analysis/letkf/gfs/mts'

frtfile = '%s/xainc.20201215_000000z.nc4' %(frtdir)
sndfile = '%s/xainc.20201215_000000z.nc4' %(snddir)

opts, args = getopt.getopt(sys.argv[1:], '', ['debug=', 'output=',
'jedifile=', 'gsifile='])
'sndfile=', 'frtfile='])
for o, a in opts:
if o in ('--debug'):
debug = int(a)
elif o in ('--output'):
output = int(a)
elif o in ('--jedifile'):
jedifile = a
elif o in ('--gsifile'):
gsifile = a
elif o in ('--sndfile'):
sndfile = a
elif o in ('--frtfile'):
frtfile = a
else:
assert False, 'unhandled option'

gp = GeneratePlot(debug=debug, output=output)

ncjedi = netcdf_dataset(jedifile)
ncgsi = netcdf_dataset(gsifile)
lats = ncjedi.variables['lat'][:]
lons = ncjedi.variables['lon'][:]
ncfrt = netcdf_dataset(frtfile)
ncsnd = netcdf_dataset(sndfile)

lats = ncfrt.variables['lat'][:]
lons = ncfrt.variables['lon'][:]

#-----------------------------------------------------------------------------------------
clevs = np.arange(-1.0, 1.01, 0.01)
Expand All @@ -171,33 +174,34 @@ def set_cmapname(self, cmapname):
gp.set_cblevs(cblevs=cblevs)

#-----------------------------------------------------------------------------------------
#jedi_varlist = ['T', 'ua', 'va', 'sphum', 'delp', 'DZ', 'o3mr']
#gsi_varlist = ['T', 'ua', 'va', 'sphum', 'delp', 'DZ', 'o3mr']
jedi_varlist = ['T', 'delp', 'sphum']
gsi_varlist = ['T', 'delp', 'sphum']

#snd_varlist = ['T', 'ua', 'va', 'sphum', 'DELP', 'DZ', 'o3mr']
#frt_varlist = ['T', 'ua', 'va', 'sphum', 'DELP', 'DZ', 'o3mr']
#unitlist = ['Unit (C)', 'Unit (m/s)', 'Unit (m/s)',
# 'Unit (kg/kg)', 'Unit (Pa)', 'Unit (m)', 'Unit (ppm)']
snd_varlist = ['T', 'ua', 'va', 'sphum', 'DELP', 'o3mr']
frt_varlist = ['T', 'ua', 'va', 'sphum', 'DELP', 'o3mr']
unitlist = ['Unit (C)', 'Unit (m/s)', 'Unit (m/s)',
'Unit (kg/kg)', 'Unit (Pa', 'Unit (m', 'Unit (ppm)']
'Unit (kg/kg)', 'Unit (Pa)', 'Unit (ppm)']

#-----------------------------------------------------------------------------------------
for n in range(len(jedi_varlist)):
jedivar = ncjedi.variables[jedi_varlist[n]][0, :, :, :]
gsivar = ncgsi.variables[gsi_varlist[n]][0,:, :, :]
for n in range(len(snd_varlist)):
sndvar = ncsnd.variables[snd_varlist[n]][0, :, :, :]
frtvar = ncfrt.variables[frt_varlist[n]][0,:, :, :]

nlev, nlat, nlon = jedivar.shape
print('jedivar.shape = ', jedivar.shape)
print('gsivar.shape = ', gsivar.shape)
nlev, nlat, nlon = sndvar.shape
print('sndvar.shape = ', sndvar.shape)
print('frtvar.shape = ', frtvar.shape)

gp.set_label(unitlist[n])

for lev in range(5, nlev, 10):
v0 = gsivar[lev,:,:]
v1 = jedivar[lev,:,:]
v0 = frtvar[lev,:,:]
v1 = sndvar[lev,:,:]
v2 = v1 - v0

data = [v0, v1, v2]

title = '%s at Level %d' %(jedi_varlist[n], lev)
title = '%s at Level %d' %(snd_varlist[n], lev)
gp.set_title(title)

print('Plotting ', title)
Expand All @@ -209,12 +213,12 @@ def set_cmapname(self, cmapname):
print('\tv1.max: %f, v1.min: %f' %(np.max(v1), np.min(v1)))
print('\tv2.max: %f, v2.min: %f' %(np.max(v2), np.min(v2)))

imagename = '%s_lev_%3.3d.png' %(jedi_varlist[n], lev)
imagename = '%s_lev_%3.3d.png' %(snd_varlist[n], lev)
gp.set_imagename(imagename)

gp.plot(lons, lats, data=data)

#-----------------------------------------------------------------------------------------
ncjedi.close()
ncgsi.close()
ncsnd.close()
ncfrt.close()

200 changes: 200 additions & 0 deletions plot-mean.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,200 @@
import getopt
import os, sys
import numpy as np

import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from mpl_toolkits.axes_grid1 import make_axes_locatable

import cartopy.crs as ccrs
from cartopy import config
from cartopy.util import add_cyclic_point
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker

from netCDF4 import Dataset as netcdf_dataset

#=========================================================================
class GeneratePlot():
def __init__(self, debug=0, output=0):
self.debug = debug
self.output = output

self.set_default()

def plot(self, lons, lats, data):
#ax.coastlines(resolution='110m')
#ax.gridlines()

nrows = 1
ncols = 1

#set up the plot
proj = ccrs.PlateCarree()

fig, axs = plt.subplots(nrows=nrows,ncols=ncols,
subplot_kw=dict(projection=proj),
figsize=(11,8.5))
axs.set_global()
pvar = data

cyclic_data, cyclic_lons = add_cyclic_point(pvar, coord=lons)

cs=axs.contourf(cyclic_lons, lats, cyclic_data, transform=proj,
levels=self.clevs, extend=self.extend,
alpha=self.alpha, cmap=self.cmapname)

axs.set_extent([-180, 180, -90, 90], crs=proj)
axs.coastlines(resolution='auto', color='k')
axs.gridlines(color='lightgrey', linestyle='-', draw_labels=True)
axs.set_title(self.runname[0])

#Adjust the location of the subplots on the page to make room for the colorbar
fig.subplots_adjust(bottom=0.1, top=0.9, left=0.05, right=0.8,
wspace=0.02, hspace=0.02)

#Add a colorbar axis at the bottom of the graph
cbar_ax = fig.add_axes([0.85, 0.1, 0.05, 0.85])

#Draw the colorbar
cbar=fig.colorbar(cs, cax=cbar_ax, pad=self.pad, ticks=self.cblevs,
orientation='vertical')

cbar.set_label(self.label, rotation=90)

#Add a big title at the top
plt.suptitle(self.title)

fig.canvas.draw()
plt.tight_layout()

if(self.output):
if(self.imagename is None):
imagename = 't_aspect.png'
else:
imagename = self.imagename
plt.savefig(imagename)
plt.close()
else:
plt.show()

def set_default(self):
self.imagename = 'sample.png'

self.runname = ['Halo', 'RR', 'RR - Halo']

#cmapname = coolwarm, bwr, rainbow, jet, seismic
#self.cmapname = 'bwr'
#self.cmapname = 'coolwarm'
self.cmapname = 'rainbow'
#self.cmapname = 'jet'

self.clevs = np.arange(-0.2, 0.21, 0.01)
self.cblevs = np.arange(-0.2, 0.3, 0.1)

self.extend = 'both'
self.alpha = 0.5
self.pad = 0.1
self.orientation = 'horizontal'
self.size = 'large'
self.weight = 'bold'
self.labelsize = 'medium'

self.label = 'Unit (C)'
self.title = 'Temperature Increment'

def set_label(self, label='Unit (C)'):
self.label = label

def set_title(self, title='Temperature Increment'):
self.title = title

def set_clevs(self, clevs=[]):
self.clevs = clevs

def set_cblevs(self, cblevs=[]):
self.cblevs = cblevs

def set_imagename(self, imagename):
self.imagename = imagename

def set_cmapname(self, cmapname):
self.cmapname = cmapname

def set_runname(self, runname):
self.runname = runname

#--------------------------------------------------------------------------------
if __name__== '__main__':
debug = 1
output = 0
#datadir = '/work2/noaa/gsienkf/weihuang/jedi/develop/build/fv3-jedi/test/Data/analysis/letkf/gfs/mts'
datadir = '/work2/noaa/gsienkf/weihuang/jedi/develop/build/fv3-jedi/test'

filename = 'letkf.mid.00020201215_000000z.nc4'

opts, args = getopt.getopt(sys.argv[1:], '', ['debug=', 'output=',
'filename=', 'datadir='])
for o, a in opts:
if o in ('--debug'):
debug = int(a)
elif o in ('--output'):
output = int(a)
elif o in ('--datadir'):
datadir = a
elif o in ('--filename'):
filename = a
else:
assert False, 'unhandled option'

gp = GeneratePlot(debug=debug, output=output)

fullpath = '%s/%s' %(datadir, filename)
ncfile = netcdf_dataset(fullpath)

lats = ncfile.variables['lat'][:]
lons = ncfile.variables['lon'][:]

#-----------------------------------------------------------------------------------------
clevs = np.arange(200.0, 306.0, 1.0)
cblevs = np.arange(200.0, 310.0, 5.0)

gp.set_clevs(clevs=clevs)
gp.set_cblevs(cblevs=cblevs)

#-----------------------------------------------------------------------------------------
varlist = ['T', 'ua', 'va', 'sphum', 'delp', 'DZ', 'o3mr']
unitlist = ['Unit (C)', 'Unit (m/s)', 'Unit (m/s)',
'Unit (kg/kg)', 'Unit (Pa', 'Unit (m', 'Unit (ppm)']

#-----------------------------------------------------------------------------------------
for n in range(len(varlist)):
var = ncfile.variables[varlist[n]][0,:, :, :]
gp.set_runname([varlist[n]])

nlev, nlat, nlon = var.shape
print('var.shape = ', var.shape)

gp.set_label(unitlist[n])

for lev in range(5, nlev, 10):
data = var[lev,:,:]

title = '%s at Level %d' %(varlist[n], lev)
gp.set_title(title)

print('Plotting ', title)
print('\tdata.shape = ', data.shape)

print('\tdata.max: %f, data.min: %f' %(np.max(data), np.min(data)))

imagename = '%s_lev_%3.3d.png' %(varlist[n], lev)
gp.set_imagename(imagename)

gp.plot(lons, lats, data)

#-----------------------------------------------------------------------------------------
ncfile.close()

Loading

0 comments on commit cc9af07

Please sign in to comment.