plot tile distance

.gitignore

@@ -5,3 +5,4 @@ regrid/
 *.txt
 *.csv
 timingit
+log.*

all-obs-jedi.py

@@ -0,0 +1,164 @@
+#=========================================================================
+import os
+import sys
+import types
+import getopt
+import netCDF4
+import matplotlib
+
+import numpy as np
+import matplotlib.pyplot
+
+from matplotlib import cm
+from mpl_toolkits.basemap import Basemap
+
+from genplot import GeneratePlot as genplot
+from readIODA2Obs import ReadIODA2Obs
+
+#=========================================================================
+class PlotObsOverlayIncrement():
+  def __init__(self, debug=0, output=0, filename=None):
+    self.debug = debug
+    self.output = output
+    self.filename = filename
+
+    if(self.debug):
+      print('debug = ', debug)
+
+    if(self.debug > 10):
+      print('self.filename = ', self.filename)
+
+  def get_var(self, varname):
+   #print('varname =', varname)
+
+    ncfile = netCDF4.Dataset(self.filename, 'r')
+    lat = ncfile.variables['lat'][:]
+    lon = ncfile.variables['lon'][:]
+   #var = ncfile.variables[varname][itime, :, :, :]
+    var = ncfile.variables[varname][:, :, :]
+    ncfile.close()
+
+    var[np.isnan(var)] = 0.0
+
+    if(self.debug):
+      msg = ('var range for variable %s: (%s, %s).' % (varname, var.min(), var.max()))
+      print(msg)
+
+    return lon, lat, var
+
+#------------------------------------------------------------------------------
+if __name__ == '__main__':
+  debug = 1
+  output = 0
+  uselogp = 1
+  varname = 'T'
+  casename = 'surf'
+  gridfile = '/work2/noaa/gsienkf/weihuang/jedi/case_study/vis/regrid/fv3latlon.nc'
+
+  obsdir = '/work2/noaa/gsienkf/weihuang/jedi/case_study/ioda_v2_data'
+  filelist = ['aircraft_tsen_obs_2020011006.nc4', 'amsua_n19_obs_2020121500_m.nc4',
+              'iasi_metop-b_obs_2020121500_m.nc4', 'satwind_obs_2020011006.nc4',
+              'scatwind_obs_2020011006.nc4', 'sfc_ps_obs_2020011006.nc4',
+              'sfcship_tsen_obs_2020011006.nc4', 'sondes_tsen_obs_2020011006.nc4',
+              'vadwind_obs_2020011006.nc4', 'windprof_obs_2020011006.nc4']
+
+  opts, args = getopt.getopt(sys.argv[1:], '', ['debug=', 'output=', 'uselogp=',
+                             'gridfile=', 'obsfile=', 'casename=', 'varname='])
+
+  for o, a in opts:
+    if o in ('--debug'):
+      debug = int(a)
+    elif o in ('--output'):
+      output = int(a)
+    elif o in ('--uselogp'):
+      uselogp = int(a)
+    elif o in ('--casename'):
+      casename = a
+    elif o in ('--varname'):
+      varname = a
+    elif o in ('--gridfile'):
+      gridfile = a
+    elif o in ('--obsfile'):
+      obsfile = a
+   #else:
+   #  assert False, 'unhandled option'
+
+  print('debug = ', debug)
+  print('output = ', output)
+
+#------------------------------------------------------------------------------
+  po = PlotObsOverlayIncrement(debug=debug, output=output, filename=gridfile)
+  lon1d, lat1d, var = po.get_var(varname)
+
+#------------------------------------------------------------------------------
+  gp = genplot(debug=debug, output=output, lat=lat1d, lon=lon1d)
+  robs = ReadIODA2Obs(debug=debug, filename=None)
+  obslist = []
+  for file in filelist:
+    filename = '%s/%s' %(obsdir, file)
+    obslat, obslon = robs.get_latlon_from_file(filename)
+    obs = [obslon, obslat]
+    obslist.append(obs)
+
+#------------------------------------------------------------------------------
+  gp.set_label('Temperature (K)')
+
+  imageprefix = '%s_%s_increment' %(casename, varname)
+  titleprefix = '%s %s Increment at' %(casename, varname)
+
+#------------------------------------------------------------------------------
+ #clevs = np.arange(-0.5, 0.51, 0.01)
+ #cblevs = np.arange(-0.5, 0.6, 0.1)
+ #clevs = np.arange(-0.2, 0.21, 0.01)
+ #cblevs = np.arange(-0.2, 0.3, 0.1)
+ #clevs = np.arange(-2., 2.1, 0.1)
+ #cblevs = np.arange(-2., 3., 1.)
+  clevs = np.arange(-5.0, 5.1, 0.1)
+  cblevs = np.arange(-5.0, 6.0, 1.0)
+  gp.set_clevs(clevs=clevs)
+  gp.set_cblevs(cblevs=cblevs)
+
+  levs = [30, 40, 50, 60]
+
+  for lev in levs:
+    pvar = var[lev,:,:]
+    imgname = '%s_lev_%d.png' %(imageprefix, lev)
+    title = '%s level %d' %(titleprefix, lev)
+    gp.set_imagename(imgname)
+    gp.set_title(title)
+    gp.plot_with_obs(pvar, obslist)
+
+  sys.exit(-1)
+
+#------------------------------------------------------------------------------
+  lons = [40, 105, 170, 270, 300]
+
+  for lon in lons:
+    pvar = var[:,:,lon]
+    title = '%s longitude %d' %(titleprefix, lon)
+    gp.set_title(title)
+
+    imgname = '%s_lon_%d_logp.png' %(imageprefix, lon)
+    gp.set_imagename(imgname)
+    gp.plot_meridional_section_logp(pvar)
+
+    imgname = '%s_lon_%d_level.png' %(imageprefix, lon)
+    gp.set_imagename(imgname)
+    gp.plot_meridional_section(pvar)
+
+#------------------------------------------------------------------------------
+  lats = [-30, 0, 45, 70]
+
+  for lat in lats:
+    pvar = var[:,90+lat,:]
+    gp.set_title(title)
+    title = '%s latitude %d' %(titleprefix, lat)
+
+    imgname = '%s_lat_%d_logp.png' %(imageprefix, lat)
+    gp.set_imagename(imgname)
+    gp.plot_zonal_section_logp(pvar)
+
+    imgname = '%s_lat_%d_level.png' %(imageprefix, lat)
+    gp.set_imagename(imgname)
+    gp.plot_zonal_section(pvar)
+

genplot.py

@@ -586,6 +586,78 @@ def plot_coast_lat_lon_line(self):
                                color=color, linewidth=linewidth,
                                dashes=dashes, fontsize=fontsize)
 
+  def plot_with_obs(self, pvar, obslist):
+    self.basemap = self.build_basemap()
+
+    self.plt = matplotlib.pyplot
+    try:
+      self.plt.close('all')
+      self.plt.clf()
+    except Exception:
+      pass
+
+    self.fig = self.plt.figure()
+    self.ax = self.plt.subplot()
+
+    msg = ('plot variable min: %s, max: %s' % (np.min(pvar), np.max(pvar)))
+    print(msg)
+    print('pvar.shape = ', pvar.shape)
+    print('pvar.size = ', pvar.size)
+
+    (self.x, self.y) = self.basemap(self.lon1d, self.lat1d)
+    v1d = np.reshape(pvar, (pvar.size, ))
+
+    contfill = self.basemap.contourf(self.x, self.y, v1d, tri=True,
+                                     levels=self.clevs, extend=self.extend,
+                                     alpha=self.alpha, cmap=self.cmapname)
+
+    cb = self.fig.colorbar(contfill, orientation=self.orientation,
+                           pad=self.pad, ticks=self.cblevs)
+
+    cb.set_label(label=self.label, size=self.size, weight=self.weight)
+
+    cb.ax.tick_params(labelsize=self.labelsize)
+    if(self.precision == 0):
+      cb.ax.set_xticklabels(['{:.0f}'.format(x) for x in self.cblevs], minor=False)
+    elif(self.precision == 1):
+      cb.ax.set_xticklabels(['{:.1f}'.format(x) for x in self.cblevs], minor=False)
+    elif(self.precision == 2):
+      cb.ax.set_xticklabels(['{:.2f}'.format(x) for x in self.cblevs], minor=False)
+    else:
+      cb.ax.set_xticklabels(['{:.3f}'.format(x) for x in self.cblevs], minor=False)
+
+    self.ax.set_title(self.title)
+
+    self.plot_coast_lat_lon_line()
+
+    self.obscolorlist = ['greenyellow', 'lawngreen', 'darkseagreen', 'palegreen',
+                         'forestgreen', 'limegreen', 'darkgreen', 'green',
+                         'lime', 'seagreen', 'springgreen', 'mediumseagreen',
+                         'lightgreen', 'turquoise', 'darkcyan', 'yellowgreen',
+                         'aquamarine', 'teal', 'cyan', 'aqua', 'deepskyblue']
+   #https://matplotlib.org/stable/api/markers_api.html
+    self.obsmarklist = ['o', 'v', '^', '<', '>', '1', '2', '3', '4',
+                        '8', 's', 'p', 'P', '*', 'h', 'H', '+', 'x', 'X',
+                        'D', 'd', '|', '-']
+
+    nm = 0
+    markersize=1
+    for obs in obslist:
+      obslon = obs[0]
+      obslat = obs[1]
+      self.add_obs_marker_at_latlon(obs[0], obs[1], marker=self.obsmarklist[nm],
+                                    size=markersize, color=self.obscolorlist[nm])
+      nm += 1
+
+    self.display(output=self.output, image_name=self.image_name)
+
+  def add_obs_marker_at_latlon(self, obslon, obslat, marker='x', size=3, color='green'):
+    if(len(obslon) > 0):
+      x, y = self.basemap(obslon, obslat)
+     #adding dotes:
+     #dotes = self.basemap.plot(x, y, 'bo', markersize=12)
+     #dotes = self.basemap.plot(x, y, 'bo', markersize=6)
+      dotes = self.basemap.scatter(x, y, marker=marker, s=size, color=color)
 # ----
 if __name__ == '__main__':
   debug = 1

plot_fv3_increment.sh

@@ -0,0 +1,23 @@
+#!/bin/bash
+
+ datadir=/work2/noaa/gsienkf/weihuang/jedi/run/run_80.40t4n_156p/analysis/increment
+ obsfile=/work2/noaa/gsienkf/weihuang/jedi/case_study/surf/ioda_v2_data/sfc_ps_obs_2020011006.nc4
+
+ output=0
+
+ gridfile=regrid/fv3latlon.nc
+ if [ ! -f ${gridfile} ]
+ then
+   cd regrid
+  #interp.sh ${case_dir}/sfc-letkf/analysis/increment/
+   interp.sh ${datadir}/
+   cd ..
+ fi
+
+ python all-obs-jedi.py \
+   --output=${output} \
+   --casename=AllObs \
+   --varname=T \
+   --gridfile=${gridfile} \
+   --obsfile=${obsfile}
+

plot_tile_orograph.py

@@ -210,26 +210,27 @@ def __init__(self, debug=0, output=0, griddir=None, gridtype='C48'):
     self.griddir = griddir
     self.gridtype = gridtype
 
-    if(griddi is None):
+    if(griddir is None):
       print('griddir not defined. Exit.')
       sys.exit(-1)
 
+    self.datalist = []
     self.lon = []
     self.lat = []
     for n in range(6):
       nt = n + 1
       datafile = '%s/%s/%s_oro_data.tile%d.nc' %(griddir, gridtype, gridtype, nt)
       if(os.path.exists(datafile)):
-       #print('File No %d: %s' %(nt, datafile))
+        print('File No %d: %s' %(nt, datafile))
         self.datalist.append(datafile)
 
       ncf = netCDF4.Dataset(datafile)
       lon = ncf.variables['geolon'][:,:]
       lat = ncf.variables['geolat'][:,:]
 
-     #print('lon.ndim=', lon.ndim)
-     #print('lon.shape=', lon.shape)
-     #print('lon.size=', lon.size)
+      print('lon.ndim=', lon.ndim)
+      print('lon.shape=', lon.shape)
+      print('lon.size=', lon.size)
 
       ny, nx = lon.shape
 
@@ -244,32 +245,23 @@ def __init__(self, debug=0, output=0, griddir=None, gridtype='C48'):
     print('len(self.lon) = ', len(self.lon))
     print('len(self.lat) = ', len(self.lat))
 
+  def get_latlon(self):
+    return np.array(self.lat), np.array(self.lon)
+
   def get_var(self, varname):
     data = []
     for datafile in self.datalist:
       ncf = netCDF4.Dataset(datafile)
       var = ncf.variables[varname][:,:]
 
-     #print('lon.ndim=', lon.ndim)
-     #print('lon.shape=', lon.shape)
-     #print('lon.size=', lon.size)
-
-      ny, nx = lon.shape
-
-      lonc1d = np.reshape(lon, (nx*ny,))
-      latc1d = np.reshape(lat, (nx*ny,))
-      var1d = np.reshape(lat, (nx*ny,))
-
-      self.lon.extend(lonc1d)
-      self.lat.extend(latc1d)
-      self.var.extend(var1d)
+      ny, nx = var.shape
+      var1d = np.reshape(var, (nx*ny,))
+      data.extend(var1d)
       ncf.close()
 
-    print('len(self.lon) = ', len(self.lon))
-    print('len(self.lat) = ', len(self.lat))
-    print('len(self.orog) = ', len(self.orog))
+    print('len(data) = ', len(data))
 
-    return np.array(self.lat), np.array(self.lon), np.array(self.orog)
+    return np.array(data)
 
 #--------------------------------------------------------------------------------
 if __name__== '__main__':
@@ -292,20 +284,14 @@ def get_var(self, varname):
     else:
       assert False, 'unhandled option'
 
-  pt = TileData(debug=debug, griddir=griddir, gridtype=gridtype)
-  pt.process()
-  lat, lon, orog = dt.get_data()
+  td = TileData(debug=debug, griddir=griddir, gridtype=gridtype)
+  lat, lon = td.get_latlon()
+  orog = td.get_var('orog_filt')
 
   gp = GeneratePlot(debug=debug, output=output)
   gp.set_grid(lat, lon)
 
-  imgname = 'ps.png'
-  gp.set_imagename(imgname)
-  title = 'Surface Pressure (Unit: hPa)'
-  gp.set_title(title)
-  gp.plot(ps)
-
-  imgname = 'orog.png'
+  imgname = 'direct_tile_orog.png'
   gp.set_imagename(imgname)
   title = 'Orograph (Unit: m)'
   gp.set_title(title)

readIODA2Obs.py

@@ -21,7 +21,8 @@ def __init__(self, debug=0, filename=None):
     self.nstring = 0
     self.nvars = 0
 
-    self.set_vardims()
+    if(filename != None):
+      self.set_vardims()
 
   def set_filename(self, filename=None):
     self.filename = filename
@@ -160,6 +161,17 @@ def get_latlon(self):
 
     return lat, lon
 
+  def get_latlon_from_file(self, filename):
+    ncfile = netCDF4.Dataset(filename, 'r')
+    ncgroup = ncfile['MetaData']
+    lat = ncgroup.variables['latitude'][:]
+    lon = ncgroup.variables['longitude'][:]
+    ncfile.close()
+
+    lon = np.where(lon > 0, lon, lon+360.0)
+
+    return lat, lon
+
   def get_latlon4var(self, varname=None):
     lat, lon = self.get_latlon()
     if(varname is None):