soca test

convertGSIhofx2ioda/createGSIobs.py

@@ -173,10 +173,10 @@ def write_var(self, ncout, variable, varname):
 #-----------------------------------------------------------------------------------------
 if __name__ == '__main__':
   debug = 1
-  dirname = '/work2/noaa/gsienkf/weihuang/jedi/case_study/surf/ioda_v2_data'
- #filename = 'sfc_ps_obs_2020011006_0000.nc4'
- #filename = 'sfcship_ps_obs_2020011006_0000.nc4'
-  filename = 'sondes_ps_obs_2020011006_0000.nc4'
+  indir = '/work2/noaa/gsienkf/weihuang/jedi/case_study/sfc-letkf/obsout'
+  infile = 'sfcship_ps_obs_2020011006_0000.nc4'
+  outdir = '/work2/noaa/gsienkf/weihuang/jedi/case_study/sfc-letkf/gsihofxbc'
+  outfile = 'sfcship_ps_obs_2020011006_0000.nc4'
 
   opts, args = getopt.getopt(sys.argv[1:], '', ['debug=', 'dirname=', 'filename='])
 

convertGSIhofx2ioda/mf.py

@@ -1,45 +1,170 @@
 #=========================================================================
+import os
 import sys
+import getopt
+import math
+import numpy as np
 import netCDF4 as nc4
 
-#-----------------------------------------------------------------------------------------
-debug = 1
-dirname = 'sfc_ps_out'
-#filename = 'sfc_ps_obs_2020011006_0000.nc4'
-#filename = 'sfcship_ps_obs_2020011006_0000.nc4'
-filename = 'sondes_ps_obs_2020011006_0000.nc4'
+#=========================================================================
+class CreateGSIobs():
+  def __init__(self, debug=0, infile=None, outfile=None):
+    self.debug = debug
+    self.infile = infile
+    self.outfile = outfile
+
+    self.ncin = nc4.Dataset(self.infile, 'r')
+    self.ncout = nc4.Dataset(self.outfile, 'w')
+
+ #Destructor
+  def __del__(self):
+   #Close the netcdf files
+    self.ncin.close()
+    self.ncout.close()
+
+  def process(self, varlist=[]):
+   #NetCDF global attributes
+    self.attrs = self.ncin.ncattrs()
+    if (self.debug):
+      print("NetCDF Global Attributes:")
+    for attr in self.attrs:
+      if (self.debug):
+        print('\t%s:' % attr, repr(self.ncin.getncattr(attr)))
+      self.ncout.setncattr(attr, self.ncin.getncattr(attr))
 
-fullpath = '%s/%s' %(dirname, filename)
+   #Dimension shape information.
+    if (self.debug):
+      print("NetCDF dimension information:")
 
-print('debug = ', debug)
-print('fullpath = ', fullpath)
+    for dim in self.ncin.dimensions:
+      if (self.debug):
+        print("\tName:", dim)
+        print("\t\tsize:", len(self.ncin.dimensions[dim]))
 
-#nco = nc4.Dataset(fullpath, 'r+')
-nco = nc4.Dataset('sfc_ps_out/sondes_ps_obs_2020011006_0000.nc', 'r+')
-#nco = nc4.Dataset('nc1.nc', 'r+')
+      dimval = self.ncin.dimensions[dim]
+      if dimval.isunlimited():
+        self.ncout.createDimension(dim, None)
+      else:
+        self.ncout.createDimension(dim, len(dimval))
 
-gsiflnm = filename.replace('_0000.nc4', '.nc4')
+   #Variable information.
+    if (self.debug):
+      print("NetCDF variable information:")
+    i = 0
+    for varname in self.ncin.variables:
+      i += 1
+      if (self.debug):
+        print('\t\tvar No %d Name: %s' %(i, varname))
+      variable = self.ncin.variables[varname]
+      self.write_var(self.ncout, variable, varname)
 
-#for varname in ['air_temperature', 'eastward_wind', 'northward_wind',
-#              'specific_humidity', 'virtual_temperature']:
-for varname in ['surface_pressure']:
-  print('\tvarname = ', varname)
+   #Group information.
+    if (self.debug):
+      print("NetCDF group information:")
 
-  for n in range(81):
-    if(0 == n):
-      grp = 'hofx_y_mean_xb0'
-      gsipath = 'ensmean/%s' %(gsiflnm)
+    n = 0
+    for grp in self.ncin.groups:
+      n += 1
+      if (self.debug):
+        print('\tgroup No. %d: %s' %(n, grp))
+
+      grpout = self.ncout.createGroup(grp)
+
+      if (self.debug):
+        print("NetCDF variable information:")
+      i = 0
+      for varname in self.ncin.groups[grp].variables:
+        i += 1
+        if (self.debug):
+          print('\t\tvar No %d Name: %s' %(i, varname))
+        variable = self.ncin.groups[grp].variables[varname]
+        if(grp == 'ombg'):
+          hofxbc = self.ncin.groups['GsiHofXBc'].variables[varname][:]
+          obsval = self.ncin.groups['ObsValue'].variables[varname][:]
+          value = obsval - hofxbc
+          self.write_var_with_value(grpout, variable, varname, value)
+        else:
+          self.write_var(grpout, variable, varname)
+
+  def write_var_with_value(self, ncout, variable, varname, value):
+    attrs = variable.ncattrs()
+    if('_FillValue' in attrs):
+      fv = variable.getncattr('_FillValue')
+      ncout.createVariable(varname, variable.datatype, variable.dimensions,
+                           fill_value=fv)
     else:
-      grp = 'hofx0_%d' %(n)
-      gsipath = 'mem%3.3d/%s' %(n, gsiflnm)
+      ncout.createVariable(varname, variable.datatype, variable.dimensions)
+
+   #NetCDF variable attributes
+    for attr in attrs:
+      if('_FillValue' == attr):
+        continue
+      attr_value = variable.getncattr(attr)
+      ncout.variables[varname].setncattr(attr, attr_value)
+
+    ncout.variables[varname][:] = value
+
+  def write_var(self, ncout, variable, varname):
+    newVarname = varname
+   #if(varname == 'dateTime'):
+   #  newVarname = 'datetime'
+   #else:
+   #  newVarname = varname
+
+    if (self.debug):
+      print('variable.datatype = ', variable.datatype)
+      print('variable.dimensions = ', variable.dimensions)
+
+    attrs = variable.ncattrs()
+    if('_FillValue' in attrs):
+      fv = variable.getncattr('_FillValue')
+      ncout.createVariable(newVarname, variable.datatype, variable.dimensions,
+                                fill_value=fv)
+    else:
+      ncout.createVariable(newVarname, variable.datatype, variable.dimensions)
+
+   #NetCDF variable attributes
+    if (self.debug):
+      print("NetCDF variable Attributes:")
+    for attr in attrs:
+      if('_FillValue' == attr):
+        continue
+      attr_value = variable.getncattr(attr)
+      if (self.debug):
+        print('\t%s:' % attr, attr_value)
+
+      ncout.variables[newVarname].setncattr(attr, attr_value)
+
+    var = variable[:]
+    ncout.variables[newVarname][:] = var
+
+#-----------------------------------------------------------------------------------------
+if __name__ == '__main__':
+  debug = 1
+ #indir = '/work2/noaa/gsienkf/weihuang/jedi/case_study/sfc-letkf/observer'
+ #outdir = '/work2/noaa/gsienkf/weihuang/jedi/case_study/sfc-letkf/gsihofxbc'
+ #infile = '%s/sfcship_ps_obs_2020011006_0000.nc4' %(indir)
+ #outfile = '%s/sfcship_ps_obs_2020011006_0000.nc4' %(outdir)
+  indir = '/work2/noaa/gsienkf/weihuang/jedi/case_study/sfc-letkf/obsout.after_observer.sfc_ps'
+  outdir = '/work2/noaa/gsienkf/weihuang/jedi/case_study/sfc-letkf/gsihofxbc.sfc_ps'
+  infile = '%s/sfc_ps_obs_2020011006_0000.nc4' %(indir)
+  outfile = '%s/sfc_ps_obs_2020011006_0000.nc4' %(outdir)
 
-    nci = nc4.Dataset(gsipath, 'r')
-    var = nci.groups['GsiHofX'].variables[varname][:]
-    nci.close()
+  opts, args = getopt.getopt(sys.argv[1:], '', ['debug=', 'infile=', 'outfile='])
 
-    print('Group No %d: %s, gsipath: %s' %(n, grp, gsipath))
+  for o, a in opts:
+    if o in ('--debug'):
+      debug = int(a)
+    elif o in ('--indir'):
+      dirname = a
+    elif o in ('--infile'):
+      infile = a
+    elif o in ('--outfile'):
+      outfile = a
+   #else:
+   #  assert False, 'unhandled option'
 
-    nco.groups[grp].variables[varname][:] = var
+  cg = CreateGSIobs(debug=debug, infile=infile, outfile=outfile)
 
-nco.close()
+  cg.process(varlist=['surface_pressure'])
 

main_timing.py

@@ -190,8 +190,10 @@ def plot(self):
     else:
      #plt.xscale('log', base=2)
      #plt.yscale('log', base=10)
-      plt.xscale('log', basex=2)
-      plt.yscale('log', basey=2)
+     #plt.xscale('log', basex=2)
+     #plt.yscale('log', basey=2)
+      plt.xscale('log', base=2)
+      plt.yscale('log', base=2)
       plt.xticks(x, xlabels)
      #plt.xticks(x, xlabels, rotation ='vertical')
       plt.yticks(yp, ylabels)
@@ -306,8 +308,8 @@ def plot(self):
 if __name__== '__main__':
   debug = 1
   casename = 'sondes'
- #workdir = '/work2/noaa/gsienkf/weihuang/jedi/case_study'
-  workdir = '/scratch2/BMC/gsienkf/Wei.Huang/tools/vis_tools'
+  workdir = '/work2/noaa/gsienkf/weihuang/jedi/case_study'
+ #workdir = '/scratch2/BMC/gsienkf/Wei.Huang/tools/vis_tools'
   corelist = [36, 78, 156, 312]
  #corelist = [36, 72, 144, 288]
   nodelist = [1, 2, 4, 8]

obs_timing.py

@@ -183,8 +183,10 @@ def plot(self):
     else:
      #plt.xscale('log', base=2)
      #plt.yscale('log', base=10)
-      plt.xscale('log', basex=2)
-      plt.yscale('log', basey=2)
+     #plt.xscale('log', basex=2)
+     #plt.yscale('log', basey=2)
+      plt.xscale('log', base=2)
+      plt.yscale('log', base=2)
       plt.xticks(x, xlabels)
      #plt.xticks(x, xlabels, rotation ='vertical')
       plt.yticks(yp, ylabels)
@@ -299,15 +301,15 @@ def plot(self):
 if __name__== '__main__':
   debug = 1
   casename = 'sondes'
- #workdir = '/work2/noaa/gsienkf/weihuang/jedi/case_study'
-  workdir = '/scratch2/BMC/gsienkf/Wei.Huang/tools/vis_tools'
+  workdir = '/work2/noaa/gsienkf/weihuang/jedi/case_study'
+ #workdir = '/scratch2/BMC/gsienkf/Wei.Huang/tools/vis_tools'
   corelist = [36, 78, 156, 312]
  #corelist = [36, 72, 144, 288]
   nodelist = [1, 2, 4, 8]
   output = 0
   linear = 1
 
-  opts, args = getopt.getopt(sys.argv[1:], '', ['debug=', 'workdir=',
+  opts, args = getopt.getopt(sys.argv[1:], '', ['debug=', 'workdir=', 'output=',
                              'corelist=', 'nodelist=', 'casename='])
 
   for o, a in opts:
@@ -322,7 +324,7 @@ def plot(self):
     elif o in ('--casename'):
       casename = a
     elif o in ('--output'):
-      output = a
+      output = int(a)
     elif o in ('--linear'):
       linear = int(a)
     else:
@@ -331,11 +333,11 @@ def plot(self):
   pr = Profiler(debug=debug, corelist=corelist, nodelist=nodelist, output=output,
                 workdir=workdir, casename=casename, linear=linear)
   pr.process()
- #for linear in [1, 0]:
-  for linear in [0]:
+  for linear in [1, 0]:
+ #for linear in [0]:
     pr.set_linear(linear=linear)
-    for output in [0, 1]:
-   #for output in [1]:
+   #for output in [0, 1]:
+    for output in [1]:
       pr.set_output(output=output)
       pr.plot()