Skip to content

Commit

Permalink
Inclusion of WAVEWATCHIII and Altimeter collocation (#8)
Browse files Browse the repository at this point in the history 
Add Model/Altimeter collocation included, modelSat_collocation.py
  • Loading branch information
@ricampos
ricampos authored May 17, 2022
1 parent 8dd7201 commit 8e2705b
Show file tree
Hide file tree
Showing 4 changed files with 436 additions and 157 deletions.
245 changes: 127 additions & 118 deletions validation/gridSatGlobal_Altimeter.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
PURPOSE:
Script to take altimeter tracks and collocate into regular
lat/lon grid (gridInfo.nc, generated by preprGridMask.py)
and hourly time interval (for different time intervals, edit atime array)
USAGE:
It runs one satellite mission, entered as argument (only the ID, sys.argv),
Expand Down Expand Up @@ -50,6 +51,7 @@
from mpl_toolkits.basemap import shiftgrid
import time
from calendar import timegm
import sys
import warnings
warnings.filterwarnings("ignore")
# netcdf format
Expand All @@ -64,9 +66,9 @@
# Maximum temporal distance (s) for pyresample weighted average
maxti=1800.
# Directory where AODN altimeter data is saved, downloaded using get_AODN_AltData.sh
dirs='/data/satellite/AODN_altm'
dirs='/work/noaa/marine/ricardo.campos/data/AODN/altimeter'
# Date interval
datemin='2000010100'; datemax='2021123123'
datemin='2021092300'; datemax='2021102500'

# Satellite missions available at AODN dataset, pick one as this code runs one satellite at a time!
s=np.int(sys.argv[1]) # argument satellite ID for satellite mission selection. s=0 is JASON3, s=1 is JASON2 etc. See list below.
Expand Down Expand Up @@ -167,122 +169,129 @@ def wf(pdist):
print(' Done reading and allocating satellite data '+sdname[s])
del ii

adatemin= np.double( (timegm( time.strptime(datemin, '%Y%m%d%H') )-float(timegm( time.strptime('1985010100', '%Y%m%d%H') ))) /(24.*3600.) )
adatemax= np.double( (timegm( time.strptime(datemax, '%Y%m%d%H') )-float(timegm( time.strptime('1985010100', '%Y%m%d%H') ))) /(24.*3600.) )

# Quality Control Check ----
indq = np.where( (aswhknstd<=max_swh_rms) & (asig0knstd<=max_sig0_rms) & (aswhknobs>=min_swh_numval[s]) & (aswhkqc<=max_swh_qc) & (ahsk>0.01) & (ahsk<hsmax) & (awnd>0.01) & (awnd<wspmax) & (ast>=float(timegm( time.strptime(datemin, '%Y%m%d%H') ))) & (ast<=float(timegm( time.strptime(datemax, '%Y%m%d%H') ))) )
del asig0knstd,aswhknobs,aswhknstd,aswhkqc
if np.any(indq):
if size(indq)>10:
ii=0
ast=np.double(np.copy(ast[indq[0]]))
ast=np.double(np.copy(ast)*24.*3600.+float(timegm( time.strptime('1985010100', '%Y%m%d%H') )))
aslat=np.copy(aslat[indq[0]]); aslon=np.copy(aslon[indq[0]])
ahsk=np.copy(ahsk[indq[0]]); ahskcal=np.copy(ahskcal[indq[0]]); ahsc=np.copy(ahsc[indq[0]])
awnd=np.copy(awnd[indq[0]]); awndcal=np.copy(awndcal[indq[0]])
# Collocated Arrays:
# final hourly time array. Combined with flatm and flonm, it represents the reference for sat collocation
atime = np.array(np.arange(np.double((np.round(ast.min()/3600.)*3600.)-3600.),np.double((np.round(ast.max()/3600.)*3600.)+3600.)+1,3600.)).astype('double')
#
ftime=np.double(np.zeros((ast.shape[0]*2),'d')); flat=np.zeros((ast.shape[0]*2),'f'); flon=np.zeros((ast.shape[0]*2),'f')
fhsk=np.zeros((ast.shape[0]*2),'f'); stdhsk=np.zeros((ast.shape[0]*2),'f'); counthsk=np.zeros((ast.shape[0]*2),'f')
fhskcal=np.zeros((ast.shape[0]*2),'f'); fhsc=np.zeros((ast.shape[0]*2),'f')
fwnd=np.zeros((ast.shape[0]*2),'f'); fwndcal=np.zeros((ast.shape[0]*2),'f')

# into the regular grid with pyresample kd tree
for t in range(0,atime.shape[0]):
indt = np.where( abs(ast[:]-atime[t]) < maxti )
if size(indt)>2:
prlon=np.copy(aslon[indt[0]]); prlon[prlon>180.]=prlon[prlon>180.]-360.
orig_def = pyresample.geometry.SwathDefinition(lons=prlon, lats=aslat[indt[0]]); del prlon
# By distance function wf
auxfhsk, auxstdhsk, auxcounthsk = pyresample.kd_tree.resample_custom(orig_def,ahsk[indt[0]],targ_def,radius_of_influence=dlim,weight_funcs=wf,fill_value=0,with_uncert=True,nprocs=npcs)
auxfhskcal = pyresample.kd_tree.resample_custom(orig_def,ahskcal[indt[0]],targ_def,radius_of_influence=dlim,weight_funcs=wf,fill_value=0,nprocs=npcs)
auxfhsc = pyresample.kd_tree.resample_custom(orig_def,ahsc[indt[0]],targ_def,radius_of_influence=dlim,weight_funcs=wf,fill_value=0,nprocs=npcs)
auxfwnd = pyresample.kd_tree.resample_custom(orig_def,awnd[indt[0]],targ_def,radius_of_influence=dlim,weight_funcs=wf,fill_value=0,nprocs=npcs)
auxfwndcal = pyresample.kd_tree.resample_custom(orig_def,awndcal[indt[0]],targ_def,radius_of_influence=dlim,weight_funcs=wf,fill_value=0,nprocs=npcs)
# print(' - ok '+repr(t))
indpqq = np.where( (auxfhskcal>0.01) & (auxfwnd>0.01) & (auxfhskcal<hsmax) & (auxfwnd<wspmax) )[0]
# allocate data into final array
ftime[ii:ii+size(indpqq)] = np.array(np.zeros(size(indpqq),'d')+atime[t]).astype('double')
flat[ii:ii+size(indpqq)] = np.array(flatm[indpqq]).astype('float')
flon[ii:ii+size(indpqq)] = np.array(flonm[indpqq]).astype('float')
fhsk[ii:ii+size(indpqq)] = np.array(auxfhsk[indpqq]).astype('float')
stdhsk[ii:ii+size(indpqq)] = np.array(auxstdhsk[indpqq]).astype('float')
counthsk[ii:ii+size(indpqq)] = np.array(auxcounthsk[indpqq]).astype('float')
fhskcal[ii:ii+size(indpqq)] = np.array(auxfhskcal[indpqq]).astype('float')
fhsc[ii:ii+size(indpqq)] = np.array(auxfhsc[indpqq]).astype('float')
fwnd[ii:ii+size(indpqq)] = np.array(auxfwnd[indpqq]).astype('float')
fwndcal[ii:ii+size(indpqq)] = np.array(auxfwndcal[indpqq]).astype('float')
ii=ii+size(indpqq)

del auxfhsk, auxstdhsk, auxcounthsk, auxfhskcal, auxfhsc, auxfwnd, auxfwndcal, indpqq
indq = np.where( (aswhknstd<=max_swh_rms) & (asig0knstd<=max_sig0_rms) & (aswhknobs>=min_swh_numval[s]) & (aswhkqc<=max_swh_qc) & (ahsk>0.01) & (ahsk<hsmax) & (awnd>0.01) & (awnd<wspmax) & (ast>=adatemin) & (ast<=adatemax) )
del asig0knstd,aswhknobs,aswhknstd,aswhkqc,adatemin,adatemax

if size(indq)>10:
ii=0
ast=np.double(np.copy(ast[indq[0]]))
ast=np.double(np.copy(ast)*24.*3600.+float(timegm( time.strptime('1985010100', '%Y%m%d%H') )))
aslat=np.copy(aslat[indq[0]]); aslon=np.copy(aslon[indq[0]])
ahsk=np.copy(ahsk[indq[0]]); ahskcal=np.copy(ahskcal[indq[0]]); ahsc=np.copy(ahsc[indq[0]])
awnd=np.copy(awnd[indq[0]]); awndcal=np.copy(awndcal[indq[0]])
# Collocated Arrays:
# final hourly time array. Combined with flatm and flonm, it represents the reference for sat collocation
atime = np.array(np.arange(np.double((np.round(ast.min()/3600.)*3600.)-3600.),np.double((np.round(ast.max()/3600.)*3600.)+3600.)+1,3600.)).astype('double')
#
ftime=np.double(np.zeros((ast.shape[0]*2),'d')); flat=np.zeros((ast.shape[0]*2),'f'); flon=np.zeros((ast.shape[0]*2),'f')
fhsk=np.zeros((ast.shape[0]*2),'f'); stdhsk=np.zeros((ast.shape[0]*2),'f'); counthsk=np.zeros((ast.shape[0]*2),'f')
fhskcal=np.zeros((ast.shape[0]*2),'f'); fhsc=np.zeros((ast.shape[0]*2),'f')
fwnd=np.zeros((ast.shape[0]*2),'f'); fwndcal=np.zeros((ast.shape[0]*2),'f')

# into the regular grid with pyresample kd tree
for t in range(0,atime.shape[0]):
indt = np.where( abs(ast[:]-atime[t]) < maxti )
if size(indt)>2:
prlon=np.copy(aslon[indt[0]]); prlon[prlon>180.]=prlon[prlon>180.]-360.
orig_def = pyresample.geometry.SwathDefinition(lons=prlon, lats=aslat[indt[0]]); del prlon
# By distance function wf
auxfhsk, auxstdhsk, auxcounthsk = pyresample.kd_tree.resample_custom(orig_def,ahsk[indt[0]],targ_def,radius_of_influence=dlim,weight_funcs=wf,fill_value=0,with_uncert=True,nprocs=npcs)
auxfhskcal = pyresample.kd_tree.resample_custom(orig_def,ahskcal[indt[0]],targ_def,radius_of_influence=dlim,weight_funcs=wf,fill_value=0,nprocs=npcs)
auxfhsc = pyresample.kd_tree.resample_custom(orig_def,ahsc[indt[0]],targ_def,radius_of_influence=dlim,weight_funcs=wf,fill_value=0,nprocs=npcs)
auxfwnd = pyresample.kd_tree.resample_custom(orig_def,awnd[indt[0]],targ_def,radius_of_influence=dlim,weight_funcs=wf,fill_value=0,nprocs=npcs)
auxfwndcal = pyresample.kd_tree.resample_custom(orig_def,awndcal[indt[0]],targ_def,radius_of_influence=dlim,weight_funcs=wf,fill_value=0,nprocs=npcs)
# print(' - ok '+repr(t))
indpqq = np.where( (auxfhskcal>0.01) & (auxfwnd>0.01) & (auxfhskcal<hsmax) & (auxfwnd<wspmax) )[0]
# allocate data into final array
ftime[ii:ii+size(indpqq)] = np.array(np.zeros(size(indpqq),'d')+atime[t]).astype('double')
flat[ii:ii+size(indpqq)] = np.array(flatm[indpqq]).astype('float')
flon[ii:ii+size(indpqq)] = np.array(flonm[indpqq]).astype('float')
fhsk[ii:ii+size(indpqq)] = np.array(auxfhsk[indpqq]).astype('float')
stdhsk[ii:ii+size(indpqq)] = np.array(auxstdhsk[indpqq]).astype('float')
counthsk[ii:ii+size(indpqq)] = np.array(auxcounthsk[indpqq]).astype('float')
fhskcal[ii:ii+size(indpqq)] = np.array(auxfhskcal[indpqq]).astype('float')
fhsc[ii:ii+size(indpqq)] = np.array(auxfhsc[indpqq]).astype('float')
fwnd[ii:ii+size(indpqq)] = np.array(auxfwnd[indpqq]).astype('float')
fwndcal[ii:ii+size(indpqq)] = np.array(auxfwndcal[indpqq]).astype('float')
ii=ii+size(indpqq)

del auxfhsk, auxstdhsk, auxcounthsk, auxfhskcal, auxfhsc, auxfwnd, auxfwndcal, indpqq

del indt
print('PyResample kdtree, hourly time, '+repr(t))

del ast, aslat, aslon, ahsk, ahskcal, ahsc, awnd, awndcal, indq, atime

print(' '); print(sdname[s]+' Done')

# Final quality control (double check)
ind=np.where( (fhsk<0.01) | (fhsk>hsmax) )
if np.any(ind):
fhsk[ind[0]]=np.nan; del ind

ind=np.where( (fhskcal<0.01) | (fhskcal>hsmax) )
if np.any(ind):
fhskcal[ind[0]]=np.nan; del ind

ind=np.where( (fhsc<0.01) | (fhsc>hsmax) )
if np.any(ind):
fhsc[ind[0]]=np.nan; del ind

ind=np.where( (fwnd<0.01) | (fwnd>wspmax) )
if np.any(ind):
fwnd[ind[0]]=np.nan; del ind

ind=np.where( (fwndcal<0.01) | (fwndcal>wspmax) )
if np.any(ind):
fwndcal[ind[0]]=np.nan; del ind

indf=np.where( (ftime>0.) & (fhsk>=0.0) )
ftime=np.array(ftime[indf[0]]).astype('double')
flat=np.array(flat[indf[0]]); flon=np.array(flon[indf[0]])
fhsk=np.array(fhsk[indf[0]]); fhskcal=np.array(fhskcal[indf[0]])
stdhsk=np.array(stdhsk[indf[0]]); counthsk=np.array(counthsk[indf[0]])
fhsc=np.array(fhsc[indf[0]])
fwnd=np.array(fwnd[indf[0]]); fwndcal=np.array(fwndcal[indf[0]])
print(sdname[s]+' . Array Size '+np.str(size(indf))); print(' ')
del indf

# Save netcdf
ncfile = nc.Dataset('AltimeterGridded_'+sdname[s]+'.nc', "w", format=fnetcdf)
ncfile.history="AODN Altimeter data on regular grid."
# create dimensions.
ncfile.createDimension('time' , ftime.shape[0] )
# create variables.
vflat = ncfile.createVariable('latitude',np.dtype('float32').char,('time'))
vflon = ncfile.createVariable('longitude',np.dtype('float32').char,('time'))
vft = ncfile.createVariable('stime',np.dtype('float64').char,('time'))
vfhsk = ncfile.createVariable('hsk',np.dtype('float32').char,('time'))
vstdhsk = ncfile.createVariable('stdhsk',np.dtype('float32').char,('time'))
vcounthsk = ncfile.createVariable('counthsk',np.dtype('float32').char,('time'))
vfhskcal = ncfile.createVariable('hskcal',np.dtype('float32').char,('time'))
vfhsc = ncfile.createVariable('hsc',np.dtype('float32').char,('time'))
vfwnd = ncfile.createVariable('wnd',np.dtype('float32').char,('time'))
vfwndcal = ncfile.createVariable('wndcal',np.dtype('float32').char,('time'))
# Assign units
vft.units = 'seconds since 1970-01-01 00:00:00'
vflat.units = 'degrees_north' ; vflon.units = 'degrees_east'
vfwnd.units = 'm/s' ; vfwndcal.units = 'm/s'
vfhsk.units = 'm'; vfhskcal.units = 'm'; vfhsc.units = 'm'
# Allocate Data
vflat[:] = flat; vflon[:] = flon; vft[:] = ftime
vfhsk[:] = fhsk; vstdhsk[:] = stdhsk; vcounthsk[:] = counthsk
vfhskcal[:] = fhskcal; vfhsc[:] = fhsc
vfwnd[:] = fwnd; vfwndcal[:] = fwndcal
ncfile.close()
print(' ')
print('netcdf ok ')

del vft, vflat, vflon, vfhsk, vstdhsk, vcounthsk, vfhskcal, vfhsc, vfwnd, vfwndcal
del indt
print('PyResample kdtree, hourly time, '+repr(t))

del ast, aslat, aslon, ahsk, ahskcal, ahsc, awnd, awndcal, indq, atime

print(' '); print(sdname[s]+' Done')

# Final quality control (double check)
ind=np.where( (fhsk<0.01) | (fhsk>hsmax) )
if np.any(ind):
fhsk[ind[0]]=np.nan; del ind

ind=np.where( (fhskcal<0.01) | (fhskcal>hsmax) )
if np.any(ind):
fhskcal[ind[0]]=np.nan; del ind

ind=np.where( (fhsc<0.01) | (fhsc>hsmax) )
if np.any(ind):
fhsc[ind[0]]=np.nan; del ind

ind=np.where( (fwnd<0.01) | (fwnd>wspmax) )
if np.any(ind):
fwnd[ind[0]]=np.nan; del ind

ind=np.where( (fwndcal<0.01) | (fwndcal>wspmax) )
if np.any(ind):
fwndcal[ind[0]]=np.nan; del ind

indf=np.where( (ftime>0.) & (fhsk>=0.0) )
ftime=np.array(ftime[indf[0]]).astype('double')
flat=np.array(flat[indf[0]]); flon=np.array(flon[indf[0]])
fhsk=np.array(fhsk[indf[0]]); fhskcal=np.array(fhskcal[indf[0]])
stdhsk=np.array(stdhsk[indf[0]]); counthsk=np.array(counthsk[indf[0]])
fhsc=np.array(fhsc[indf[0]])
fwnd=np.array(fwnd[indf[0]]); fwndcal=np.array(fwndcal[indf[0]])
print(sdname[s]+' . Array Size '+np.str(size(indf))); print(' ')
del indf

# Save netcdf
ncfile = nc.Dataset('AltimeterGridded_'+sdname[s]+'.nc', "w", format=fnetcdf)
ncfile.history="AODN Altimeter data on regular grid."
# create dimensions.
ncfile.createDimension('time' , ftime.shape[0] )
# create variables.
vflat = ncfile.createVariable('latitude',np.dtype('float32').char,('time'))
vflon = ncfile.createVariable('longitude',np.dtype('float32').char,('time'))
vft = ncfile.createVariable('stime',np.dtype('float64').char,('time'))
vfhsk = ncfile.createVariable('hsk',np.dtype('float32').char,('time'))
vstdhsk = ncfile.createVariable('stdhsk',np.dtype('float32').char,('time'))
vcounthsk = ncfile.createVariable('counthsk',np.dtype('float32').char,('time'))
vfhskcal = ncfile.createVariable('hskcal',np.dtype('float32').char,('time'))
vfhsc = ncfile.createVariable('hsc',np.dtype('float32').char,('time'))
vfwnd = ncfile.createVariable('wnd',np.dtype('float32').char,('time'))
vfwndcal = ncfile.createVariable('wndcal',np.dtype('float32').char,('time'))
# Assign units
vft.units = 'seconds since 1970-01-01 00:00:00'
vflat.units = 'degrees_north' ; vflon.units = 'degrees_east'
vfwnd.units = 'm/s' ; vfwndcal.units = 'm/s'
vfhsk.units = 'm'; vfhskcal.units = 'm'; vfhsc.units = 'm'
# Allocate Data
vflat[:] = flat; vflon[:] = flon; vft[:] = ftime
vfhsk[:] = fhsk; vstdhsk[:] = stdhsk; vcounthsk[:] = counthsk
vfhskcal[:] = fhskcal; vfhsc[:] = fhsc
vfwnd[:] = fwnd; vfwndcal[:] = fwndcal
ncfile.close()
print(' ')
print('netcdf ok ')

del vft, vflat, vflon, vfhsk, vstdhsk, vcounthsk, vfhskcal, vfhsc, vfwnd, vfwndcal

else:
sys.exit(' No satellite records within the given time/date range and quality control parameters selected.')


Loading

0 comments on commit 8e2705b

Please sign in to comment.