Skip to content

Commit

Permalink
add
Browse files Browse the repository at this point in the history
  • Loading branch information
@weihuang-jedi
weihuang-jedi committed Feb 10, 2024
1 parent 1cf3166 commit 5608f73
Show file tree
Hide file tree
Showing 7 changed files with 1,435 additions and 113 deletions.
223 changes: 223 additions & 0 deletions diff2obs.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,223 @@
#########################################################################
#$Id: bld.py 28 2021-01-21 15:10:31Z whuang $
#$Revision: 28 $
#$HeadURL: file:///Users/whuang/.wei_svn_repository/trunk/jedi-build-tools/bld.py $
#$Date: 2021-01-21 08:10:31 -0700 (Thu, 21 Jan 2021) $
#$Author: whuang $
#########################################################################

import getopt
import os, sys
import types
import time
import datetime
import subprocess
import netCDF4

import numpy as np


from readIODA2Obs import ReadIODA2Obs

#------------------------------------------------------------------
class Compare2ObsFiles():
def __init__(self, debug=0):
self.debug = debug

self.delt = 1.0e-10

#------------------------------------------------------------------
def setup(self, datadir=None, filename=None):
if(self.debug):
print('datadir = ', datadir)
print('filename = ', filename)

basefile = '%s/hofx.nl/%s' %(datadir, filename)
if(os.path.exists(basefile)):
self.baseio = ReadIODA2Obs(debug=debug, filename=basefile)
else:
print('basefile: <%s> does not exist. Stop.' %(basefile))
sys.exit(-1)

casefile = '%s/hofx.lo/%s' %(datadir, filename)
if(os.path.exists(casefile)):
self.caseio = ReadIODA2Obs(debug=debug, filename=casefile)
else:
print('casefile: <%s> does not exist. Stop.' %(casefile))
sys.exit(-1)

#------------------------------------------------------------------
def get_groups(self):
basegrps = self.baseio.get_groups()
casegrps = self.caseio.get_groups()
grps = []
for grp in casegrps:
if(grp in basegrps):
if(grp != 'MetaData'):
grps.append(grp)

return grps

#------------------------------------------------------------------
def get_variables(self, grp):
return self.caseio.groups[grp].variables

#------------------------------------------------------------------
def compare_int_variable(self, varname):
base_var = self.baseio.get_var(varname)
case_var = self.caseio.get_var(varname)

#print('len(case_var) = ', len(case_var))
#print('case_var) = ', case_var)
#print('len(base_var) = ', len(base_var))
#print('base_var) = ', base_var)

diff = case_var - base_var
nd = 0
for n in range(len(case_var)):
if(diff[n] != 0):
nd += 1

if(nd):
print('Difference of %s\n' %(varname))
print('\tLen(var): %d, nd: %d, mindif: %f, maxdif: %f' %(len(case_var), nd, np.min(diff), np.max(diff)))
for n in range(len(case_var)):
if(diff[n] != 0):
print('\tNo %d: c %d, b %d, diff: %d' %(n+1, case_var[n], base_var[n], diff[n]))
else:
print('var: %s has no difference' %(varname))

#------------------------------------------------------------------
def compare_float_variable(varname):
base_var = self.baseio.get_var(varname)
case_var = self.caseio.get_var(varname)

#print('len(case_var) = ', len(case_var))
#print('case_var) = ', case_var)
#print('len(base_var) = ', len(base_var))
#print('base_var) = ', base_var)

diff = case_var - base_var

nd = 0
for n in range(len(case_var)):
if(abs(diff[n]) > self.delt):
nd += 1

if(nd):
print('Difference of %s\n' %(varname))
print('\tLen(var): %d, nd: %d, mindif: %f, maxdif: %f' %(len(case_var), nd, np.min(diff), np.max(diff)))
for n in range(len(case_var)):
if(abs(diff[n]) > self.delt):
print('\tNo %d: c %f, b %f, diff: %f' %(n+1, case_var[n], base_var[n], diff[n]))
else:
print('var: %s has no difference' %(varname))

#------------------------------------------------------------------
def compare_2d_int_variable(varname):
base_var = self.baseio.get_2d_var(varname)
case_var = self.caseio.get_2d_var(varname)

#print('len(case_var) = ', len(case_var))
#print('case_var) = ', case_var)
#print('len(base_var) = ', len(base_var))
#print('base_var) = ', base_var)

diff = case_var - base_var
nd = 0
for j in range(len(case_var)):
for i in range(len(case_var[0])):
if(diff[j,i] != 0):
nd += 1

if(nd):
print('Difference of %s\n' %(varname))
print('\tLen(var): %d, nd: %d, mindif: %f, maxdif: %f' %(len(case_var), nd, np.min(diff), np.max(diff)))
for j in range(len(case_var)):
for i in range(len(case_var[0])):
if(diff[j,i] != 0):
print('\tJ: %d, I: %d: c %d, b %d, diff: %d' %(j, i, case_var[j,i], base_var[j,i], diff[j,i]))
else:
print('var: %s has no difference' %(varname))

#------------------------------------------------------------------
def compare_2d_float_variable(datadir, filename, varname):
base_var = self.baseio.get_2d_var(varname)
case_var = self.caseio.get_2d_var(varname)

diff = case_var - base_var

#print('len(case_var) = ', len(case_var))
#print('case_var) = ', case_var)
#print('len(base_var) = ', len(base_var))
#print('base_var) = ', base_var)

#print('\tbase_var.max: %f, base_var.min: %f' %(np.max(base_var[:,0]), np.min(base_var[:,0])))
#print('\tcase_var.max: %f, case_var.min: %f' %(np.max(case_var[:,0]), np.min(case_var[:,0])))
#print('\diffvar.max: %f, diffvar.min: %f' %(np.max(diffvar), np.min(diffvar)))

nd = 0
for j in range(len(case_var)):
for i in range(len(case_var[0])):
if(abs(diff[j,i]) > delt):
nd += 1

if(nd):
print('Difference of %s\n' %(varname))
print('\tbase_var.max: %f, base_var.min: %f' %(np.max(base_var[:,0]), np.min(base_var[:,0])))
print('\tcase_var.max: %f, case_var.min: %f' %(np.max(case_var[:,0]), np.min(case_var[:,0])))
print('\tLen(var): %d, nd: %d, mindif: %f, maxdif: %f' %(len(case_var), nd, np.min(diff), np.max(diff)))
for j in range(len(case_var)):
for i in range(len(case_var[0])):
if(abs(diff[j,i]) > delt):
print('\tJ: %d, I: %d: c %d, b %d, diff: %d' %(j, i, case_var[j,i], base_var[j,i], diff[j,i]))
else:
print('var: %s has no difference' %(varname))

#------------------------------------------------------------------
if __name__== '__main__':
debug = 0
datadir = '/scratch2/BMC/gsienkf/Wei.Huang/jedi/dev/build/intel/fv3-jedi/test/Data'
#filename = 'sfc_letkf-gfs_2020121500_m.nc4'
filename = 'sfc_letkf-gfs_2020121500_m.nc4'

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

for o, a in opts:
if o in ('--debug'):
debug = int(a)
elif o in ('--filename'):
filename = a
else:
print('option: ', a)
assert False, 'unhandled option'

#------------------------------------------------------------------
cof = Compare2ObsFiles()

#------------------------------------------------------------------
for filename in ['sfc_letkf-gfs_2020121500_m.nc4',
'aircraft_letkf-gfs_2020121500_m.nc4',
'scatwind_letkf-gfs_2020121500_m.nc4']:
cof.setup(datadir=datadir, filename=filename)
varlist = cof.get_variables('EffectiveQC0')
for grp in ctof.get_groups():
for var in varlist:
varname = '/%s/%s' %(grp, var)
if(grp in ['EffectiveQC0', 'EffectiveQC2', 'ObsType', 'PreQC', 'PreUseFlag']):
cof.compare_int_variable(varname)
else:
cof.compare_float_variable(varname)

#------------------------------------------------------------------
filename = 'amsua_n19_letkf-gfs_2020121500_m.nc4'
cof.setup(datadir=datadir, filename=filename)
varlist = cof.get_variables('EffectiveQC0')
for grp in ctof.get_groups():
for var in varlist:
varname = '/%s/%s' %(grp, var)
if(grp in ['EffectiveQC0', 'EffectiveQC2', 'ObsType', 'PreQC', 'PreUseFlag']):
cof.compare_2d_int_variable(varname)
else:
cof.compare_2d_float_variable(varname)

66 changes: 3 additions & 63 deletions obs+fv3-plot.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,6 @@
from matplotlib import cm
from mpl_toolkits.basemap import Basemap

from scipy_regridder import RegridFV3 as regridder
from readIODA2Obs import ReadIODA2Obs

sys.path.append('plot-utils')
Expand All @@ -22,14 +21,9 @@
debug = 1
output = 0
uselogp = 0
#casename = 'amsua'
#casename = 'surf'
casename = 'sondes'
casedir = '/work2/noaa/gsienkf/weihuang/jedi/case_study/%s/run_80.40t1n_36p' %(casename)
datadir = '%s/analysis/increment' %(casedir)
#obsfile = '%s/manmade-amsua-obs/amsua_n19_obs_2020121500_m.nc4' %(casedir)
#obsfile = '%s/ioda_v2_data/sfc_ps_obs_2020011006.nc4' %(casedir)
obsfile = '%s/ioda_v2_data/sondes_q_obs_2020011006.nc4' %(casedir)
casename = 'surf'
datadir = '/scratch2/BMC/gsienkf/Wei.Huang/jedi/dev/build/intel/fv3-jedi/test/Data/hofx.lo'
obsfile = '%s/sfc_letkf-gfs_2020121500_m.nc4' %(datadir)

opts, args = getopt.getopt(sys.argv[1:], '', ['debug=', 'output=',
'casename=', 'obsfile=', 'datadir=', 'uselogp='])
Expand All @@ -55,20 +49,6 @@
print('casename = ', casename)

#------------------------------------------------------------------------------
#griddir = '/work/noaa/gsienkf/weihuang/UFS-RNR-tools/JEDI.FV3-increments/grid/C96/'
griddir = '/work/noaa/gsienkf/weihuang/UFS-RNR-tools/JEDI.FV3-increments/grid/C48/'

datafiles = []
gridspecfiles = []
for ntile in range(1,7,1):
gridfile = '%sC48_grid.tile%s.nc' %(griddir, ntile)
gridspecfiles.append(gridfile)

datafile = '%s/20200110.030000.fv_core.res.tile%s.nc' %(datadir, ntile)
datafiles.append(datafile)

rg = regridder(debug=debug, datafiles=datafiles, gridspecfiles=gridspecfiles)

nlon = 360
nlat = nlon/2 + 1
varname = 'T'
Expand Down Expand Up @@ -131,43 +111,3 @@
#pt.plot(pvar, addmark=1, marker='x', size=1, color='green')
pt.plot(pvar, addmark=1, marker='+', size=5, color='red')

sys.exit(-1)
#------------------------------------------------------------------------------
#clevs = np.arange(-0.5, 0.51, 0.01)
#cblevs = np.arange(-0.5, 0.6, 0.1)
pt.set_clevs(clevs=clevs)
pt.set_cblevs(cblevs=cblevs)

#lons = [60, 200]
lons = [35, 75, 160, 205, 210]
#lons = [170]
for lon in lons:
pvar = var[:,:,lon]
title = '%s longitude %d' %(title_preix, lon)
pt.set_title(title)

imgname = '%s_lon_%d_logp.png' %(image_prefix, lon)
pt.set_imagename(imgname)
pt.plot_meridional_section_logp(pvar)

imgname = '%s_lon_%d_level.png' %(image_prefix, lon)
pt.set_imagename(imgname)
pt.plot_meridional_section(pvar)

#------------------------------------------------------------------------------
#lats = [50, 55]
lats = [-49, -31, -27, 12, 40]
#lats = [-41, -23, 0, 22, 41]
for lat in lats:
pvar = var[:,90+lat,:]
title = '%s latitude %d' %(title_preix, lat)
pt.set_title(title)

imgname = '%s_lat_%d_logp.png' %(image_prefix, lat)
pt.set_imagename(imgname)
pt.plot_zonal_section_logp(pvar)

imgname = '%s_lat_%d_level.png' %(image_prefix, lat)
pt.set_imagename(imgname)
pt.plot_zonal_section(pvar)

Loading

0 comments on commit 5608f73

Please sign in to comment.