From aba047bb0e5b6b2c86135a7adbe4cb981f791ada Mon Sep 17 00:00:00 2001 From: Ming Chen Date: Tue, 12 May 2026 18:52:12 +0000 Subject: [PATCH 1/7] Use argparse and add configurable observation paths --- ww3tools/modelBuoy_collocation.py | 104 +++++++++++++++++++++--------- 1 file changed, 75 insertions(+), 29 deletions(-) diff --git a/ww3tools/modelBuoy_collocation.py b/ww3tools/modelBuoy_collocation.py index 3002b46..5a035cc 100755 --- a/ww3tools/modelBuoy_collocation.py +++ b/ww3tools/modelBuoy_collocation.py @@ -103,39 +103,85 @@ from time import strptime from calendar import timegm import wread +import argparse # netcdf format fnetcdf="NETCDF4" -# Paths -# ndbcp="/data/buoys/NDBC/wparam" -ndbcp="/work/noaa/marine/ricardo.campos/data/buoys/NDBC/ncformat/wparam" -# Copernicus buoys -# copernp="/data/buoys/Copernicus/wtimeseries" -copernp="/work/noaa/marine/ricardo.campos/data/buoys/Copernicus/wtimeseries" -print(' ') +# INPUT PARSING SECTION +def parse_args(): + parser = argparse.ArgumentParser( + description="Collocate WW3 point output with buoy observations." + ) -# Options of including grid and cyclone information -gridinfo=int(0); cyclonemap=int(0); wlist=[]; ftag=''; forecastds=0 -if len(sys.argv) < 2 : - sys.exit(' At least one argument (list of ww3 files) must be informed.') -if len(sys.argv) >= 2 : - # # import os; os.system("ls -d $PWD/*tab.nc > ww3list.txt &") - wlist=np.atleast_1d(np.loadtxt(sys.argv[1],dtype=str)) - ftag=str(sys.argv[1]).split('list')[1].split('.txt')[0] - print(' Reading ww3 list '+str(sys.argv[1])) - print(' Tag '+ftag) -if len(sys.argv) >= 3: - forecastds=int(sys.argv[2]) - if forecastds>0: - print(' Forecast-type data structure') -if len(sys.argv) >= 4: - gridinfo=str(sys.argv[3]) - print(' Using gridInfo '+gridinfo) -if len(sys.argv) >= 5: - cyclonemap=str(sys.argv[4]) - print(' Using cyclone map '+cyclonemap) -if len(sys.argv) > 5: - sys.exit(' Too many inputs') + parser.add_argument( + "ww3_list", + help="Text file containing WW3 point-output file paths." + ) + + parser.add_argument( + "-f", "--forecastds", + type=int, + default=0, + help="Forecast data structure flag. 0 = continuous/hindcast; >0 = forecast cycle/lead-time structure." + ) + + parser.add_argument( + "-g", "--gridinfo", + default="", + help="Optional gridInfo NetCDF file generated by prepGridMask.py." + ) + + parser.add_argument( + "-c", "--cyclonemap", + default="", + help="Optional CycloneMap NetCDF file generated by procyclmap.py. Requires --gridinfo." + ) + + parser.add_argument( + "--ndbcp", + default="/work/noaa/marine/ricardo.campos/data/buoys/NDBC/ncformat/wparam", + help="Directory containing NDBC buoy wave-parameter NetCDF files." + ) + + parser.add_argument( + "--copernp", + default="/work/noaa/marine/ricardo.campos/data/buoys/Copernicus/wtimeseries", + help="Directory containing Copernicus buoy time-series NetCDF files." + ) + + return parser.parse_args() + +args = parse_args() + +wlist = np.atleast_1d(np.loadtxt(args.ww3_list, dtype=str)) +ftag = str(args.ww3_list).split('list')[1].split('.txt')[0] + +forecastds = args.forecastds + +gridinfo = args.gridinfo if args.gridinfo else 0 +cyclonemap = args.cyclonemap if args.cyclonemap else 0 + +ndbcp = args.ndbcp +copernp = args.copernp + +# basic check +print(" Reading ww3 list " + str(args.ww3_list)) +print(" Tag " + str(ftag)) + +if forecastds > 0: + print(" Forecast-type data structure") + +if gridinfo != 0: + print(" Using gridInfo " + str(gridinfo)) + +if cyclonemap != 0 and gridinfo == 0: + sys.exit("ERROR: --cyclonemap requires --gridinfo.") + +if cyclonemap != 0: + print(" Using cyclone map " + str(cyclonemap)) + +print(" NDBC path : " + str(ndbcp)) +print(" Copernicus path : " + str(copernp)) # READ DATA print(" ") From 53687634bbf8caec4d685fab6130d6f5938e7b52 Mon Sep 17 00:00:00 2001 From: Ming Chen Date: Thu, 14 May 2026 18:36:26 +0000 Subject: [PATCH 2/7] fix NumPy year string conversion for buoy file paths --- ww3tools/modelBuoy_collocation.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/ww3tools/modelBuoy_collocation.py b/ww3tools/modelBuoy_collocation.py index 5a035cc..0b9c5af 100755 --- a/ww3tools/modelBuoy_collocation.py +++ b/ww3tools/modelBuoy_collocation.py @@ -95,6 +95,7 @@ import warnings; warnings.filterwarnings("ignore") import numpy as np +import os from matplotlib.mlab import * from pylab import * import xarray as xr @@ -417,7 +418,9 @@ def parse_args(): ahs=[];atm=[];atp=[];adm=[];atime=[] for y in yrange: - f=nc.Dataset(ndbcp+"/"+stname[b]+"h"+repr(y)+".nc") + fname = os.path.join(ndbcp, f'{stname[b]}h{y}.nc') + f = nc.Dataset(fname) + if 'wave_height' in f.variables.keys(): ahs = np.append(ahs,f.variables['wave_height'][:,0,0]) elif 'hs' in f.variables.keys(): From f1d5502991846aa968dfd60fc932fa0e1f24679a Mon Sep 17 00:00:00 2001 From: Ming Chen Date: Fri, 15 May 2026 15:55:08 +0000 Subject: [PATCH 3/7] remove unnecessary space and tab --- ww3tools/modelBuoy_collocation.py | 90 +++++++++++++++---------------- 1 file changed, 45 insertions(+), 45 deletions(-) diff --git a/ww3tools/modelBuoy_collocation.py b/ww3tools/modelBuoy_collocation.py index 0b9c5af..f6b2d46 100755 --- a/ww3tools/modelBuoy_collocation.py +++ b/ww3tools/modelBuoy_collocation.py @@ -17,26 +17,26 @@ Matchups of ww3 results and buoy data are generated for the same points (lat/lon) and time. Additional information can be collocated as well, such as water depth, - distance to the nearest coast, ocean names, forecast zones, and + distance to the nearest coast, ocean names, forecast zones, and cyclone information. - This code is designed for ww3 hindcasts or forecast with + This code is designed for ww3 hindcasts or forecast with consecutive cycles (overlapped time). The ww3 file(s) are not inserted directly but it is informed through a list, ww3list.txt (or any other - name), that is read as an argument. Multiple file names can be writen - in the list, ww3 results will be appended, depending on the data - structure selected (hindcast or forecast). The default is hindcast, so - arrays are directly appended in time. By entering a second argument - (any value greater than zero), the program assumes it is a forecast data - structure, i.e., the list contains consecutive cycles (each file is - one cycle) interpreted as another dimension (cycle time and + name), that is read as an argument. Multiple file names can be writen + in the list, ww3 results will be appended, depending on the data + structure selected (hindcast or forecast). The default is hindcast, so + arrays are directly appended in time. By entering a second argument + (any value greater than zero), the program assumes it is a forecast data + structure, i.e., the list contains consecutive cycles (each file is + one cycle) interpreted as another dimension (cycle time and forecast lead time). USAGE: Input WW3 point outputs are utilized (not grids), with option of - reading different formats: netcdf, text (bull and ts), + reading different formats: netcdf, text (bull and ts), or tar files with multiple bull or ts files. - For the observations, it uses two public buoy databases, - NDBC and Copernicus, which (at least one) must have been previously + For the observations, it uses two public buoy databases, + NDBC and Copernicus, which (at least one) must have been previously downloaded. See wfetchbuoy.py at ww3tools/downloadobs Users must edit ndbcp and copernp paths below. Python code can be run directly. At least one argument is needed, with @@ -44,8 +44,8 @@ can be written, for example: ls -d -1 $PWD/ww3.*.nc >> ww3list.txt The second argument (optional) is an integer with any value greater than zero for the program to assume it is a forecast data structure. - In addition to ww3list.txt and forecast data-shape, users can enter - two extra arguments, gridInfo and CycloneMap, generated by + In addition to ww3list.txt and forecast data-shape, users can enter + two extra arguments, gridInfo and CycloneMap, generated by prepGridMask.py and procyclmap.py, where the information of buoy's position (nearest grid point) will be extracted and included in the output netcdf file. @@ -84,8 +84,8 @@ included. Allow forecast array dimension (cycle-time and lead-time) for bull_tar 01/31/2023: Ricardo M. Campos, fix reshape array for forecast data (two time - dimensions), and check if variable names exist in the netcdf file (buoy - and ww3) to maximize the amount of matchups even when one variable is + dimensions), and check if variable names exist in the netcdf file (buoy + and ww3) to maximize the amount of matchups even when one variable is not available. PERSON OF CONTACT: @@ -192,7 +192,7 @@ def parse_args(): mlat=gridmask['latitude']; mlon=gridmask['longitude'] mask=gridmask['mask']; distcoast=gridmask['distcoast']; depth=gridmask['depth'] oni=gridmask['GlobalOceansSeas']; ocnames=gridmask['names_GlobalOceansSeas'] - hsmz=gridmask['HighSeasMarineZones']; hsmznames=gridmask['names_HighSeasMarineZones'] + hsmz=gridmask['HighSeasMarineZones']; hsmznames=gridmask['names_HighSeasMarineZones'] print(" GridInfo Ok. "+gridinfo) # Cyclone Information @@ -201,7 +201,7 @@ def parse_args(): clat=cycloneinfo['latitude']; clon=cycloneinfo['longitude'] cmap=cycloneinfo['cmap']; ctime=cycloneinfo['time'] cinfo=np.array(cycloneinfo['info'].split(':')[1].split(';')) - if np.array_equal(clat,mlat)==True & np.array_equal(clon,mlon)==True: + if np.array_equal(clat,mlat)==True & np.array_equal(clon,mlon)==True: print(" CycloneMap Ok. "+cyclonemap) else: sys.exit(' Error: Cyclone grid and Mask grid are different.') @@ -232,7 +232,7 @@ def parse_args(): if (mhs.shape[0]==result['hs'].shape[0]) and (np.size(stname)==np.size(result['station_name'])): if (stname==result['station_name']).all(): mtime=np.append(mtime,at) - mfcycle=np.append(mfcycle,fcycle) + mfcycle=np.append(mfcycle,fcycle) mhs=np.append(mhs,result['hs'],axis=1) mtp=np.append(mtp,result['tp'],axis=1) if 'dp' in result.keys(): @@ -312,7 +312,7 @@ def parse_args(): elif str(wlist[0]).split('/')[-1].split('.')[-1]=='nc': print(" Using ww3 netcdf point output format") - # netcdf point output file + # netcdf point output file for t in range(0,np.size(wlist)): try: f=nc.Dataset(str(wlist[t])) @@ -347,14 +347,14 @@ def parse_args(): if 'th1p' in f.variables.keys(): adp = np.array(f.variables['th1p'][:,:]).T else: - adp = np.array(np.copy(ahs*nan)) + adp = np.array(np.copy(ahs*nan)) - if 'tr' in f.variables.keys(): + if 'tr' in f.variables.keys(): atm = np.array(f.variables['tr'][:,:]).T else: - atm = np.array(np.copy(ahs*nan)) - - if 'fp' in f.variables.keys(): + atm = np.array(np.copy(ahs*nan)) + + if 'fp' in f.variables.keys(): auxtp = np.array(f.variables['fp'][:,:]).T indtp=np.where(auxtp>0.) if np.size(indtp)>0: @@ -425,19 +425,19 @@ def parse_args(): ahs = np.append(ahs,f.variables['wave_height'][:,0,0]) elif 'hs' in f.variables.keys(): ahs = np.append(ahs,f.variables['hs'][:,0,0]) - elif 'swh' in f.variables.keys(): + elif 'swh' in f.variables.keys(): ahs = np.append(ahs,f.variables['swh'][:,0,0]) - if 'average_wpd' in f.variables.keys(): + if 'average_wpd' in f.variables.keys(): atm = np.append(atm,f.variables['average_wpd'][:,0,0]) else: - atm = np.array(np.copy(ahs*nan)) + atm = np.array(np.copy(ahs*nan)) if 'dominant_wpd' in f.variables.keys(): atp = np.append(atp,f.variables['dominant_wpd'][:,0,0]) else: - atp = np.array(np.copy(ahs*nan)) - + atp = np.array(np.copy(ahs*nan)) + if 'mean_wave_dir' in f.variables.keys(): adm = np.append(adm,f.variables['mean_wave_dir'][:,0,0]) else: @@ -450,9 +450,9 @@ def parse_args(): else: lat[b] = nan - if 'longitude' in f.variables.keys(): + if 'longitude' in f.variables.keys(): lon[b] = f.variables['longitude'][:] - elif 'LONGITUDE' in f.variables.keys(): + elif 'LONGITUDE' in f.variables.keys(): lon[b] = f.variables['LONGITUDE'][:] else: lon[b] = nan @@ -470,8 +470,8 @@ def parse_args(): ahs = np.nanmean(f.variables['VHM0'][:,:],axis=1) elif 'VAVH' in f.variables.keys(): ahs = np.nanmean(f.variables['VAVH'][:,:],axis=1) - elif 'VGHS' in f.variables.keys(): - ahs = np.nanmean(f.variables['VGHS'][:,:],axis=1) + elif 'VGHS' in f.variables.keys(): + ahs = np.nanmean(f.variables['VGHS'][:,:],axis=1) elif 'significant_swell_wave_height' in f.variables.keys(): ahs = np.nanmean(f.variables['significant_swell_wave_height'][:,:],axis=1) elif 'sea_surface_significant_wave_height' in f.variables.keys(): @@ -499,7 +499,7 @@ def parse_args(): else: atp = ahs*nan - if 'VMDR' in f.variables.keys(): + if 'VMDR' in f.variables.keys(): adm = np.nanmean(f.variables['VMDR'][:,:],axis=1) else: adm = ahs*nan @@ -534,14 +534,14 @@ def parse_args(): if 'TIME' in f.variables.keys(): atime = np.array(f.variables['TIME'][:]*24*3600 + timegm( strptime('195001010000', '%Y%m%d%H%M') )).astype('double') - elif 'time' in f.variables.keys(): + elif 'time' in f.variables.keys(): atime = np.array(f.variables['time'][:]*24*3600 + timegm( strptime('195001010000', '%Y%m%d%H%M') )).astype('double') f.close(); del f except: ahs=[] - + if np.size(ahs)>0: # First layer of simple quality-control @@ -617,7 +617,7 @@ def parse_args(): ind=np.where((mtm>40.)|(mtm<0.0)) if np.size(ind)>0: mtm[ind]=np.nan; del ind - + ind=np.where((mtp>40.)|(mtp<0.0)) if np.size(ind)>0: mtp[ind]=np.nan; del ind @@ -743,7 +743,7 @@ def parse_args(): nmtime=np.zeros((unt.shape[0],mxsz),'double')*np.nan if cyclonemap!=0: nfcmap=np.zeros((mhs.shape[0],unt.shape[0],mxsz),'f')*np.nan - + nmtime[i,0:np.size(ind)]=np.array(mtime[ind]).astype('double') nmhs[:,i,:][:,0:np.size(ind)]=np.array(mhs[:,ind]) nmtm[:,i,:][:,0:np.size(ind)]=np.array(mtm[:,ind]) @@ -755,7 +755,7 @@ def parse_args(): nbtp[:,i,:][:,0:np.size(ind)]=np.array(btp[:,ind]) nbdm[:,i,:][:,0:np.size(ind)]=np.array(bdm[:,ind]) nbdp[:,i,:][:,0:np.size(ind)]=np.array(bdp[:,ind]) - if cyclonemap!=0: + if cyclonemap!=0: nfcmap[:,i,:][:,0:np.size(ind)]=np.array(fcmap[:,ind]) ind=np.where( (nmhs>0.0) & (nbhs>0.0) ) @@ -767,7 +767,7 @@ def parse_args(): if np.size(ind)>0: print(' Total amount of matchups model/buoy: '+repr(np.size(ind))) - # Save netcdf output file + # Save netcdf output file lon[lon>180.]=lon[lon>180.]-360. initime=str(time.gmtime(mtime.min())[0])+str(time.gmtime(mtime.min())[1]).zfill(2)+str(time.gmtime(mtime.min())[2]).zfill(2)+str(time.gmtime(mtime.min())[3]).zfill(2) fintime=str(time.gmtime(mtime.max())[0])+str(time.gmtime(mtime.max())[1]).zfill(2)+str(time.gmtime(mtime.max())[2]).zfill(2)+str(time.gmtime(mtime.max())[3]).zfill(2) @@ -815,11 +815,11 @@ def parse_args(): vbdp = ncfile.createVariable('obs_dp',np.dtype('float32').char,('buoypoints','time')) if gridinfo!=0: - vpdistcoast = ncfile.createVariable('distcoast',np.dtype('float32').char,('buoypoints')) - vpdepth = ncfile.createVariable('depth',np.dtype('float32').char,('buoypoints')) + vpdistcoast = ncfile.createVariable('distcoast',np.dtype('float32').char,('buoypoints')) + vpdepth = ncfile.createVariable('depth',np.dtype('float32').char,('buoypoints')) vponi = ncfile.createVariable('GlobalOceansSeas',np.dtype('float32').char,('buoypoints')) vocnames = ncfile.createVariable('names_GlobalOceansSeas',dtype('a25'),('GlobalOceansSeas')) - vphsmz = ncfile.createVariable('HighSeasMarineZones',np.dtype('float32').char,('buoypoints')) + vphsmz = ncfile.createVariable('HighSeasMarineZones',np.dtype('float32').char,('buoypoints')) vhsmznames = ncfile.createVariable('names_HighSeasMarineZones',dtype('a25'),('HighSeasMarineZones')) if cyclonemap!=0: if forecastds>0: @@ -873,7 +873,7 @@ def parse_args(): if cyclonemap!=0: vcinfo[:] = cinfo[:] if forecastds>0: - vcmap[:,:,:]=nfcmap[:,:,:] + vcmap[:,:,:]=nfcmap[:,:,:] else: vcmap[:,:]=fcmap[:,:] From 147eeb325763269828b4b376619a33b9879c5f2a Mon Sep 17 00:00:00 2001 From: Ming Chen Date: Sat, 16 May 2026 00:16:12 +0000 Subject: [PATCH 4/7] load wind speed and direction from NDBC NetCDF files and add to collocation output NetCDF files --- ww3tools/modelBuoy_collocation.py | 82 ++++++++++++++++++++++++++++++- 1 file changed, 81 insertions(+), 1 deletion(-) diff --git a/ww3tools/modelBuoy_collocation.py b/ww3tools/modelBuoy_collocation.py index f6b2d46..6ef08be 100755 --- a/ww3tools/modelBuoy_collocation.py +++ b/ww3tools/modelBuoy_collocation.py @@ -406,6 +406,8 @@ def parse_args(): btp=np.zeros((np.size(stname),np.size(mtime)),'f')*np.nan bdm=np.zeros((np.size(stname),np.size(mtime)),'f')*np.nan bdp=np.zeros((np.size(stname),np.size(mtime)),'f')*np.nan +bwsp=np.zeros((np.size(stname),np.size(mtime)),'f')*np.nan # wind speed +bwdir=np.zeros((np.size(stname),np.size(mtime)),'f')*np.nan # wind direction lat=np.zeros(np.size(stname),'f')*np.nan; lon=np.zeros(np.size(stname),'f')*np.nan # help reading NDBC buoys, divided by year yrange=np.array(np.arange(time.gmtime(mtime.min())[0],time.gmtime(mtime.min())[0]+1,1)).astype('int') @@ -415,7 +417,7 @@ def parse_args(): ahs=[] try: - ahs=[];atm=[];atp=[];adm=[];atime=[] + ahs=[];atm=[];atp=[];adm=[];atime=[];awdir=[];awsp=[] for y in yrange: fname = os.path.join(ndbcp, f'{stname[b]}h{y}.nc') @@ -443,6 +445,22 @@ def parse_args(): else: adm = np.array(np.copy(ahs*nan)) + # Wind speed + if 'wind_spd' in f.variables.keys(): + awsp = np.append(awsp, f.variables['wind_spd'][:,0,0]) + elif 'wspd' in f.variables.keys(): + awsp = np.append(awsp, f.variables['wspd'][:,0,0]) + else: + awsp = np.append(awsp, np.zeros(f.variables['time'][:].shape, 'f') * np.nan) + + # Wind direction + if 'wind_dir' in f.variables.keys(): + awdir = np.append(awdir, f.variables['wind_dir'][:,0,0]) + elif 'wdir' in f.variables.keys(): + awdir = np.append(awdir, f.variables['wdir'][:,0,0]) + else: + awdir = np.append(awdir, np.zeros(f.variables['time'][:].shape, 'f') * np.nan) + if 'latitude' in f.variables.keys(): lat[b] = f.variables['latitude'][:] elif 'LATITUDE' in f.variables.keys(): @@ -565,6 +583,16 @@ def parse_args(): if np.size(indq)>0: adp[indq]=np.nan; del indq + # Wind speed quality control + indq=np.where((awsp>100.)|(awsp<0.0)) + if np.size(indq)>0: + awsp[indq]=np.nan; del indq + + # Wind direction quality control + indq=np.where((awdir>360.)|(awdir<-180.)) + if np.size(indq)>0: + awdir[indq]=np.nan; del indq + c=0 for t in range(0,np.size(mtime)): indt=np.where(np.abs(atime-mtime[t])<1800.) @@ -580,6 +608,10 @@ def parse_args(): bdm[b,t] = np.nanmean(adm[indt[0]][adm[indt[0]].mask==False]) if np.any(adp[indt[0]].mask==False): bdp[b,t] = np.nanmean(adp[indt[0]][adp[indt[0]].mask==False]) + if np.any(awsp[indt[0]].mask==False): + bwsp[b,t] = np.nanmean(awsp[indt[0]][awsp[indt[0]].mask==False]) + if np.any(awdir[indt[0]].mask==False): + bwdir[b,t] = np.nanmean(awdir[indt[0]][awdir[indt[0]].mask==False]) del indt @@ -610,6 +642,14 @@ def parse_args(): if np.size(ind)>0: bdp[ind]=np.nan; del ind +ind=np.where((bwsp>100.)|(bwsp<0.0)) +if np.size(ind)>0: + bwsp[ind]=np.nan; del ind + +ind=np.where((bwdir>360.)|(bwdir<-180.)) +if np.size(ind)>0: + bwdir[ind]=np.nan; del ind + ind=np.where((mhs>30.)|(mhs<0.0)) if np.size(ind)>0: mhs[ind]=np.nan; del ind @@ -646,6 +686,8 @@ def parse_args(): btp=np.array(btp[ind[0],:]) bdm=np.array(bdm[ind[0],:]) bdp=np.array(bdp[ind[0],:]) + bwsp=np.array(bwsp[ind[0],:]) + bwdir=np.array(bwdir[ind[0],:]) else: sys.exit(' Error: No matchups Model/Buoy available.') @@ -690,6 +732,8 @@ def parse_args(): btp=np.array(btp[ind[0],:]) bdm=np.array(bdm[ind[0],:]) bdp=np.array(bdp[ind[0],:]) + bwsp=np.array(bwsp[ind[0],:]) + bwdir=np.array(bwdir[ind[0],:]) pdistcoast=np.array(pdistcoast[ind[0]]) pdepth=np.array(pdepth[ind[0]]) poni=np.array(poni[ind[0]]) @@ -740,6 +784,8 @@ def parse_args(): nbtp=np.zeros((mhs.shape[0],unt.shape[0],mxsz),'f')*np.nan nbdm=np.zeros((mhs.shape[0],unt.shape[0],mxsz),'f')*np.nan nbdp=np.zeros((mhs.shape[0],unt.shape[0],mxsz),'f')*np.nan + nbwsp=np.zeros((mhs.shape[0],unt.shape[0],mxsz),'f')*np.nan + nbwdir=np.zeros((mhs.shape[0],unt.shape[0],mxsz),'f')*np.nan nmtime=np.zeros((unt.shape[0],mxsz),'double')*np.nan if cyclonemap!=0: nfcmap=np.zeros((mhs.shape[0],unt.shape[0],mxsz),'f')*np.nan @@ -754,6 +800,8 @@ def parse_args(): nbtm[:,i,:][:,0:np.size(ind)]=np.array(btm[:,ind]) nbtp[:,i,:][:,0:np.size(ind)]=np.array(btp[:,ind]) nbdm[:,i,:][:,0:np.size(ind)]=np.array(bdm[:,ind]) + nbwsp[:,i,:][:,0:np.size(ind)]=np.array(bwsp[:,ind]) + nbwdir[:,i,:][:,0:np.size(ind)]=np.array(bwdir[:,ind]) nbdp[:,i,:][:,0:np.size(ind)]=np.array(bdp[:,ind]) if cyclonemap!=0: nfcmap[:,i,:][:,0:np.size(ind)]=np.array(fcmap[:,ind]) @@ -800,6 +848,8 @@ def parse_args(): vbtp = ncfile.createVariable('obs_tp',np.dtype('float32').char,('buoypoints','fcycle','time')) vbdm = ncfile.createVariable('obs_dm',np.dtype('float32').char,('buoypoints','fcycle','time')) vbdp = ncfile.createVariable('obs_dp',np.dtype('float32').char,('buoypoints','fcycle','time')) + vbwsp = ncfile.createVariable('obs_wsp', np.dtype('float32').char, ('buoypoints','fcycle','time')) + vbwdir = ncfile.createVariable('obs_wdir', np.dtype('float32').char, ('buoypoints','fcycle','time')) else: ncfile.createDimension('time', bhs.shape[1] ) vt = ncfile.createVariable('time',np.dtype('float64').char,('time')) @@ -813,6 +863,8 @@ def parse_args(): vbtp = ncfile.createVariable('obs_tp',np.dtype('float32').char,('buoypoints','time')) vbdm = ncfile.createVariable('obs_dm',np.dtype('float32').char,('buoypoints','time')) vbdp = ncfile.createVariable('obs_dp',np.dtype('float32').char,('buoypoints','time')) + vbwsp = ncfile.createVariable('obs_wsp', np.dtype('float32').char, ('buoypoints','time')) + vbwdir = ncfile.createVariable('obs_wdir', np.dtype('float32').char, ('buoypoints','time')) if gridinfo!=0: vpdistcoast = ncfile.createVariable('distcoast',np.dtype('float32').char,('buoypoints')) @@ -835,9 +887,33 @@ def parse_args(): vmtp.units='s'; vbtp.units='s' vmdm.units='degrees'; vbdm.units='degrees' vmdp.units='degrees'; vbdp.units='degrees' + vbwsp.units = 'm s-1' + vbwdir.units = 'degrees' if gridinfo!=0: vpdepth.units='m'; vpdistcoast.units='km' + # Assign long names + vstname.long_name = 'Buoy Station ID' + + vlat.long_name = 'Latitude' + vlon.long_name = 'Longitude' + + vt.long_name = 'Valid Time' + + vmhs.long_name = 'Model Significant Wave Height' + vmtm.long_name = 'Model Mean Wave Period' + vmtp.long_name = 'Model Peak Wave Period' + vmdm.long_name = 'Model Mean Wave Direction' + vmdp.long_name = 'Model Peak Wave Direction' + + vbhs.long_name = 'Observed Significant Wave Height' + vbtm.long_name = 'Observed Mean Wave Period' + vbtp.long_name = 'Observed Peak Wave Period' + vbdm.long_name = 'Observed Mean Wave Direction' + vbdp.long_name = 'Observed Peak Wave Direction' + vbwsp.long_name = 'Observed Wind Speed' + vbwdir.long_name = 'Observed Wind Direction' + # Allocate Data vstname[:]=stname[:]; vlat[:] = lat[:]; vlon[:] = lon[:] if forecastds>0: @@ -852,6 +928,8 @@ def parse_args(): vbtp[:,:,:]=nbtp[:,:,:] vbdm[:,:,:]=nbdm[:,:,:] vbdp[:,:,:]=nbdp[:,:,:] + vbwsp[:,:,:]=nbwsp[:,:,:] + vbwdir[:,:,:]=nbwdir[:,:,:] else: vt[:]=mtime[:] vmhs[:,:]=mhs[:,:] @@ -864,6 +942,8 @@ def parse_args(): vbtp[:,:]=btp[:,:] vbdm[:,:]=bdm[:,:] vbdp[:,:]=bdp[:,:] + vbwsp[:,:]=bwsp[:,:] + vbwdir[:,:]=bwdir[:,:] if gridinfo!=0: vpdistcoast[:]=pdistcoast[:] From 4451b596d80ac827c5347ac3a8b00ba34f26ffa6 Mon Sep 17 00:00:00 2001 From: Ming Chen Date: Sat, 16 May 2026 00:50:52 +0000 Subject: [PATCH 5/7] clear all unnecessary space and tab in wread.py --- ww3tools/wread.py | 132 +++++++++++++++++++++++----------------------- 1 file changed, 66 insertions(+), 66 deletions(-) diff --git a/ww3tools/wread.py b/ww3tools/wread.py index 97dcaf6..0f79ad7 100755 --- a/ww3tools/wread.py +++ b/ww3tools/wread.py @@ -100,7 +100,7 @@ def readconfig(fname): the full path '/home/user/ww3tools.yaml' """ - try: + try: with open(fname, 'r') as file: wconfig = yaml.safe_load(file) except: @@ -161,11 +161,11 @@ def mask(*args): if 'distcoast' in f.variables.keys(): result['distcoast'] = np.array(f.variables['distcoast'][:,:]) if 'depth' in f.variables.keys(): - result['depth'] = np.array(f.variables['depth'][:,:]) + result['depth'] = np.array(f.variables['depth'][:,:]) if 'GlobalOceansSeas' in f.variables.keys(): - result['GlobalOceansSeas'] = np.array(f.variables['GlobalOceansSeas'][:,:]) + result['GlobalOceansSeas'] = np.array(f.variables['GlobalOceansSeas'][:,:]) if 'HighSeasMarineZones' in f.variables.keys(): - result['HighSeasMarineZones'] = np.array(f.variables['HighSeasMarineZones'][:,:]) + result['HighSeasMarineZones'] = np.array(f.variables['HighSeasMarineZones'][:,:]) if 'names_GlobalOceansSeas' in f.variables.keys(): result['names_GlobalOceansSeas'] = f.variables['names_GlobalOceansSeas'][:] if 'names_HighSeasMarineZones' in f.variables.keys(): @@ -192,7 +192,7 @@ def cyclonemap(*args): f=nc.MFDataset(fname, aggdim='time') at=f.variables['time'][:]; adate=[] for j in range(0,at.shape[0]): - adate=np.append(adate,date2num(datetime.datetime(time.gmtime(at[j])[0],time.gmtime(at[j])[1],time.gmtime(at[j])[2],time.gmtime(at[j])[3],time.gmtime(at[j])[4]))) + adate=np.append(adate,date2num(datetime.datetime(time.gmtime(at[j])[0],time.gmtime(at[j])[1],time.gmtime(at[j])[2],time.gmtime(at[j])[3],time.gmtime(at[j])[4]))) # -------- # build dictionary result={'latitude':np.array(f.variables['lat'][:]),'longitude':np.array(f.variables['lon'][:]), @@ -207,7 +207,7 @@ def cyclonemap(*args): del result -# ================= OBSERVATIONS ================= +# ================= OBSERVATIONS ================= # --- Buoys --- # Observations NDBC, netcdf format def tseriesnc_ndbc(fname=None,anh=None): @@ -215,7 +215,7 @@ def tseriesnc_ndbc(fname=None,anh=None): Observations NDBC, time series/table, netcdf format one file per buoy and year Input: file name (example: 46047h2016.nc), and anemometer height (optional) - Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, + Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, and arrays sst,mslp,dwp,tmp,gst(10-m height),wsp(10-m height),wdir,hs,tm,tp,dm ''' if fname==None: @@ -283,8 +283,8 @@ def tseriesnc_ndbc(fname=None,anh=None): result={'latitude':blat,'longitude':blon, 'time':btime,'date':ds['time'].values[:], 'sst':bsst, 'mslp':bmslp, 'dewpt_temp':bdwp, - 'air_temp':btmp, 'gust':bgst, 'wind_spd':bwsp, - 'wind_dir':bwdir, 'hs':bhs, 'tm':btm, + 'air_temp':btmp, 'gust':bgst, 'wind_spd':bwsp, + 'wind_dir':bwdir, 'hs':bhs, 'tm':btm, 'tp':btp, 'dm':bdm, 'tm':btm} return result @@ -314,7 +314,7 @@ def tseriestxt_ndbc(fname=None,anh=None): ds['date']=pd.to_datetime(ds['date'],format='%Y %m %d %H') for i in range(0,btime.shape[0]): - btime[i]=double(ds['date'][i].timestamp()) + btime[i]=double(ds['date'][i].timestamp()) except: sys.exit(" Cannot open "+fname) @@ -349,7 +349,7 @@ def tseriestxt_ndbc(fname=None,anh=None): if 'W' in auxlatlon: blon=-float(auxlatlon[8:16]) else: - blon=float(auxlatlon[8:16]) + blon=float(auxlatlon[8:16]) except: if anh==None: @@ -388,8 +388,8 @@ def tseriestxt_ndbc(fname=None,anh=None): result={'latitude':blat,'longitude':blon, 'time':btime,'date':ds['date'].values[:], 'sst':bsst, 'mslp':bmslp, 'dewpt_temp':bdwp, - 'air_temp':btmp, 'gust':bgst, 'wind_spd':bwsp, - 'wind_dir':bwdir, 'hs':bhs, 'tm':btm, + 'air_temp':btmp, 'gust':bgst, 'wind_spd':bwsp, + 'wind_dir':bwdir, 'hs':bhs, 'tm':btm, 'tp':btp, 'dm':bdm, 'tm':btm} return result @@ -401,7 +401,7 @@ def tseriesnc_copernicus(*args): Observations NDBC, time series/table, netcdf format one file per buoy Input: file name (example: GL_TS_MO_41004.nc) - Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, + Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, and arrays with the environmental variables available. ''' if len(args) == 1: @@ -422,8 +422,8 @@ def tseriesnc_copernicus(*args): result={'latitude':np.array(blat),'longitude':np.array(blon), 'time':btime,'date':ds['TIME'].values[:]} - if 'DEPH' in ds.keys(): - bdepth = np.nanmean(ds['DEPH'].values[:,:],axis=1) # Depth + if 'DEPH' in ds.keys(): + bdepth = np.nanmean(ds['DEPH'].values[:,:],axis=1) # Depth result['depth']=np.array(bdepth) if 'VHM0' in ds.keys(): @@ -435,7 +435,7 @@ def tseriesnc_copernicus(*args): bhs[(bhs<0.1)|(bhs>20)]=np.nan result['hs']=np.array(bhs) - if 'VAVH' in ds.keys(): + if 'VAVH' in ds.keys(): bvavh = np.nanmean(ds['VAVH'].values[:,:],axis=1) # H 1/3 vavh bvavh[(bvavh<0.1)|(bvavh>20)]=np.nan result['hs_vavh']=np.array(bvavh) @@ -459,8 +459,8 @@ def tseriesnc_copernicus(*args): btp[(btp<1)|(btp>30)]=np.nan result['tp']=np.array(btp) - if 'TEMP' in ds.keys(): - bsst = np.nanmean(ds['TEMP'].values[:,:],axis=1) # SST + if 'TEMP' in ds.keys(): + bsst = np.nanmean(ds['TEMP'].values[:,:],axis=1) # SST bsst[np.abs(bsst)>70]=np.nan result['sst']=np.array(bsst) @@ -486,13 +486,13 @@ def tseriesnc_copernicus(*args): result['gust']=np.array(bgst) if 'WSPD' in ds.keys(): - bwsp = np.nanmean(ds['WSPD'].values[:,:],axis=1) # wind speed + bwsp = np.nanmean(ds['WSPD'].values[:,:],axis=1) # wind speed bwsp=np.copy(((10./4.0)**(0.12))*bwsp) # conversion to 10m, approximation DNVGL C-205 Table 2-1 bwsp[(bwsp<0)|(bwsp>150)]=np.nan result['wind_spd']=np.array(bwsp) if 'WDIR' in ds.keys(): - bwdir = np.nanmean(ds['WDIR'].values[:,:],axis=1) # wind direction + bwdir = np.nanmean(ds['WDIR'].values[:,:],axis=1) # wind direction bwdir[(bwdir<-180)|(bwdir>360)]=np.nan result['wind_dir']=np.array(bwdir) @@ -501,7 +501,7 @@ def tseriesnc_copernicus(*args): bhcmax[(bhcmax<0)|(bhcmax>30)]=np.nan result['hc_max']=np.array(bhcmax) - if 'VMDR' in ds.keys(): + if 'VMDR' in ds.keys(): bdm = np.nanmean(ds['VMDR'].values[:,:],axis=1) # Mean direction bdm[(bdm<-180)|(bdm>360)]=np.nan result['dm']=np.array(bdm) @@ -520,7 +520,7 @@ def tseriesnc_cdip(*args): Observations CDIP, time series/table, netcdf format one file per buoy Input: file name (example: CDIP_buoy_144_historic.nc) - Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, + Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, and arrays hs,tm,tp,tz,dp ''' @@ -569,7 +569,7 @@ def tseriesnc_microswift(*args): Observations microSWIFT, time series/table, netcdf format one file per buoy Input: file name (example: microSWIFT041_HurricaneLee_Sep2023.nc) - Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, + Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, and arrays hs,tp,dp,sst ''' @@ -614,7 +614,7 @@ def tseries_spotter(*args): multiple buoys in the same file Input: file name (example: campos_hurricane_spotters_2022_spectra_with_direction_coefficients.pkl) and station ID. - Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, + Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, and arrays hs,tm,tp,dm,dp ''' @@ -684,7 +684,7 @@ def tseriesnc_dwsd(*args): multiple buoys in the same file Input: file name (example: LDL_AtlanticHurricane2021_4a61_8818_8454.nc) and station ID. - Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, + Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, and arrays wsp (10-m height),wdir,slp,sst,hs,tm,tp,dp ''' @@ -751,7 +751,7 @@ def tseriesnc_dwsd(*args): result={'latitude':blat,'longitude':blon, 'time':btime,'date':bdate, - 'wind_spd':bwsp, 'wind_dir':bwdir, 'slp':bslp, 'sst':bsst, + 'wind_spd':bwsp, 'wind_dir':bwdir, 'slp':bslp, 'sst':bsst, 'hs':bhs, 'tm':btm, 'tp':btp, 'dp':bdp} return result @@ -764,7 +764,7 @@ def tseriesnc_saildrone(*args): Observations saildrones, time series/table, netcdf format one file per saildrone Input: file name (example: sd1031_hurricane_2024_af17_91e6_4d06.nc) - Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, + Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, and arrays wsp (10-m height),wdir,sst,slp,rh,tmp,hs,tp ''' @@ -800,7 +800,7 @@ def tseriesnc_saildrone(*args): # dictionary result={'latitude':blat,'longitude':blon, 'time':btime,'date':ds['time'].values[:], - 'sst':bsst, 'mlp':bslp, 'rh':brh, 'air_temp':btmp, + 'sst':bsst, 'mlp':bslp, 'rh':brh, 'air_temp':btmp, 'hs':bhs, 'tp':btp} if 'WIND_SPEED_MEAN' in ds.keys(): @@ -832,7 +832,7 @@ def tseriesnc_wsra(*args): ''' Observations WSRA_L4, time series/table, netcdf format Input: file name (example: WSRA-L4-20220924H1.nc) - Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, + Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, and arrays wsp(10-m height),wdir,sst,slp,rh,tmp,hs,tp ''' @@ -885,7 +885,7 @@ def tseriesnc_wsra(*args): # rdp[(rdp<-180)|(rdp>360)]=np.nan # rrflr[(rrflr<0)|(rrflr>200)]=np.nan # rrflrm[(rrflrm<0)|(rrflrm>50)]=np.nan - # + # # pralt[(pralt<1000)|(pralt>4000)]=np.nan # pseed[(pseed<80)|(pseed>250)]=np.nan # wcroll[(wcroll<-2.5)|(wcroll>2.5)]=np.nan @@ -893,9 +893,9 @@ def tseriesnc_wsra(*args): # dictionary result={'latitude':rlat,'longitude':rlon, 'time':rtime,'date':ds['time'].values[:], - 'wind_spd':rwsp, 'wind_dir':rwdir, + 'wind_spd':rwsp, 'wind_dir':rwdir, 'hs':rhs, 'dhs':rdwh, 'dp':rdp, 'rainfall_rate':rrflr, 'rainfall_rate_median':rrflrm, - 'porient':porient, 'wcroll':wcroll, 'pralt':pralt, 'pseed':pseed, + 'porient':porient, 'wcroll':wcroll, 'pralt':pralt, 'pseed':pseed, 'hurricane_eye_distance':rhed, 'diff_pcourse':diffcourse } return result @@ -1072,7 +1072,7 @@ def tseriestxt_ww3(*args): WAVEWATCH III, time series/table, text tab format This file format has all point outputs (results) in the same file (not divided by point/buoy). Input: file name (example: tab50.ww3), and number of point ouputs (example: 4) - Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, + Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, and arrays with the wave variables available. Inside the dictionary, the arrays of wave variables have dimension (point_outputs, time). ''' @@ -1091,7 +1091,7 @@ def tseriestxt_ww3(*args): tt = int(np.size(mcontent)/(7+tnb)+1) myear = []; mmonth = [] ; mday = [] ; mhour = []; mmin = [] - mlon = np.zeros((tnb,tt),'f'); mlat = np.zeros((tnb,tt),'f'); mhs = np.zeros((tnb,tt),'f'); mL = np.zeros((tnb,tt),'f') + mlon = np.zeros((tnb,tt),'f'); mlat = np.zeros((tnb,tt),'f'); mhs = np.zeros((tnb,tt),'f'); mL = np.zeros((tnb,tt),'f') mtm = np.zeros((tnb,tt),'f'); mdm = np.zeros((tnb,tt),'f'); mspr = np.zeros((tnb,tt),'f') atp = np.zeros((tnb,tt),'f'); mdp = np.zeros((tnb,tt),'f'); mpspr = np.zeros((tnb,tt),'f') for i in range(0,tt): @@ -1114,7 +1114,7 @@ def tseriestxt_ww3(*args): mpspr[k,i] = mcontent[j+tnb+1+k].strip().split()[9] mtp = np.zeros((atp.shape[0],atp.shape[1]),'f')*np.nan - for i in range(0,mtp.shape[0]): + for i in range(0,mtp.shape[0]): #mtp[i,atp[i,:]>0.0] = 1./atp[i,atp[i,:]>0.0] indtp=np.where(atp[i,:]>0.0) if np.size(indtp)>0: @@ -1138,7 +1138,7 @@ def tseriesnc_ww3(*args): ''' WAVEWATCH III, time series/table, netcdf format Input: file name (example: ww3gefs.20160928_tab.nc), and station name (example: 41002) - Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, + Output: dictionary containing the arrays: time(seconds since 1970),time(datetime64),lat,lon, and arrays with the wave variables available. ''' if len(args) == 2: @@ -1155,7 +1155,7 @@ def tseriesnc_ww3(*args): else: mtime = np.array(f.variables['time'][:]*24*3600 + timegm( strptime(str(f.variables['time'].units).split(' ')[2][0:4]+'01010000', '%Y%m%d%H%M') )).astype('double') f.close(); del f - + auxstationname=ds['station_name'].values[:,:]; stationname=[] for i in range(0,auxstationname.shape[0]): stationname=np.append(stationname,"".join(np.array(auxstationname[i,:]).astype('str'))) @@ -1170,7 +1170,7 @@ def tseriesnc_ww3(*args): mlon = np.nanmean(ds['longitude'].values[:,inds]) # dictionary result={'latitude':np.array(mlat),'longitude':np.array(mlon), - 'time':mtime,'date':ds['time'].values[:]} + 'time':mtime,'date':ds['time'].values[:]} if 'hs' in ds.keys(): mhs = ds['hs'].values[:,inds] @@ -1266,7 +1266,7 @@ def bull(*args): del hour,day,month,year for j in range(0,at.shape[0]): - adate=np.append(adate,date2num(datetime.datetime(time.gmtime(at[j])[0],time.gmtime(at[j])[1],time.gmtime(at[j])[2],time.gmtime(at[j])[3],time.gmtime(at[j])[4]))) + adate=np.append(adate,date2num(datetime.datetime(time.gmtime(at[j])[0],time.gmtime(at[j])[1],time.gmtime(at[j])[2],time.gmtime(at[j])[3],time.gmtime(at[j])[4]))) # -------- ahs=[]; atp=[] @@ -1280,7 +1280,7 @@ def bull(*args): # build dictionary result={'time':np.array(at).astype('double'),'date':np.array(adate).astype('double'), - 'latitude':alat,'longitude':alon,'station_name':stname, + 'latitude':alat,'longitude':alon,'station_name':stname, 'hs':np.array(ahs),'tp':np.array(atp)} # GFS, HAFS, and other formats @@ -1317,7 +1317,7 @@ def bull(*args): del hour,day,month,year for j in range(0,at.shape[0]): - adate=np.append(adate,date2num(datetime.datetime(time.gmtime(at[j])[0],time.gmtime(at[j])[1],time.gmtime(at[j])[2],time.gmtime(at[j])[3],time.gmtime(at[j])[4]))) + adate=np.append(adate,date2num(datetime.datetime(time.gmtime(at[j])[0],time.gmtime(at[j])[1],time.gmtime(at[j])[2],time.gmtime(at[j])[3],time.gmtime(at[j])[4]))) # -------- for j in range(7,np.size(lines)-8): @@ -1326,21 +1326,21 @@ def bull(*args): # aux... is organizing the partitions. Ready for future versions (not included yet) auxhs=[] - for k in range(0,4): + for k in range(0,4): if len(str(lines[j][int(iauxhs[0]+18*k):int(iauxhs[1]+18*k)]).replace(' ', '')): auxhs=np.append(auxhs,float(lines[j][int(iauxhs[0]+18*k):int(iauxhs[1]+18*k)])) else: auxhs=np.append(auxhs,np.nan) auxtp=[] - for k in range(0,4): + for k in range(0,4): if len(str(lines[j][int(iauxtp[0]+18*k):int(iauxtp[1]+18*k)]).replace(' ', '')): auxtp=np.append(auxtp,float(lines[j][int(iauxtp[0]+18*k):int(iauxtp[1]+18*k)])) else: auxtp=np.append(auxtp,np.nan) auxdp=[] - for k in range(0,4): + for k in range(0,4): if len(str(lines[j][int(iauxdp[0]+18*k):int(iauxdp[1]+18*k)]).replace(' ', '')): auxdp=np.append(auxdp,float(lines[j][int(iauxdp[0]+18*k):int(iauxdp[1]+18*k)])) else: @@ -1402,7 +1402,7 @@ def bull_tar(*args): # GEFS specific format if 'gefs' in str(fname).split('/')[-1]: iauxhs=[10,15];iauxtp=[28,33] - + for t in range(0,np.size(tar.getmembers())): # station names stname=np.append(stname,str(str(tar.getmembers()[t].name).split('/')[-1]).split('/')[-1].split('.')[-2]) @@ -1460,10 +1460,10 @@ def bull_tar(*args): atp[t,:]=np.array(auxtp) else: print(" Time duration of "+tar.getmembers()[t]+" (in "+fname+") do not match the other stations. Mantained NaN.") - + del auxhs,auxtp,tfile,lines - # build dictionary + # build dictionary result={'time':np.array(at).astype('double'),'date':np.array(adate).astype('double'), 'latitude':np.array(alat),'longitude':np.array(alon),'station_name':np.array(stname), 'hs':np.array(ahs),'tp':np.array(atp)} @@ -1494,7 +1494,7 @@ def bull_tar(*args): else: alon=np.append(alon,-1.*float(auxpos[7:13])) - if t==0: + if t==0: # time array ---- auxdate = str(lines[2]).split(':')[1].split('UTC')[0][1::] auxt = np.double(timegm( strptime( auxdate[0:8]+' '+auxdate[9:11]+'00', '%Y%m%d %H%M') )) @@ -1527,21 +1527,21 @@ def bull_tar(*args): if len(auxlines[10:15].replace(' ',''))>0: auxhs=np.append(auxhs,float(auxlines[10:15])) fuxhs=[] - for k in range(0,4): + for k in range(0,4): if len(str(auxlines[int(iauxhs[0]+18*k):int(iauxhs[1]+18*k)]).replace(' ', '')): fuxhs=np.append(fuxhs,float(auxlines[int(iauxhs[0]+18*k):int(iauxhs[1]+18*k)])) else: fuxhs=np.append(fuxhs,np.nan) fuxtp=[] - for k in range(0,4): + for k in range(0,4): if len(str(auxlines[int(iauxtp[0]+18*k):int(iauxtp[1]+18*k)]).replace(' ', '')): fuxtp=np.append(fuxtp,float(auxlines[int(iauxtp[0]+18*k):int(iauxtp[1]+18*k)])) else: fuxtp=np.append(fuxtp,np.nan) fuxdp=[] - for k in range(0,4): + for k in range(0,4): if len(str(auxlines[int(iauxdp[0]+18*k):int(iauxdp[1]+18*k)]).replace(' ', '')): fuxdp=np.append(fuxdp,float(auxlines[int(iauxdp[0]+18*k):int(iauxdp[1]+18*k)])) else: @@ -1569,7 +1569,7 @@ def bull_tar(*args): del auxhs,auxtp,auxdp,tfile,lines - # build dictionary + # build dictionary result={'time':np.array(at).astype('double'),'date':np.array(adate).astype('double'), 'latitude':np.array(alat),'longitude':np.array(alon),'station_name':np.array(stname), 'hs':np.array(ahs),'tp':np.array(atp),'dp':np.array(adp)} @@ -1621,7 +1621,7 @@ def ts(*args): ahspr=np.append(ahspr,np.nan) atp=np.append(atp,np.nan) - # build dictionary + # build dictionary result={'time':np.array(at).astype('double'),'date':np.array(adate).astype('double'), 'station_name':np.array(stname),'hs':np.array(ahs),'hs_spr':np.array(ahspr),'tp':np.array(atp)} @@ -1657,7 +1657,7 @@ def ts(*args): atp[atp<0.01]=np.nan; atp=1./atp - # build dictionary + # build dictionary result={'time':np.array(at).astype('double'),'date':np.array(adate).astype('double'), 'station_name':np.array(stname), 'hs':np.array(ahs),'l':np.array(al), @@ -1738,7 +1738,7 @@ def station_tar(*args): del auxhs,auxhspr,auxtp,tfile,lines - # build dictionary + # build dictionary result={'time':np.array(at).astype('double'),'date':np.array(adate).astype('double'), 'station_name':np.array(stname),'hs':np.array(ahs),'hs_spr':np.array(ahspr),'tp':np.array(atp)} @@ -1784,7 +1784,7 @@ def spec_ndbc(*args): freq = ds['frequency'].values[:] pspec = ds['spectral_wave_density'].values[::sk,:,0,0] dmspec = ds['mean_wave_dir'][::sk,:,0,0] - dpspec = ds['principal_wave_dir'][::sk,:,0,0] + dpspec = ds['principal_wave_dir'][::sk,:,0,0] r1spec = ds['wave_spectrum_r1'][::sk,:,0,0] r2spec = ds['wave_spectrum_r2'][::sk,:,0,0] ds.close(); del ds @@ -1801,10 +1801,10 @@ def spec_ndbc(*args): # final directional wave spectrum (frequency X direction) dirspec = np.zeros((btime.shape[0],freq.shape[0],theta.shape[0]),'f') for t in range(0,btime.shape[0]): - dirspec[t,:,:] = np.array([pspec[t,:]]).T * (1/pi)*(0.5+ np.array([r1spec[t,:]]).T * cos(np.array( np.array([theta])-np.array([dmspec[t,:]]).T )*(pi/180)) + dirspec[t,:,:] = np.array([pspec[t,:]]).T * (1/pi)*(0.5+ np.array([r1spec[t,:]]).T * cos(np.array( np.array([theta])-np.array([dmspec[t,:]]).T )*(pi/180)) + np.array([r2spec[t,:]]).T*cos(2*np.array( np.array([theta]) - np.array([dpspec[t,:]]).T )*(pi/180))) - # build dictionary + # build dictionary result={'time':btime,'date':bdate,'latitude':blat,'longitude':blon, 'freq':freq,'deltafreq':dfreq,'pspec':pspec,'dmspec':dmspec,'dpspec':dpspec, 'theta':theta,'dirspec':dirspec} @@ -1874,8 +1874,8 @@ def spec_ww3(*args): # water depth (constant in time) depth=np.nanmean(ds['dpt'].values[::sk,inds],axis=0) lon=np.array(np.nanmean(ds['longitude'].values[::sk,inds],axis=0)) - lat=np.array(np.nanmean(ds['latitude'].values[::sk,inds],axis=0)) - + lat=np.array(np.nanmean(ds['latitude'].values[::sk,inds],axis=0)) + ds.close(); del ds, auxstationname, inds, stationname freq1=freq; freq2=freq @@ -1889,7 +1889,7 @@ def spec_ww3(*args): cabc=fp.readline(); cabc=cabc.strip().split() nf=int(cabc[3]) # number of frequencies nd=int(cabc[4]) # number of directions - npo=int(cabc[5]) # number of point outputs + npo=int(cabc[5]) # number of point outputs freq=zeros(nf,'f');dire=zeros(nd,'f') dspec=zeros((nt,nf,nd),'f') @@ -1912,7 +1912,7 @@ def spec_ww3(*args): line=line.strip().split() for i in range(0,rncf): freq[k]=float(line[i]) - k=k+1 + k=k+1 # DF in frequency (dfreq) dfreq=np.zeros(freq.shape[0],'f') @@ -1944,7 +1944,7 @@ def spec_ww3(*args): wnds=np.zeros((nt),'f');wndd=np.zeros((nt),'f') for t in range(0,nt): - + cabc=fp.readline(); cabc.strip().split()[0] mtime[t] = np.double(timegm( strptime(cabc.strip().split()[0]+cabc.strip().split()[1][0:2], '%Y%m%d%H') )) cabc=fp.readline(); cabc=cabc.strip().split() @@ -1997,7 +1997,7 @@ def spec_ww3(*args): adspec[t,:,nd-(inddire+1):nd]=dspec[t,:,0:(inddire+1)] for i in range(0,nd): dspec[t,:,i]=adspec[t,:,nd-i-1] - + adspec[t,:,0:int(nd/2)]=dspec[t,:,int(nd/2):nd] adspec[t,:,int(nd/2):nd]=dspec[t,:,0:int(nd/2)] dspec[t,:,:]=adspec[t,:,:] @@ -2007,7 +2007,7 @@ def spec_ww3(*args): # 1D directional spectrum d1sp=np.zeros((dspec.shape[0],nf),'f') for t in range(0,dspec.shape[0]): - for il in range(0,nf): + for il in range(0,nf): a = np.sum(dspec[t,il,:] * np.array(np.sin((pi*dire)/180.)/np.sum(dspec[t,il,:])) ) b = np.sum(dspec[t,il,:] * np.array(np.cos((pi*dire)/180.)/np.sum(dspec[t,il,:])) ) aux = math.atan2(a,b)*(180./pi) From e1ec414988d5d8c3cf0dcf7512f3750204e1747b Mon Sep 17 00:00:00 2001 From: Ming Chen Date: Sun, 17 May 2026 17:41:51 +0000 Subject: [PATCH 6/7] add spec_tar.gz load function in wread --- ww3tools/wread.py | 133 +++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 131 insertions(+), 2 deletions(-) diff --git a/ww3tools/wread.py b/ww3tools/wread.py index 0f79ad7..21f25a3 100755 --- a/ww3tools/wread.py +++ b/ww3tools/wread.py @@ -47,6 +47,7 @@ station_tar spec_ndbc spec_ww3 + spec_tar Explanation for each function is contained in the headers OUTPUT: @@ -65,9 +66,10 @@ 03/07/2025: Ricardo M. Campos, new functions included: tseriestxt_ndbc, tseriesnc_cdip, tseriesnc_microswift, tseries_spotter, tseriesnc_dwsd, tseriesnc_saildrone, tseriesnc_wsra, tseriestxt_ww3. New satellite missions added to AODN altimeter data reading - + 05/15/2026: Ming Chen, new functions, spec_tar to read spec_tar.gz operational point output PERSON OF CONTACT: Ricardo M Campos: ricardo.campos@noaa.gov + Ming Chen: ming.chen1@noaa.gov """ @@ -90,7 +92,9 @@ # import pickle import sys import warnings; warnings.filterwarnings("ignore") - +import gzip +import io +import tarfile def readconfig(fname): """ @@ -2025,4 +2029,129 @@ def spec_ww3(*args): return result del mtime,mdate,lat,lon,wnds,wndd,freq,freq1,freq2,dfreq,pwst,dire,d1sp,dspec +# WAVEWATCH III spectra output for wind speed and direction +def spec_tar(*args): + ''' + WAVEWATCH III, spec_tar.gz operational point output. + Input: file name (example: gfswave.t00z.spec_tar.gz) + Output: dictionary containing: + time(seconds since 1970),lat,lon,station names; Arrays: wind_spd, wind_dir + ''' + + if len(args) == 1: + fname = str(args[0]) + else: + sys.exit(' One input is required: spec_tar.gz file name') + + print(" reading ww3 spec_tar.gz file for wind ...") + + # Open gzip -> tar + try: + with gzip.open(fname, "rb") as gz: + tar_bytes = gz.read() + + tar = tarfile.open(fileobj=io.BytesIO(tar_bytes), mode="r:") + + except: + sys.exit(' Cannot open ' + fname) + + members = [m for m in tar.getmembers() if m.isfile()] + + station_name_all = [] + wsp_all = [] + wdir_all = [] + time_all = None + + # Loop station files + for member in members: + + try: + tfile = tar.extractfile(member) + lines = tfile.readlines() + except: + print(" Cannot read " + member.name + ". Skipped.") + continue + + at = [] + awsp = [] + awdir = [] + stname = None + + # Parse file + for j in range(len(lines) - 1): + + line = lines[j].decode("utf-8", errors="ignore").strip() + parts = line.split() + + # Time line follows the pattern: YYYYMMDD HHMMSS + if ( + len(parts) == 2 + and len(parts[0]) == 8 + and len(parts[1]) == 6 + and parts[0].isdigit() + and parts[1].isdigit() + ): + + info_line = lines[j+1].decode( + "utf-8", + errors="ignore" + ).strip() + + info = info_line.replace("'", "").split() + + # Expected line contains: + # station lat lon depth wspd wdir current cdir + if len(info) >= 6: + + try: + tsec = np.double( + timegm( + strptime( + parts[0] + parts[1], + '%Y%m%d%H%M%S' + ) + ) + ) + + stname = str(info[0]) + + at.append(tsec) + awsp.append(float(info[4])) + awdir.append(float(info[5])) + + except: + continue + + if len(at) == 0: + print(" No wind records found in " + member.name) + continue + + station_name_all.append(stname) + wsp_all.append(awsp) + wdir_all.append(awdir) + + # Use first station as reference time + if time_all is None: + time_all = np.array(at).astype('double') + + tar.close() + + # Convert to arrays + wind_spd = np.array(wsp_all).astype('float') + wind_dir = np.array(wdir_all).astype('float') + + # Basic QC + wind_spd[(wind_spd < 0.0) | (wind_spd > 100.0)] = np.nan + wind_dir[(wind_dir < 0.0) | (wind_dir > 360.0)] = np.nan + + result = { + 'station_name': np.array(station_name_all).astype('str'), + 'time': np.array(time_all).astype('double'), + 'wind_spd': wind_spd, + 'wind_dir': wind_dir + } + + print(" ww3 spec_tar.gz wind file OK. " + fname) + + return result From bef45db77db5c0f3dcd8d5af4d352b83396fd748 Mon Sep 17 00:00:00 2001 From: Ming Chen Date: Mon, 18 May 2026 17:51:18 +0000 Subject: [PATCH 7/7] Add model wind speed and direction from spec_tar.gz to buoy collocation output --- ww3tools/modelBuoy_collocation.py | 108 ++++++++++++++++++++++++++++-- 1 file changed, 104 insertions(+), 4 deletions(-) diff --git a/ww3tools/modelBuoy_collocation.py b/ww3tools/modelBuoy_collocation.py index 6ef08be..10a5f7d 100755 --- a/ww3tools/modelBuoy_collocation.py +++ b/ww3tools/modelBuoy_collocation.py @@ -87,7 +87,8 @@ dimensions), and check if variable names exist in the netcdf file (buoy and ww3) to maximize the amount of matchups even when one variable is not available. - + 05/18/2026: Ming Chen, add wind speed and direction loading and writing when + using spec_tar.gz (GFS) PERSON OF CONTACT: Ricardo M Campos: ricardo.campos@noaa.gov @@ -212,8 +213,20 @@ def parse_args(): for i in range(0,np.size(wlist)): if str(wlist[i]).split('/')[-1].split('.')[-1]=='bull_tar': result = wread.bull_tar(wlist[i]) + + # Read model wind from matching spec_tar.gz file + # assume spec_tar.gz files in the same folder of bull_tar (valid for GFS) + specfile = str(wlist[i]).replace('.bull_tar', '.spec_tar.gz') + # check if spec_tar.gz files exist + try: + specresult = wread.spec_tar(specfile) + except: + print(" Warning: Cannot read "+specfile+". Model wind set to NaN.") + specresult = None + if str(wlist[i]).split('/')[-1].split('.')[-1]=='station_tar': result = wread.bull_tar(wlist[i]) + specresult = None # will implement wind when using station_tar if necessary at=result['time'] fcycle = np.array(np.zeros((at.shape[0]),'d')+at[0]).astype('double') @@ -223,6 +236,37 @@ def parse_args(): mfcycle=np.copy(fcycle) mhs=np.copy(result['hs']) mtp=np.copy(result['tp']) + + # allocate wind parameters if wind parameters exist + if specresult is not None: + # reorder specresult to match bull_tar stations + bull_stname = np.array([str(s).strip() for s in result['station_name']]) + spec_stname = np.array([str(s).strip() for s in specresult['station_name']]) + spec_order = [] + for s in bull_stname: + ind = np.where(spec_stname == s)[0] + if np.size(ind) > 0: + spec_order.append(ind[0]) + else: + spec_order.append(-1) + + spec_order = np.array(spec_order) + + # before allocation, check if station and time are matched + if ( + np.array_equal( + np.array([str(s).strip() for s in stname]), + np.array([str(s).strip() for s in specresult['station_name'][spec_order]]) + ) + and np.array_equal(at,specresult['time']) + ): + mwsp=np.copy(specresult['wind_spd'][spec_order, :]) + mwdir=np.copy(specresult['wind_dir'][spec_order, :]) + else: + print(" Warning: spec_tar station/time does not match bull_tar. Model wind set to NaN for "+str(wlist[i])) + mwsp=np.copy(mhs)*np.nan + mwdir=np.copy(mhs)*np.nan + if 'dp' in result.keys(): mdp=np.copy(result['dp']) else: @@ -235,6 +279,36 @@ def parse_args(): mfcycle=np.append(mfcycle,fcycle) mhs=np.append(mhs,result['hs'],axis=1) mtp=np.append(mtp,result['tp'],axis=1) + + if specresult is not None: + bull_stname = np.array([str(s).strip() for s in result['station_name']]) + spec_stname = np.array([str(s).strip() for s in specresult['station_name']]) + spec_order = [] + for s in bull_stname: + ind = np.where(spec_stname == s)[0] + if np.size(ind) > 0: + spec_order.append(ind[0]) + else: + spec_order.append(-1) + spec_order = np.array(spec_order) + + if ( + np.array_equal( + np.array([str(s).strip() for s in stname]), + np.array([str(s).strip() for s in specresult['station_name'][spec_order]]) + ) + and np.array_equal(at,specresult['time']) + ): + mwsp=np.append(mwsp,specresult['wind_spd'][spec_order, :],axis=1) + mwdir=np.append(mwdir,specresult['wind_dir'][spec_order, :],axis=1) + else: + print(" Warning: spec_tar station/time does not match bull_tar. Model wind set to NaN for "+str(wlist[i])) + mwsp=np.append(mwsp,np.copy(result['hs'])*np.nan,axis=1) + mwdir=np.append(mwdir,np.copy(result['hs'])*np.nan,axis=1) + else: + mwsp=np.append(mwsp,np.copy(result['hs'])*np.nan,axis=1) + mwdir=np.append(mwdir,np.copy(result['hs'])*np.nan,axis=1) + if 'dp' in result.keys(): mdp=np.append(mdp,result['dp'],axis=1) else: @@ -243,7 +317,7 @@ def parse_args(): else: print(" Stations in "+wlist[i]+" do not match the other tar files. Skipped "+wlist[i]) - del result,at,fcycle + del result,at,fcycle,specresult,spec_order,bull_stname,spec_stname mdm=np.copy(mhs)*np.nan; mtm=np.copy(mhs)*np.nan # not saved in this file format print(" ww3 file "+wlist[i]+" OK") @@ -670,6 +744,14 @@ def parse_args(): if np.size(ind)>0: mdp[ind]=np.nan; del ind +ind=np.where((mwsp>100.)|(mwsp<0.0)) +if np.size(ind)>0: + mwsp[ind]=np.nan; del ind + +ind=np.where((mwdir>360.)|(mwdir<-180.)) +if np.size(ind)>0: + mwdir[ind]=np.nan; del ind + # Clean data excluding some stations. Select matchups only when model and buoy are available. ind=np.where( (np.isnan(lat)==False) & (np.isnan(lon)==False) & (np.isnan(np.nanmean(mhs,axis=1))==False) & (np.isnan(np.nanmean(bhs,axis=1))==False) ) if np.size(ind)>0: @@ -681,6 +763,8 @@ def parse_args(): mtp=np.array(mtp[ind[0],:]) mdm=np.array(mdm[ind[0],:]) mdp=np.array(mdp[ind[0],:]) + mwsp=np.array(mwsp[ind[0],:]) + mwdir=np.array(mwdir[ind[0],:]) bhs=np.array(bhs[ind[0],:]) btm=np.array(btm[ind[0],:]) btp=np.array(btp[ind[0],:]) @@ -727,6 +811,8 @@ def parse_args(): mtp=np.array(mtp[ind[0],:]) mdm=np.array(mdm[ind[0],:]) mdp=np.array(mdp[ind[0],:]) + mwsp=np.array(mwsp[ind[0],:]) + mwdir=np.array(mwdir[ind[0],:]) bhs=np.array(bhs[ind[0],:]) btm=np.array(btm[ind[0],:]) btp=np.array(btp[ind[0],:]) @@ -779,6 +865,8 @@ def parse_args(): nmtp=np.zeros((mhs.shape[0],unt.shape[0],mxsz),'f')*np.nan nmdm=np.zeros((mhs.shape[0],unt.shape[0],mxsz),'f')*np.nan nmdp=np.zeros((mhs.shape[0],unt.shape[0],mxsz),'f')*np.nan + nmwsp=np.zeros((mhs.shape[0],unt.shape[0],mxsz),'f')*np.nan + nmwdir=np.zeros((mhs.shape[0],unt.shape[0],mxsz),'f')*np.nan nbhs=np.zeros((mhs.shape[0],unt.shape[0],mxsz),'f')*np.nan nbtm=np.zeros((mhs.shape[0],unt.shape[0],mxsz),'f')*np.nan nbtp=np.zeros((mhs.shape[0],unt.shape[0],mxsz),'f')*np.nan @@ -796,6 +884,8 @@ def parse_args(): nmtp[:,i,:][:,0:np.size(ind)]=np.array(mtp[:,ind]) nmdm[:,i,:][:,0:np.size(ind)]=np.array(mdm[:,ind]) nmdp[:,i,:][:,0:np.size(ind)]=np.array(mdp[:,ind]) + nmwsp[:,i,:][:,0:np.size(ind)]=np.array(mwsp[:,ind]) + nmwdir[:,i,:][:,0:np.size(ind)]=np.array(mwdir[:,ind]) nbhs[:,i,:][:,0:np.size(ind)]=np.array(bhs[:,ind]) nbtm[:,i,:][:,0:np.size(ind)]=np.array(btm[:,ind]) nbtp[:,i,:][:,0:np.size(ind)]=np.array(btp[:,ind]) @@ -843,6 +933,8 @@ def parse_args(): vmtp = ncfile.createVariable('model_tp',np.dtype('float32').char,('buoypoints','fcycle','time')) vmdm = ncfile.createVariable('model_dm',np.dtype('float32').char,('buoypoints','fcycle','time')) vmdp = ncfile.createVariable('model_dp',np.dtype('float32').char,('buoypoints','fcycle','time')) + vmwsp = ncfile.createVariable('model_wsp', np.dtype('float32').char, ('buoypoints','fcycle','time')) + vmwdir = ncfile.createVariable('model_wdir', np.dtype('float32').char, ('buoypoints','fcycle','time')) vbhs = ncfile.createVariable('obs_hs',np.dtype('float32').char,('buoypoints','fcycle','time')) vbtm = ncfile.createVariable('obs_tm',np.dtype('float32').char,('buoypoints','fcycle','time')) vbtp = ncfile.createVariable('obs_tp',np.dtype('float32').char,('buoypoints','fcycle','time')) @@ -858,6 +950,8 @@ def parse_args(): vmtp = ncfile.createVariable('model_tp',np.dtype('float32').char,('buoypoints','time')) vmdm = ncfile.createVariable('model_dm',np.dtype('float32').char,('buoypoints','time')) vmdp = ncfile.createVariable('model_dp',np.dtype('float32').char,('buoypoints','time')) + vmwsp = ncfile.createVariable('model_wsp', np.dtype('float32').char, ('buoypoints','time')) + vmwdir = ncfile.createVariable('model_wdir', np.dtype('float32').char, ('buoypoints','time')) vbhs = ncfile.createVariable('obs_hs',np.dtype('float32').char,('buoypoints','time')) vbtm = ncfile.createVariable('obs_tm',np.dtype('float32').char,('buoypoints','time')) vbtp = ncfile.createVariable('obs_tp',np.dtype('float32').char,('buoypoints','time')) @@ -887,8 +981,8 @@ def parse_args(): vmtp.units='s'; vbtp.units='s' vmdm.units='degrees'; vbdm.units='degrees' vmdp.units='degrees'; vbdp.units='degrees' - vbwsp.units = 'm s-1' - vbwdir.units = 'degrees' + vmwsp.units = 'm s-1'; vbwsp.units = 'm s-1' + vmwdir.units = 'degrees'; vbwdir.units = 'degrees' if gridinfo!=0: vpdepth.units='m'; vpdistcoast.units='km' @@ -905,6 +999,8 @@ def parse_args(): vmtp.long_name = 'Model Peak Wave Period' vmdm.long_name = 'Model Mean Wave Direction' vmdp.long_name = 'Model Peak Wave Direction' + vmwsp.long_name = 'Model Wind Speed' + vmwdir.long_name = 'Model Wind Direction' vbhs.long_name = 'Observed Significant Wave Height' vbtm.long_name = 'Observed Mean Wave Period' @@ -923,6 +1019,8 @@ def parse_args(): vmtp[:,:,:]=nmtp[:,:,:] vmdm[:,:,:]=nmdm[:,:,:] vmdp[:,:,:]=nmdp[:,:,:] + vmwsp[:,:,:]=nmwsp[:,:,:] + vmwdir[:,:,:]=nmwdir[:,:,:] vbhs[:,:,:]=nbhs[:,:,:] vbtm[:,:,:]=nbtm[:,:,:] vbtp[:,:,:]=nbtp[:,:,:] @@ -937,6 +1035,8 @@ def parse_args(): vmtp[:,:]=mtp[:,:] vmdm[:,:]=mdm[:,:] vmdp[:,:]=mdp[:,:] + vmwsp[:,:]=mwsp[:,:] + vmwdir[:,:]=mwdir[:,:] vbhs[:,:]=bhs[:,:] vbtm[:,:]=btm[:,:] vbtp[:,:]=btp[:,:]