diff --git a/ww3tools/modelBuoy_collocation.py b/ww3tools/modelBuoy_collocation.py old mode 100755 new mode 100644 index 3002b46..ff9fb8a --- a/ww3tools/modelBuoy_collocation.py +++ b/ww3tools/modelBuoy_collocation.py @@ -11,6 +11,7 @@ v1.3 11/17/2022 v1.4 12/08/2022 v1.5 01/31/2023 + v1.6 05/21/2024 PURPOSE: Collocation/pairing ww3 point output results with wave buoys. @@ -33,7 +34,7 @@ USAGE: Input WW3 point outputs are utilized (not grids), with option of - reading different formats: netcdf, text (bull and ts), +modelBuoy_collocation.py 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 @@ -55,6 +56,37 @@ multiple forecast data files: nohup python3 modelBuoy_collocation.py ww3list.gfs-d36.GSE1.5.txt 2 gridInfo.nc CycloneMap_2021.nc >> nohup_modelBuoy_collocation.out 2>&1 & +- From 05/21/24 if you want to run this you have to define the (buoy_path) variable in the command line or in the job script. This recent does not accept the +hard coded ndbcp path anymore. + +USAGE (spec.gz): + In wread.py the ww3_spec1 is responsible for reading and accepting the spec.gz files like,ex:gfswave.t00z.spec_tar.gz + sort of files. + to run the code one should have to define these values in the job script or the command line: + +# Define the arguments + input_gz_file="./gfswave.t00z.spec_tar.gz" + output_directory="./" (Where the users wants to have the output saved) + buoy_path="/scratch2/NCEPDEV/marine/Matthew.Masarik/dat/buoys/NDBC/ncformat/wparam" (The path should be changed based on the users buoy path directory) + model_name="HR3a" # Define the model name + forecast_ds="1" # Indicator for forecast data structure, set to 1 or 0 based on requirements + (If it is 1 or greater than 1, it does the forecast structure) + +# Process data for each date + python3 modelBuoy_collocation.py spec.gz "$input_gz_file" "$output_directory" "$buoy_path" "$model_name " "$forecast_ds" + + +Note: Updated way for submitting formats(bulltar) other than spec.gz: + + +buoy_path="/scratch2/NCEPDEV/marine/Matthew.Masarik/dat/buoys/NDBC/ncformat/wparam" +The buoy path should be defined in the job script. + +# Process data for each date +python3 modelBuoy_collocation.py ww3list.txt(any text file name that user defined) 2 $buoy_path +(The buoy path should be the last item) + + OUTPUT: netcdf file WW3.Buoy*.nc containing matchups of buoy and ww3 data, for the stations (lat/lon) where both data sources are available. @@ -87,748 +119,1058 @@ 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/21/2024 : Ghazal Mohammadpour , developed the code to be able to process + the spec.gz files. + PERSON OF CONTACT: Ricardo M Campos: ricardo.campos@noaa.gov """ -import warnings; warnings.filterwarnings("ignore") +import sys +import os +import gzip +import tarfile +import warnings import numpy as np -from matplotlib.mlab import * -from pylab import * -import xarray as xr import netCDF4 as nc import time from time import strptime from calendar import timegm import wread -# 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(' ') - -# 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') - -# READ DATA -print(" ") -if gridinfo!=0: - # Grid Information - gridmask = wread.mask(gridinfo) - 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'] - print(" GridInfo Ok. "+gridinfo) - - # Cyclone Information - if cyclonemap!=0: - cycloneinfo = wread.cyclonemap(cyclonemap) - 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: - print(" CycloneMap Ok. "+cyclonemap) - else: - sys.exit(' Error: Cyclone grid and Mask grid are different.') - -# MODEL Point output, search among possible formats ------------------ - -if (str(wlist[0]).split('/')[-1].split('.')[-1]=='bull_tar') or (str(wlist[0]).split('/')[-1].split('.')[-1]=='station_tar'): - for i in range(0,np.size(wlist)): - if str(wlist[i]).split('/')[-1].split('.')[-1]=='bull_tar': - result = wread.bull_tar(wlist[i]) - if str(wlist[i]).split('/')[-1].split('.')[-1]=='station_tar': - result = wread.bull_tar(wlist[i]) - - at=result['time'] - fcycle = np.array(np.zeros((at.shape[0]),'d')+at[0]).astype('double') - if i==0: - stname=result['station_name'] - mtime=np.copy(at) - mfcycle=np.copy(fcycle) - mhs=np.copy(result['hs']) - mtp=np.copy(result['tp']) - if 'dp' in result.keys(): - mdp=np.copy(result['dp']) - else: - mdp=np.copy(mhs)*np.nan - - else: - 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) - mhs=np.append(mhs,result['hs'],axis=1) - mtp=np.append(mtp,result['tp'],axis=1) - if 'dp' in result.keys(): - mdp=np.append(mdp,result['dp'],axis=1) - else: - mdp=np.append(mdp,np.copy(result['hs'])*np.nan,axis=1) - - else: - print(" Stations in "+wlist[i]+" do not match the other tar files. Skipped "+wlist[i]) - - del result,at,fcycle - mdm=np.copy(mhs)*np.nan; mtm=np.copy(mhs)*np.nan # not saved in this file format - print(" ww3 file "+wlist[i]+" OK") - -elif (str(wlist[0]).split('/')[-1].split('.')[-1]=='bull') or (str(wlist[0]).split('/')[-1].split('.')[-1]=='ts'): - if forecastds>0: - forecastds=0 # no 2D time array possible. - print(" Warning: no 2D time array possible. \ - Use multiple tar files (one per forecast cycle), station_tar or bull_tar, containing text files for each station.") - - for i in range(0,np.size(wlist)): - # re-check each file in the list respects the same format - if str(wlist[i]).split('/')[-1].split('.')[-1]=='ts': - result = wread.ts(wlist[i]) - elif str(wlist[i]).split('/')[-1].split('.')[-1]=='bull': - result = wread.bull(wlist[i]) - - at=result['time'] - if i==0: - mfcycle = np.array(np.zeros((at.shape[0]),'d')+at[0]).astype('double') - stname=np.atleast_1d(np.array(result['station_name'])) - mtime=np.copy(at) - mhs=np.copy([result['hs']]) - mtp=np.copy([result['tp']]) - if 'dp' in result.keys(): - mdp=np.copy([result['dp']]) - else: - mdp=np.copy(mhs)*np.nan - - if 'dm' in result.keys(): - mdm=np.copy([result['dm']]) - else: - mdm=np.copy(mhs)*np.nan - - if 'tm' in result.keys(): - mtm=np.copy([result['tm']]) - else: - mtm=np.copy(mhs)*np.nan - - else: - if (mhs.shape[1]==result['hs'].shape[0]): - stname=np.append(stname,np.atleast_1d(np.array(result['station_name']))) - mtime=np.append(mtime,at) - mhs=np.append(mhs,[result['hs']],axis=0) - mtp=np.append(mtp,[result['tp']],axis=0) - if 'dp' in result.keys(): - mdp=np.append(mdp,[result['dp']],axis=0) - else: - mdp=np.append(mdp,[np.copy(result['hs'])*np.nan],axis=0) - - if 'dm' in result.keys(): - mdm=np.append(mdp,[result['dm']],axis=0) - else: - mdm=np.append(mdm,[np.copy(result['hs'])*np.nan],axis=0) - - if 'tm' in result.keys(): - mtm=np.append(mtm,[result['tm']],axis=0) - else: - mtm=np.append(mtm,[np.copy(result['hs'])*np.nan],axis=0) - - else: - print(" Stations in "+wlist[i]+" do not match the other tar files. Skipped "+wlist[i]) - - del result,at - print(" ww3 file "+wlist[i]+" OK") - - -elif str(wlist[0]).split('/')[-1].split('.')[-1]=='nc': - print(" Using ww3 netcdf point output format") - # netcdf point output file - for t in range(0,np.size(wlist)): - try: - f=nc.Dataset(str(wlist[t])) - except: - print(" Cannot open "+wlist[t]) - else: - # list of station/buoy names - if t==0: - auxstationname=f.variables['station_name'][:,:]; stname=[] - for i in range(0,auxstationname.shape[0]): - astname="".join(np.array(auxstationname[i,:]).astype('str')) - if '\t' in astname: - astname=str(astname).replace("\t","") - - stname=np.append(stname,astname); del astname - - funits=f.variables['time'].units - if str(funits).split(' ')[0] == 'seconds': - tincr=1 - elif str(funits).split(' ')[0] == 'hours': - tincr=3600 - elif str(funits).split(' ')[0] == 'days': - tincr=24*3600 - - ahs = np.array(f.variables['hs'][:,:]).T - - if 'th1m' in f.variables.keys(): - adm = np.array(f.variables['th1m'][:,:]).T - else: - adm = np.array(np.copy(ahs*nan)) - - if 'th1p' in f.variables.keys(): - adp = np.array(f.variables['th1p'][:,:]).T - else: - adp = np.array(np.copy(ahs*nan)) - - 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(): - auxtp = np.array(f.variables['fp'][:,:]).T - indtp=np.where(auxtp>0.) - if np.size(indtp)>0: - atp=np.copy(auxtp) - atp[indtp]=np.copy(1./atp[indtp]) - del indtp - - del auxtp - - else: - atm = np.array(np.copy(ahs*nan)) - - ftunits=str(f.variables['time'].units).split('since')[1][1::].replace('T',' ').replace('+00:00','') - at = np.array(f.variables['time'][:]*tincr + timegm( strptime(ftunits,'%Y-%m-%d %H:%M:%S') )).astype('double') - - f.close(); del f - fcycle = np.array(np.zeros((at.shape[0]),'d')+at[0]).astype('double') - if t==0: - mhs=np.copy(ahs) - mtm=np.copy(atm) - mtp=np.copy(atp) - mdm=np.copy(adm) - mdp=np.copy(adp) - mtime=np.copy(at) - mfcycle=np.copy(fcycle) - else: - mhs=np.append(mhs,ahs,axis=1) - mtm=np.append(mtm,atm,axis=1) - mtp=np.append(mtp,atp,axis=1) - mdm=np.append(mdm,adm,axis=1) - mdp=np.append(mdp,adp,axis=1) - mtime=np.append(mtime,at) - mfcycle=np.append(mfcycle,fcycle) - - del ahs,atm,atp,adm,adp,at,fcycle - - print(" Read WW3 data OK."); print(' ') +import glob +import shutil + +# Suppressing warnings +warnings.filterwarnings("ignore") + +# netCDF format specification +fnetcdf = "NETCDF4" + +gridinfo = None +cyclonemap = None +wlist = [] +ftag = '' + +# Function to unzip and untar files +def unzip_and_untar(gz_file, output_dir): + try: + file_name = os.path.splitext(os.path.basename(gz_file))[0] + base_name = file_name.rsplit('.', 1)[0] + output_file = os.path.join(output_dir, file_name) + + with gzip.open(gz_file, 'rb') as f_in: + with open(output_file, 'wb') as f_out: + f_out.write(f_in.read()) + + extracted_folder = os.path.join(output_dir, base_name) + + with tarfile.open(output_file, 'r') as tar: + tar.extractall(extracted_folder) + + # Creating a list of the extracted files + list_file = os.path.join(output_dir, f'{base_name}_contents.txt') + with open(list_file, 'w') as f: + for root, dirs, files in os.walk(extracted_folder): + for file in files: + f.write(os.path.join(root, file) + '\n') + + # Creating a list of ids between two dots in the filenames + id_file = os.path.join(output_dir, f'{base_name}_id.txt') + with open(id_file, 'w') as f: + for root, dirs, files in os.walk(extracted_folder): + for file in files: + file_parts = file.split('.') + if len(file_parts) >= 3: + id_between_dots = file_parts[1] + f.write(id_between_dots + '\n') + + # Return the path to the extracted folder + return extracted_folder, list_file, id_file + except Exception as e: + print(f"Error in unzip_and_untar: {e}") + sys.exit(1) + +# Main part of the script +if __name__ == "__main__": + mode = sys.argv[1] if len(sys.argv) > 1 else "" + + if mode == "spec.gz": + if len(sys.argv) != 7: + sys.exit("Usage: python script.py spec.gz input_gz_file output_directory buoy_path model_name forecastds") + + gz_file = sys.argv[2] + output_dir = sys.argv[3] + buoy_path = sys.argv[4] + model_name = sys.argv[5] + forecastds = int(sys.argv[6]) # Convert forecast indicator to integer + + # Debugging lines + print(f"gz_file: {gz_file}") + print(f"output_dir: {output_dir}") + print(f"buoy_path: {buoy_path}") + print(f"model_name: {model_name}") + print(f"forecastds: {forecastds}") + + print(f"Unzipping and untarring {gz_file} to {output_dir}") + ndbcp = buoy_path + extracted_folder, list_file, id_file = unzip_and_untar(gz_file, output_dir) + + print("Extraction completed. Extracted files are in:", extracted_folder) + print("Contents list file:", list_file) + print("ID list file:", id_file) + + # Read file names from the extracted folder + file_names = [] + for root, dirs, files in os.walk(extracted_folder): + for file in files: + file_names.append(os.path.join(root, file)) + + if not file_names: + print("No files found in the extracted folder.") + sys.exit(1) + + print("List of files to process:", file_names) + + # Read id_files (if needed) + id_files = [] + with open(id_file, 'r') as f: + for line in f: + id_files.append(line.strip()) + + if not id_files: + print("No IDs found in the ID list file.") + sys.exit(1) + + print("List of IDs:", id_files) + + # Call the spec_ww3 function with the file names and ID files + try: + print("Calling wread.spec1_ww3...") + results =wread. spec1_ww3(file_names, id_files) + print("wread.spec1_ww3 call completed.") + except Exception as e: + print(f"Error in wread.spec1_ww3: {e}") + sys.exit(1) + + # Additional processing (if needed) + mfcycle = None + stname = None + mtime = None + mhs = None + mtp = None + mwn = None + mwd = None + mfreq = None + mtm = None + mdm = None + mdp = None + + try: + for i, result in enumerate(results): + print(f"Processing result {i+1}/{len(results)}...") + at = result['time'] + + if i == 0: + mfcycle = np.array(np.zeros((at.shape[0]), 'd') + at[0]).astype('double') + stname = np.atleast_1d(np.array(result['station_name'])) + mtime = np.copy(at) + mhs = np.copy([result['Hs']]) + mtp = np.copy([result['Tp']]) + + if 'wind_spd' in result.keys(): + mwn = np.copy([result['wind_spd']]) + else: + mwn = np.copy(mhs) * np.nan + + if 'wind_dir' in result.keys(): + mwd = np.copy([result['wind_dir']]) + else: + mwd = np.copy(mhs) * np.nan + + if 'freq' in result.keys(): + mfreq = np.copy([result['freq']]) + else: + mfreq = np.copy(mhs) * np.nan + + if 'tm' in result.keys(): + mtm = np.copy([result['tm']]) + else: + mtm = np.copy(mhs) * np.nan + + if 'dp' in result.keys(): + mdp = np.copy([result['dp']]) + else: + mdp = np.copy(mhs) * np.nan + + if 'dm' in result.keys(): + mdm = np.copy([result['dm']]) + else: + mdm = np.copy(mhs) * np.nan + + else: + if mhs.shape[1] == result['Hs'].shape[0]: + stname = np.append(stname, np.atleast_1d(np.array(result['station_name']))) + mhs = np.append(mhs, [result['Hs']], axis=0) + mtp = np.append(mtp, [result['Tp']], axis=0) + + if 'wind_spd' in result.keys(): + mwn = np.append(mwn, [result['wind_spd']], axis=0) + else: + mwn = np.append(mwn, [np.copy(result['Hs']) * np.nan], axis=0) + + if 'wind_dir' in result.keys(): + mwd = np.append(mwd, [result['wind_dir']], axis=0) + else: + mwd = np.append(mwd, [np.copy(result['Hs']) * np.nan], axis=0) + + if 'freq' in result.keys(): + mfreq = np.append(mfreq, [result['freq']], axis=0) + else: + mfreq = np.append(mfreq, [np.copy(result['Hs']) * np.nan], axis=0) + + if 'tm' in result.keys(): + mtm = np.append(mtm, [result['tm']], axis=0) + else: + mtm = np.append(mtm, [np.copy(result['Hs']) * np.nan], axis=0) + + if 'dp' in result.keys(): + mdp = np.append(mdp, [result['dp']], axis=0) + else: + mdp = np.append(mdp, [np.copy(result['Hs']) * np.nan], axis=0) + + if 'dm' in result.keys(): + mdm = np.append(mdp, [result['dm']], axis=0) + else: + mdm = np.append(mdm, [np.copy(result['Hs']) * np.nan], axis=0) + + else: + print(" Stations in file do not match the other tar files. Skipped") + + print("Completed processing results.") + + # Add additional debugging statements here + print("Checking final data structures...") + print(f"mfcycle: {mfcycle}") + print(f"stname: {stname}") + print(f"mtime: {mtime}") + print(f"mhs: {mhs}") + print(f"mtp: {mtp}") + print(f"mwn: {mwn}") + print(f"mwd: {mwd}") + print(f"mfreq: {mfreq}") + print(f"mtm: {mtm}") + print(f"mdm: {mdm}") + print(f"mdp: {mdp}") + + except Exception as e: + print(f"Error during processing: {e}") + sys.exit(1) + + + else: + wlist = [] + ftag = '' + forecastds = 0 + + if len(sys.argv) < 3: + sys.exit('At least one argument (list of ww3 files) must be informed.') + 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) + + forecastds = int(sys.argv[2]) if len(sys.argv) >= 4 else 0 + if forecastds > 0: + print('Forecast-type data structure') + + # Determine if gridinfo and cyclonemap are provided + if len(sys.argv) >= 6: + ndbcp = sys.argv[3] + gridinfo = sys.argv[4] + cyclonemap = sys.argv[5] + print('Using ndbcp path:', ndbcp) + print('Using gridinfo:', gridinfo) + print('Using cyclonemap:', cyclonemap) + else: + # If gridinfo and cyclonemap are not both provided, use the last argument as buoy_path + buoy_path = sys.argv[-1] + ndbcp = buoy_path + # Check if gridinfo is provided without cyclonemap or vice versa + if 'gridinfo' in sys.argv: + gridinfo = sys.argv[3] + print('Using gridinfo:', gridinfo) + + if 'cyclonemap' in sys.argv: + cyclonemap = sys.argv[3] + print('Using cyclonemap:', cyclonemap) + + print('Using ndbcp path:', buoy_path) + + # READ DATA + print(" ") + if gridinfo: + # Grid Information + gridmask = wread.mask(gridinfo) + 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'] + print(" GridInfo Ok. " + gridinfo) + + # Cyclone Information + if cyclonemap: + cycloneinfo = wread.cyclonemap(cyclonemap) + 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) & np.array_equal(clon, mlon): + print(" CycloneMap Ok. " + cyclonemap) + else: + sys.exit(' Error: Cyclone grid and Mask grid are different.') + + # MODEL Point output, search among possible formats + if (str(wlist[0]).split('/')[-1].split('.')[-1] == 'bull_tar') or ( + str(wlist[0]).split('/')[-1].split('.')[-1] == 'station_tar'): + for i in range(0, np.size(wlist)): + if str(wlist[i]).split('/')[-1].split('.')[-1] == 'bull_tar': + result = wread.bull_tar(wlist[i]) + if str(wlist[i]).split('/')[-1].split('.')[-1] == 'station_tar': + result = wread.bull_tar(wlist[i]) + + at = result['time'] + fcycle = np.array(np.zeros((at.shape[0]), 'd') + at[0]).astype('double') + if i == 0: + stname = result['station_name'] + mtime = np.copy(at) + mfcycle = np.copy(fcycle) + mhs = np.copy(result['hs']) + mtp = np.copy(result['tp']) + if 'dp' in result.keys(): + mdp = np.copy(result['dp']) + else: + mdp = np.copy(mhs) * np.nan + + if 'wind_spd' in result.keys(): + mwn = np.copy(result['wind_spd']) + else: + mwn = np.copy(mhs) * np.nan + + if 'wind_dir' in result.keys(): + mwd = np.copy(result['wind_dir']) + else: + mwd = np.copy(mhs) * np.nan + + else: + 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) + mhs = np.append(mhs, result['hs'], axis=1) + mtp = np.append(mtp, result['tp'], axis=1) + if 'dp' in result.keys(): + mdp = np.append(mdp, result['dp'], axis=1) + else: + mdp = np.append(mdp, np.copy(result['hs']) * np.nan, axis=1) + + if 'wind_spd' in result.keys(): + mwn = np.append(mwn, result['wind_spd'], axis=1) + else: + mwn = np.append(mwn, np.copy(result['hs']) * np.nan, axis=1) + + if 'wind_dir' in result.keys(): + mwd = np.append(mwd, result['wind_dir'], axis=1) + else: + mwd = np.append(mwd, np.copy(result['hs']) * np.nan, axis=1) + + else: + print(" Stations in " + wlist[i] + " do not match the other tar files. Skipped " + wlist[i]) + + del result, at, fcycle + mdm = np.copy(mhs) * np.nan + mtm = np.copy(mhs) * np.nan + print(" ww3 file " + wlist[i] + " OK") + + elif (str(wlist[0]).split('/')[-1].split('.')[-1] == 'bull') or ( + str(wlist[0]).split('/')[-1].split('.')[-1] == 'ts'): + if forecastds > 0: + forecastds = 0 # no 2D time array possible. + print( + " Warning: no 2D time array possible. Use multiple tar files (one per forecast cycle), station_tar or bull_tar, containing text files for each station.") + + for i in range(0, np.size(wlist)): + # re-check each file in the list respects the same format + if str(wlist[i]).split('/')[-1].split('.')[-1] == 'ts': + result = wread.ts(wlist[i]) + elif str(wlist[i]).split('/')[-1].split('.')[-1] == 'bull': + result = wread.bull(wlist[i]) + + at = result['time'] + if i == 0: + mfcycle = np.array(np.zeros((at.shape[0]), 'd') + at[0]).astype('double') + stname = np.atleast_1d(np.array(result['station_name'])) + mtime = np.copy(at) + mhs = np.copy([result['hs']]) + mtp = np.copy([result['tp']]) + if 'dp' in result.keys(): + mdp = np.copy([result['dp']]) + else: + mdp = np.copy(mhs) * np.nan + + if 'dm' in result.keys(): + mdm = np.copy([result['dm']]) + else: + mdm = np.copy(mhs) * np.nan + + if 'tm' in result.keys(): + mtm = np.copy([result['tm']]) + else: + mtm = np.copy(mhs) * np.nan + + else: + if (mhs.shape[1] == result['hs'].shape[0]): + stname = np.append(stname, np.atleast_1d(np.array(result['station_name']))) + mtime = np.append(mtime, at) + mhs = np.append(mhs, [result['hs']], axis=0) + mtp = np.append(mtp, [result['tp']], axis=0) + if 'dp' in result.keys(): + mdp = np.append(mdp, [result['dp']], axis=0) + else: + mdp = np.append(mdp, [np.copy(result['hs']) * np.nan], axis=0) + + if 'dm' in result.keys(): + mdm = np.append(mdp, [result['dm']], axis=0) + else: + mdm = np.append(mdm, [np.copy(result['hs']) * np.nan], axis=0) + + if 'tm' in result.keys(): + mtm = np.append(mtm, [result['tm']], axis=0) + else: + mtm = np.append(mtm, [np.copy(result['hs']) * np.nan], axis=0) + + else: + print(" Stations in " + wlist[i] + " do not match the other tar files. Skipped " + wlist[i]) + + del result, at + print(" ww3 file " + wlist[i] + " OK") + + elif str(wlist[0]).split('/')[-1].split('.')[-1] == 'nc': + print(" Using ww3 netcdf point output format") + # netcdf point output file + for t in range(0, np.size(wlist)): + try: + f = nc.Dataset(str(wlist[t])) + except: + print(" Cannot open " + wlist[t]) + else: + # list of station/buoy names + if t == 0: + auxstationname = f.variables['station_name'][:, :] + stname = [] + for i in range(0, auxstationname.shape[0]): + astname = "".join(np.array(auxstationname[i, :]).astype('str')) + if '\t' in astname: + astname = str(astname).replace("\t", "") + + stname = np.append(stname, astname) + del astname + + funits = f.variables['time'].units + if str(funits).split(' ')[0] == 'seconds': + tincr = 1 + elif str(funits).split(' ')[0] == 'hours': + tincr = 3600 + elif str(funits).split(' ')[0] == 'days': + tincr = 24 * 3600 + + ahs = np.array(f.variables['hs'][:, :]).T + + if 'th1m' in f.variables.keys(): + adm = np.array(f.variables['th1m'][:, :]).T + else: + adm = np.array(np.copy(ahs * np.nan)) + + if 'th1p' in f.variables.keys(): + adp = np.array(f.variables['th1p'][:, :]).T + else: + adp = np.array(np.copy(ahs * np.nan)) + + if 'tr' in f.variables.keys(): + atm = np.array(f.variables['tr'][:, :]).T + else: + atm = np.array(np.copy(ahs * np.nan)) + + if 'fp' in f.variables.keys(): + auxtp = np.array(f.variables['fp'][:, :]).T + indtp = np.where(auxtp > 0.) + if np.size(indtp) > 0: + atp = np.copy(auxtp) + atp[indtp] = np.copy(1. / atp[indtp]) + del indtp + + del auxtp + + else: + atm = np.array(np.copy(ahs * np.nan)) + + ftunits = str(f.variables['time'].units).split('since')[1][1::].replace('T', ' ').replace('+00:00', '') + at = np.array( + f.variables['time'][:] * tincr + timegm(strptime(ftunits, '%Y-%m-%d %H:%M:%S'))).astype('double') + + f.close() + del f + fcycle = np.array(np.zeros((at.shape[0]), 'd') + at[0]).astype('double') + if t == 0: + mhs = np.copy(ahs) + mtm = np.copy(atm) + mtp = np.copy(atp) + mdm = np.copy(adm) + mdp = np.copy(adp) + mwn = np.copy(awm) + mtime = np.copy(at) + mfcycle = np.copy(fcycle) + else: + mhs = np.append(mhs, ahs, axis=1) + mtm = np.append(mtm, atm, axis=1) + mtp = np.append(mtp, atp, axis=1) + mdm = np.append(mdm, adm, axis=1) + mdp = np.append(mdp, adp, axis=1) + mwn = np.append(mwn, awm, axis=1) + mtime = np.append(mtime, at) + mfcycle = np.append(mfcycle, fcycle) + + del ahs, atm, atp, adm, adp, at, fcycle + + print(" Read WW3 data OK.") + print(' ') + + else: + sys.exit(' Point output file format not recognized: only bull, bull_tar, and .nc implemented.') + # include other text formats: tab50, .ts + + print(" Start building the matchups model/buoy ...") + print(' ') + + + +# BUOYS +bwind = np.zeros((np.size(stname), np.size(mtime)), 'f') * np.nan +bhs = np.zeros((np.size(stname), np.size(mtime)), 'f') * np.nan +btm = np.zeros((np.size(stname), np.size(mtime)), 'f') * np.nan +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 +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') + +# Loop buoys +for b in range(0, np.size(stname)): + ahs = [] + try: + awm = [] + ahs = [] + atm = [] + atp = [] + adm = [] + atime = [] + for y in yrange: + f = nc.Dataset(ndbcp + "/" + stname[b] + "h" + repr(y) + ".nc") + if 'wave_height' in f.variables.keys(): + 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(): + ahs = np.append(ahs, f.variables['swh'][:, 0, 0]) + + if 'wind_spd' in f.variables.keys(): + awm = np.append(awm, f.variables['wind_spd'][:, 0, 0]) + else: + awm = np.array(np.copy(ahs * np.nan)) + + if 'average_wpd' in f.variables.keys(): + atm = np.append(atm, f.variables['average_wpd'][:, 0, 0]) + else: + atm = np.array(np.copy(ahs * np.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 * np.nan)) + + if 'mean_wave_dir' in f.variables.keys(): + adm = np.append(adm, f.variables['mean_wave_dir'][:, 0, 0]) + else: + adm = np.array(np.copy(ahs * np.nan)) + + if 'latitude' in f.variables.keys(): + lat[b] = f.variables['latitude'][:] + elif 'LATITUDE' in f.variables.keys(): + lat[b] = f.variables['LATITUDE'][:] + else: + lat[b] = np.nan + + if 'longitude' in f.variables.keys(): + lon[b] = f.variables['longitude'][:] + elif 'LONGITUDE' in f.variables.keys(): + lon[b] = f.variables['LONGITUDE'][:] + else: + lon[b] = np.nan + + atime = np.append(atime, np.array(f.variables['time'][:]).astype('double')) + + f.close() + del f + + adp = adm * np.nan # no peak direction available in this format + + if np.size(ahs) > 0: + # First layer of simple quality-control + indq = np.where((ahs > 30.) | (ahs < 0.0)) + if np.size(indq) > 0: + ahs[indq] = np.nan + del indq + + indq = np.where((atm > 40.) | (atm < 0.0)) + if np.size(indq) > 0: + atm[indq] = np.nan + del indq + + indq = np.where((atp > 40.) | (atp < 0.0)) + if np.size(indq) > 0: + atp[indq] = np.nan + del indq + + indq = np.where((adm > 360.) | (adm < -180.)) + if np.size(indq) > 0: + adm[indq] = np.nan + del indq + + indq = np.where((adp > 360.) | (adp < -180.)) + if np.size(indq) > 0: + adp[indq] = np.nan + del indq + + indq = np.where((awm > 50.) | (awm < 0.0)) + if np.size(indq) > 0: + awm[indq] = np.nan + del indq + + c = 0 + for t in mtime: + indt = np.where(np.abs(atime - t) < 1800.)[0] + if np.size(indt) > 0: + if np.any(ahs[indt].mask == False): + bhs[b, c] = np.nanmean(ahs[indt][ahs[indt].mask == False]) + if np.any(atm[indt].mask == False): + btm[b, c] = np.nanmean(atm[indt][atm[indt].mask == False]) + if np.any(atp[indt].mask == False): + btp[b, c] = np.nanmean(atp[indt][atp[indt].mask == False]) + if np.any(adm[indt].mask == False): + bdm[b, c] = np.nanmean(adm[indt][adm[indt].mask == False]) + if np.any(adp[indt].mask == False): + bdp[b, c] = np.nanmean(adp[indt][adp[indt].mask == False]) + if np.any(awm[indt].mask == False): + bwind[b, c] = np.nanmean(awm[indt][awm[indt].mask == False]) + c += 1 + + print(" station " + stname[b] + " ok") + + except Exception as e: + print("Error occurred while processing station", stname[b]) + print(e) -else: - sys.exit(' Point output file format not recognized: only bull, bull_tar, and .nc implemented.') - # include other text formats: tab50, .ts - - -print(" Start building the matchups model/buoy ..."); print(' ') - -# BUOYS ------------------ -bhs=np.zeros((np.size(stname),np.size(mtime)),'f')*np.nan -btm=np.zeros((np.size(stname),np.size(mtime)),'f')*np.nan -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 -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') -# loop buoys -for b in range(0,np.size(stname)): - - ahs=[] - try: - - ahs=[];atm=[];atp=[];adm=[];atime=[] - for y in yrange: - - f=nc.Dataset(ndbcp+"/"+stname[b]+"h"+repr(y)+".nc") - if 'wave_height' in f.variables.keys(): - 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(): - ahs = np.append(ahs,f.variables['swh'][:,0,0]) - - 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)) - - 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)) - - if 'mean_wave_dir' in f.variables.keys(): - adm = np.append(adm,f.variables['mean_wave_dir'][:,0,0]) - else: - adm = np.array(np.copy(ahs*nan)) - - if 'latitude' in f.variables.keys(): - lat[b] = f.variables['latitude'][:] - elif 'LATITUDE' in f.variables.keys(): - lat[b] = f.variables['LATITUDE'][:] - else: - lat[b] = nan - - if 'longitude' in f.variables.keys(): - lon[b] = f.variables['longitude'][:] - elif 'LONGITUDE' in f.variables.keys(): - lon[b] = f.variables['LONGITUDE'][:] - else: - lon[b] = nan - - atime = np.append(atime,np.array(f.variables['time'][:]).astype('double')) - - f.close(); del f - - adp = adm*np.nan # no peak direction available in this format - - except: - try: - f=nc.Dataset(copernp+"/GL_TS_MO_"+stname[b]+".nc") - if 'VHM0' in f.variables.keys(): - 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 '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(): - ahs = np.nanmean(f.variables['sea_surface_significant_wave_height'][:,:],axis=1) - elif 'SWHT' in f.variables.keys(): - ahs = np.nanmean(f.variables['SWHT'][:,:],axis=1) - elif 'wave_height_h1d3' in f.variables.keys(): - ahs = np.nanmean(f.variables['wave_height_h1d3'][:,:],axis=1) - elif 'spectral_significant_wave_height' in f.variables.keys(): - ahs = np.nanmean(f.variables['spectral_significant_wave_height'][:,:],axis=1) - - if 'VTM02' in f.variables.keys(): - atm = np.nanmean(f.variables['VTM02'][:,:],axis=1) - elif 'VGTA' in f.variables.keys(): - atm = np.nanmean(f.variables['VGTA'][:,:],axis=1) - else: - atm = ahs*nan - - if 'VTPK' in f.variables.keys(): - atp = np.nanmean(f.variables['VTPK'][:,:],axis=1) - elif 'dominant_wave_period' in f.variables.keys(): - atp = np.nanmean(f.variables['dominant_wave_period'][:,:],axis=1) - elif 'sea_surface_wave_period_at_spectral_density_maximum' in f.variables.keys(): - atp = np.nanmean(f.variables['sea_surface_wave_period_at_spectral_density_maximum'][:,:],axis=1) - else: - atp = ahs*nan - - if 'VMDR' in f.variables.keys(): - adm = np.nanmean(f.variables['VMDR'][:,:],axis=1) - else: - adm = ahs*nan - - if 'LATITUDE' in f.variables.keys(): - lat[b] = np.nanmean(f.variables['LATITUDE'][:]) - elif 'latitude' in f.variables.keys(): - lat[b] = np.nanmean(f.variables['latitude'][:]) - elif 'LAT' in f.variables.keys(): - lat[b] = np.nanmean(f.variables['LAT'][:]) - elif 'lat' in f.variables.keys(): - lat[b] = np.nanmean(f.variables['lat'][:]) - else: - lat[b] = nan - - if 'LONGITUDE' in f.variables.keys(): - lon[b] = np.nanmean(f.variables['LONGITUDE'][:]) - elif 'LON' in f.variables.keys(): - lon[b] = np.nanmean(f.variables['LON'][:]) - elif 'LONG' in f.variables.keys(): - lon[b] = np.nanmean(f.variables['LONG'][:]) - elif 'longitude' in f.variables.keys(): - lon[b] = np.nanmean(f.variables['longitude'][:]) - elif 'lon' in f.variables.keys(): - lon[b] = np.nanmean(f.variables['lon'][:]) - elif 'long' in f.variables.keys(): - lon[b] = np.nanmean(f.variables['long'][:]) - else: - lon[b] = nan - - adp = ahs*nan # no peak direction available in this format - - 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(): - 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 - indq=np.where((ahs>30.)|(ahs<0.0)) - if np.size(indq)>0: - ahs[indq]=np.nan; del indq - - indq=np.where((atm>40.)|(atm<0.0)) - if np.size(indq)>0: - atm[indq]=np.nan; del indq - - indq=np.where((atp>40.)|(atp<0.0)) - if np.size(indq)>0: - atp[indq]=np.nan; del indq - - indq=np.where((adm>360.)|(adm<-180.)) - if np.size(indq)>0: - adm[indq]=np.nan; del indq - - indq=np.where((adp>360.)|(adp<-180.)) - if np.size(indq)>0: - adp[indq]=np.nan; del indq - - c=0 - for t in range(0,np.size(mtime)): - indt=np.where(np.abs(atime-mtime[t])<1800.) - if np.size(indt)>0: - if np.any(ahs[indt[0]].mask==False): - bhs[b,t] = np.nanmean(ahs[indt[0]][ahs[indt[0]].mask==False]) - c=c+1 - if np.any(atm[indt[0]].mask==False): - btm[b,t] = np.nanmean(atm[indt[0]][atm[indt[0]].mask==False]) - if np.any(atp[indt[0]].mask==False): - btp[b,t] = np.nanmean(atp[indt[0]][atp[indt[0]].mask==False]) - if np.any(adm[indt[0]].mask==False): - 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]) - - del indt - - # print("counted "+repr(c)+" at "+stname[b]) - - print(" station "+stname[b]+" ok") - del ahs - -print(' ') # Simple quality-control (range) -ind=np.where((bhs>30.)|(bhs<0.0)) -if np.size(ind)>0: - bhs[ind]=np.nan; del ind - -ind=np.where((btm>40.)|(btm<0.0)) -if np.size(ind)>0: - btm[ind]=np.nan; del ind - -ind=np.where((btp>40.)|(btp<0.0)) -if np.size(ind)>0: - btp[ind]=np.nan; del ind - -ind=np.where((bdm>360.)|(bdm<-180.)) -if np.size(ind)>0: - bdm[ind]=np.nan; del ind - -ind=np.where((bdp>360.)|(bdp<-180.)) -if np.size(ind)>0: - bdp[ind]=np.nan; del ind - -ind=np.where((mhs>30.)|(mhs<0.0)) -if np.size(ind)>0: - mhs[ind]=np.nan; del ind - -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 - -ind=np.where((mdm>360.)|(mdm<-180.)) -if np.size(ind)>0: - mdm[ind]=np.nan; del ind - -ind=np.where((mdp>360.)|(mdp<-180.)) -if np.size(ind)>0: - mdp[ind]=np.nan; del ind +ind = np.where((bhs > 30.) | (bhs < 0.0)) +if np.size(ind) > 0: + bhs[ind] = np.nan + del ind + +ind = np.where((btm > 40.) | (btm < 0.0)) +if np.size(ind) > 0: + btm[ind] = np.nan + del ind + +ind = np.where((btp > 40.) | (btp < 0.0)) +if np.size(ind) > 0: + btp[ind] = np.nan + del ind + +ind = np.where((bdm > 360.) | (bdm < -180.)) +if np.size(ind) > 0: + bdm[ind] = np.nan + del ind + +ind = np.where((bdp > 360.) | (bdp < -180.)) +if np.size(ind) > 0: + bdp[ind] = np.nan + del ind + +ind = np.where((bwind > 50.0) | (bwind < 0.0)) +if np.size(ind) > 0: + bwind[ind] = np.nan + del ind + +ind = np.where((mhs > 30.) | (mhs < 0.0)) +if np.size(ind) > 0: + mhs[ind] = np.nan + del ind + +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 + +ind = np.where((mdm > 360.) | (mdm < -180.)) +if np.size(ind) > 0: + mdm[ind] = np.nan + del ind + +ind = np.where((mdp > 360.) | (mdp < -180.)) +if np.size(ind) > 0: + mdp[ind] = np.nan + del ind + +ind = np.where((mwn > 50.) | (mwn < 0.0)) +if np.size(ind) > 0: + mwn[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: - stname=np.array(stname[ind[0]]) - lat=np.array(lat[ind[0]]) - lon=np.array(lon[ind[0]]) - mhs=np.array(mhs[ind[0],:]) - mtm=np.array(mtm[ind[0],:]) - mtp=np.array(mtp[ind[0],:]) - mdm=np.array(mdm[ind[0],:]) - mdp=np.array(mdp[ind[0],:]) - bhs=np.array(bhs[ind[0],:]) - btm=np.array(btm[ind[0],:]) - btp=np.array(btp[ind[0],:]) - bdm=np.array(bdm[ind[0],:]) - bdp=np.array(bdp[ind[0],:]) +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: + stname = np.array(stname[ind[0]]) + lat = np.array(lat[ind[0]]) + lon = np.array(lon[ind[0]]) + mhs = np.array(mhs[ind[0], :]) + mtm = np.array(mtm[ind[0], :]) + mtp = np.array(mtp[ind[0], :]) + mdm = np.array(mdm[ind[0], :]) + mdp = np.array(mdp[ind[0], :]) + mwn = np.array(mwn[ind[0], :]) + bhs = np.array(bhs[ind[0], :]) + btm = np.array(btm[ind[0], :]) + btp = np.array(btp[ind[0], :]) + bdm = np.array(bdm[ind[0], :]) + bdp = np.array(bdp[ind[0], :]) + bwind = np.array(bwind[ind[0], :]) else: - sys.exit(' Error: No matchups Model/Buoy available.') + sys.exit(' Error: No matchups Model/Buoy available.') -print(" Matchups model/buoy complete. Total of "+repr(np.size(ind))+" stations/buoys avaliable."); del ind +print(" Matchups model/buoy complete. Total of " + repr(np.size(ind)) + " stations/buoys available.") +del ind # Processing grid and/or cyclone information -if gridinfo!=0: - print(" Adding extra information ... ") - alon=np.copy(lon); alon[alon<0]=alon[alon<0]+360. - indgplat=[]; indgplon=[] - for i in range(0,lat.shape[0]): - # indexes nearest point. - indgplat = np.append(indgplat,np.where( abs(mlat-lat[i])==abs(mlat-lat[i]).min() )[0][0]) - indgplon = np.append(indgplon,np.where( abs(mlon-alon[i])==abs(mlon-alon[i]).min() )[0][0]) - - indgplat=np.array(indgplat).astype('int'); indgplon=np.array(indgplon).astype('int') - pdistcoast=np.zeros(lat.shape[0],'f')*np.nan - pdepth=np.zeros(lat.shape[0],'f')*np.nan - poni=np.zeros(lat.shape[0],'f')*np.nan - phsmz=np.zeros(lat.shape[0],'f')*np.nan - for i in range(0,lat.shape[0]): - pdistcoast[i]=distcoast[indgplat[i],indgplon[i]] - pdepth[i]=depth[indgplat[i],indgplon[i]] - poni[i]=oni[indgplat[i],indgplon[i]] - phsmz[i]=hsmz[indgplat[i],indgplon[i]] - - print(" Grid Information Included.") - - # Excluding shallow water points too close to the coast (mask information not accurate) - ind=np.where( (np.isnan(pdistcoast)==False) & (np.isnan(pdepth)==False) ) - if np.size(ind)>0: - stname=np.array(stname[ind[0]]) - lat=np.array(lat[ind[0]]) - lon=np.array(lon[ind[0]]) - mhs=np.array(mhs[ind[0],:]) - mtm=np.array(mtm[ind[0],:]) - mtp=np.array(mtp[ind[0],:]) - mdm=np.array(mdm[ind[0],:]) - mdp=np.array(mdp[ind[0],:]) - bhs=np.array(bhs[ind[0],:]) - btm=np.array(btm[ind[0],:]) - btp=np.array(btp[ind[0],:]) - bdm=np.array(bdm[ind[0],:]) - bdp=np.array(bdp[ind[0],:]) - pdistcoast=np.array(pdistcoast[ind[0]]) - pdepth=np.array(pdepth[ind[0]]) - poni=np.array(poni[ind[0]]) - phsmz=np.array(phsmz[ind[0]]) - else: - sys.exit(' Error: No matchups Model/Buoy available after using grid mask.') - - del ind - - if cyclonemap!=0: - fcmap=np.zeros((lat.shape[0],mtime.shape[0]),'f')*np.nan - for t in range(0,np.size(mtime)): - # search for cyclone time index and cyclone map - indt=np.where(np.abs(ctime-mtime[t])<5400.) - if np.size(indt)>0: - for i in range(0,lat.shape[0]): - fcmap[i,t] = np.array(cmap[indt[0][0],indgplat[i],indgplon[i]]) - - del indt - else: - print(' - No cyclone information for this time step: '+repr(t)) - - # print(' Done cyclone analysis at step: '+repr(t)) - - ind=np.where(fcmap<0) - if np.size(ind)>0: - fcmap[ind]=np.nan - - print(" Cyclone Information Included.") +if gridinfo: + print(" Adding extra information ... ") + alon = np.copy(lon) + alon[alon < 0] = alon[alon < 0] + 360. + indgplat = [] + indgplon = [] + for i in range(0, lat.shape[0]): + # indexes nearest point. + indgplat = np.append(indgplat, np.where(abs(mlat - lat[i]) == abs(mlat - lat[i]).min())[0][0]) + indgplon = np.append(indgplon, np.where(abs(mlon - alon[i]) == abs(mlon - alon[i]).min())[0][0]) + + indgplat = np.array(indgplat).astype('int') + indgplon = np.array(indgplon).astype('int') + pdistcoast = np.zeros(lat.shape[0], 'f') * np.nan + pdepth = np.zeros(lat.shape[0], 'f') * np.nan + poni = np.zeros(lat.shape[0], 'f') * np.nan + phsmz = np.zeros(lat.shape[0], 'f') * np.nan + for i in range(0, lat.shape[0]): + pdistcoast[i] = distcoast[indgplat[i], indgplon[i]] + pdepth[i] = depth[indgplat[i], indgplon[i]] + poni[i] = oni[indgplat[i], indgplon[i]] + phsmz[i] = hsmz[indgplat[i], indgplon[i]] + + print(" Grid Information Included.") + + # Excluding shallow water points too close to the coast (mask information not accurate) + ind = np.where((np.isnan(pdistcoast) == False) & (np.isnan(pdepth) == False)) + if np.size(ind) > 0: + stname = np.array(stname[ind[0]]) + lat = np.array(lat[ind[0]]) + lon = np.array(lon[ind[0]]) + mhs = np.array(mhs[ind[0], :]) + mtm = np.array(mtm[ind[0], :]) + mtp = np.array(mtp[ind[0], :]) + mdm = np.array(mdm[ind[0], :]) + mdp = np.array(mdp[ind[0], :]) + bhs = np.array(bhs[ind[0], :]) + btm = np.array(btm[ind[0], :]) + btp = np.array(btp[ind[0], :]) + bdm = np.array(bdm[ind[0], :]) + bdp = np.array(bdp[ind[0], :]) + pdistcoast = np.array(pdistcoast[ind[0]]) + pdepth = np.array(pdepth[ind[0]]) + poni = np.array(poni[ind[0]]) + phsmz = np.array(phsmz[ind[0]]) + else: + sys.exit(' Error: No matchups Model/Buoy available after using grid mask.') + + del ind + + if cyclonemap: + fcmap = np.zeros((lat.shape[0], mtime.shape[0]), 'f') * np.nan + for t in range(0, np.size(mtime)): + # search for cyclone time index and cyclone map + indt = np.where(np.abs(ctime - mtime[t]) < 5400.) + if np.size(indt) > 0: + for i in range(0, lat.shape[0]): + fcmap[i, t] = np.array(cmap[indt[0][0], indgplat[i], indgplon[i]]) + + del indt + else: + print(' - No cyclone information for this time step: ' + repr(t)) + + ind = np.where(fcmap < 0) + if np.size(ind) > 0: + fcmap[ind] = np.nan + + print(" Cyclone Information Included.") # Edit format if this is forecast model data. Reshape and allocate -if forecastds>0: - unt=np.unique(mfcycle); mxsz=1 - for i in range(0,unt.shape[0]): - ind=np.where(mfcycle==unt[i])[0] - mxsz=np.max([mxsz,np.size(ind)]) - - for i in range(0,unt.shape[0]): - ind=np.where(mfcycle==unt[i])[0] - if i==0: - nmhs=np.zeros((mhs.shape[0],unt.shape[0],mxsz),'f')*np.nan - nmtm=np.zeros((mhs.shape[0],unt.shape[0],mxsz),'f')*np.nan - 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 - 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 - 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 - 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]) - 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]) - 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]) - nbdm[:,i,:][:,0:np.size(ind)]=np.array(bdm[:,ind]) - nbdp[:,i,:][:,0:np.size(ind)]=np.array(bdp[:,ind]) - if cyclonemap!=0: - nfcmap[:,i,:][:,0:np.size(ind)]=np.array(fcmap[:,ind]) - - ind=np.where( (nmhs>0.0) & (nbhs>0.0) ) +if forecastds > 0: + unt = np.unique(mfcycle) + mxsz = 1 + + for i in range(0, unt.shape[0]): + ind = np.where(mfcycle == unt[i])[0] + mxsz = np.max([mxsz, np.size(ind)]) + + for i in range(0, unt.shape[0]): + ind = np.where(mfcycle == unt[i])[0] + if i == 0: + nmhs = np.zeros((mhs.shape[0], unt.shape[0], mxsz), 'f') * np.nan + nmtm = np.zeros((mhs.shape[0], unt.shape[0], mxsz), 'f') * np.nan + 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 + nmwn = 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 + 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 + nbwind = np.zeros((mhs.shape[0], unt.shape[0], mxsz), 'f') * np.nan + nmtime = np.zeros((unt.shape[0], mxsz), 'double') * np.nan + if cyclonemap: + 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]) + 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]) + nmwn[:, i, :][:, 0:np.size(ind)] = np.array(mwn[:, 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]) + nbdm[:, i, :][:, 0:np.size(ind)] = np.array(bdm[:, ind]) + nbdp[:, i, :][:, 0:np.size(ind)] = np.array(bdp[:, ind]) + nbwind[:, i, :][:, 0:np.size(ind)] = np.array(bwind[:, ind]) + if cyclonemap: + nfcmap[:, i, :][:, 0:np.size(ind)] = np.array(fcmap[:, ind]) + + ind = np.where((nmhs > 0.0) & (nbhs > 0.0)) else: - ind=np.where( (mhs>0.0) & (bhs>0.0) ) - - -if np.size(ind)>0: - print(' Total amount of matchups model/buoy: '+repr(np.size(ind))) - - # 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) - ncfile = nc.Dataset('WW3.Buoy'+ftag+'_'+initime+'to'+fintime+'.nc', "w", format=fnetcdf) - ncfile.history="Matchups of WAVEWATCHIII point output (table) and NDBC and Copernicus Buoys. Total of "+repr(bhs[bhs>0.].shape[0])+" observations or pairs model/observation." - # create dimensions - ncfile.createDimension('buoypoints', bhs.shape[0] ) - if gridinfo!=0: - ncfile.createDimension('GlobalOceansSeas', ocnames.shape[0] ) - ncfile.createDimension('HighSeasMarineZones', hsmznames.shape[0] ) - if cyclonemap!=0: - ncfile.createDimension('cycloneinfo', cinfo.shape[0] ) - vcinfo = ncfile.createVariable('cycloneinfo',dtype('a25'),('cycloneinfo')) - # create variables. - vstname = ncfile.createVariable('buoyID',dtype('a25'),('buoypoints')) - vlat = ncfile.createVariable('latitude',np.dtype('float32').char,('buoypoints')) - vlon = ncfile.createVariable('longitude',np.dtype('float32').char,('buoypoints')) - - if forecastds>0: - ncfile.createDimension('time', nmhs.shape[2] ) - ncfile.createDimension('fcycle', unt.shape[0] ) - vt = ncfile.createVariable('time',np.dtype('float64').char,('fcycle','time')) - vmhs = ncfile.createVariable('model_hs',np.dtype('float32').char,('buoypoints','fcycle','time')) - vmtm = ncfile.createVariable('model_tm',np.dtype('float32').char,('buoypoints','fcycle','time')) - 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')) - 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')) - vbdm = ncfile.createVariable('obs_dm',np.dtype('float32').char,('buoypoints','fcycle','time')) - vbdp = ncfile.createVariable('obs_dp',np.dtype('float32').char,('buoypoints','fcycle','time')) - else: - ncfile.createDimension('time', bhs.shape[1] ) - vt = ncfile.createVariable('time',np.dtype('float64').char,('time')) - vmhs = ncfile.createVariable('model_hs',np.dtype('float32').char,('buoypoints','time')) - vmtm = ncfile.createVariable('model_tm',np.dtype('float32').char,('buoypoints','time')) - 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')) - 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')) - vbdm = ncfile.createVariable('obs_dm',np.dtype('float32').char,('buoypoints','time')) - 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')) - 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')) - vhsmznames = ncfile.createVariable('names_HighSeasMarineZones',dtype('a25'),('HighSeasMarineZones')) - if cyclonemap!=0: - if forecastds>0: - vcmap = ncfile.createVariable('cyclone',np.dtype('float32').char,('buoypoints','fcycle','time')) - else: - vcmap = ncfile.createVariable('cyclone',np.dtype('float32').char,('buoypoints','time')) - - # Assign units - vlat.units = 'degrees_north' ; vlon.units = 'degrees_east' - vt.units = 'seconds since 1970-01-01T00:00:00+00:00' - vmhs.units='m'; vbhs.units='m' - vmtm.units='s'; vbtm.units='s' - vmtp.units='s'; vbtp.units='s' - vmdm.units='degrees'; vbdm.units='degrees' - vmdp.units='degrees'; vbdp.units='degrees' - if gridinfo!=0: - vpdepth.units='m'; vpdistcoast.units='km' - - # Allocate Data - vstname[:]=stname[:]; vlat[:] = lat[:]; vlon[:] = lon[:] - if forecastds>0: - vt[:,:]=nmtime[:,:] - vmhs[:,:,:]=nmhs[:,:,:] - vmtm[:,:,:]=nmtm[:,:,:] - vmtp[:,:,:]=nmtp[:,:,:] - vmdm[:,:,:]=nmdm[:,:,:] - vmdp[:,:,:]=nmdp[:,:,:] - vbhs[:,:,:]=nbhs[:,:,:] - vbtm[:,:,:]=nbtm[:,:,:] - vbtp[:,:,:]=nbtp[:,:,:] - vbdm[:,:,:]=nbdm[:,:,:] - vbdp[:,:,:]=nbdp[:,:,:] - else: - vt[:]=mtime[:] - vmhs[:,:]=mhs[:,:] - vmtm[:,:]=mtm[:,:] - vmtp[:,:]=mtp[:,:] - vmdm[:,:]=mdm[:,:] - vmdp[:,:]=mdp[:,:] - vbhs[:,:]=bhs[:,:] - vbtm[:,:]=btm[:,:] - vbtp[:,:]=btp[:,:] - vbdm[:,:]=bdm[:,:] - vbdp[:,:]=bdp[:,:] - - if gridinfo!=0: - vpdistcoast[:]=pdistcoast[:] - vpdepth[:]=pdepth[:] - vponi[:]=poni[:]; vocnames[:] = ocnames[:] - vphsmz[:]=phsmz[:]; vhsmznames[:] = hsmznames[:] - if cyclonemap!=0: - vcinfo[:] = cinfo[:] - if forecastds>0: - vcmap[:,:,:]=nfcmap[:,:,:] - else: - vcmap[:,:]=fcmap[:,:] - - ncfile.close() - print(' ') - print('Done. Netcdf ok. New file saved: WW3.Buoy'+ftag+'_'+initime+'to'+fintime+'.nc') + ind = np.where((mhs > 0.0) & (bhs > 0.0)) + +if np.size(ind) > 0: + print(' Total amount of matchups model/buoy: ' + repr(np.size(ind))) + + # 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) + + # Ensure model_name is defined + if 'model_name' not in locals(): + model_name = '' + + ncfile = nc.Dataset(f'WW3.{model_name}Buoy{ftag}_{initime}to{fintime}.nc', "w", format=fnetcdf) + print(f"Model Name: {model_name}, Tag: {ftag}, Start Time: {initime}, End Time: {fintime}") + ncfile.history = "Matchups of WAVEWATCHIII point output (table) and NDBC and Copernicus Buoys. Total of " + repr(bhs[bhs > 0.].shape[0]) + " observations or pairs model/observation." + ncfile.initial_condition = initime + ncfile.time_units = "seconds since 1970-01-01T00:00:00+00:00" + + # Create dimensions + ncfile.createDimension('buoypoints', bhs.shape[0]) + if gridinfo: + ncfile.createDimension('GlobalOceansSeas', ocnames.shape[0]) + ncfile.createDimension('HighSeasMarineZones', hsmznames.shape[0]) + if cyclonemap: + ncfile.createDimension('cycloneinfo', cinfo.shape[0]) + vcinfo = ncfile.createVariable('cycloneinfo', 'S1', ('cycloneinfo')) + + # Create variables + vstname = ncfile.createVariable('buoyID', 'S1', ('buoypoints')) + vlat = ncfile.createVariable('latitude', 'f4', ('buoypoints')) + vlon = ncfile.createVariable('longitude', 'f4', ('buoypoints')) + + if forecastds > 0: + ncfile.createDimension('time', nmhs.shape[2]) + ncfile.createDimension('fcycle', unt.shape[0]) + vt = ncfile.createVariable('time', 'f8', ('fcycle', 'time')) + vmhs = ncfile.createVariable('model_hs', 'f4', ('buoypoints', 'fcycle', 'time')) + vmtm = ncfile.createVariable('model_tm', 'f4', ('buoypoints', 'fcycle', 'time')) + vmtp = ncfile.createVariable('model_tp', 'f4', ('buoypoints', 'fcycle', 'time')) + vmdm = ncfile.createVariable('model_dm', 'f4', ('buoypoints', 'fcycle', 'time')) + vmdp = ncfile.createVariable('model_dp', 'f4', ('buoypoints', 'fcycle', 'time')) + vmwn = ncfile.createVariable('model_wind', 'f4', ('buoypoints', 'fcycle', 'time')) + + vbhs = ncfile.createVariable('obs_hs', 'f4', ('buoypoints', 'fcycle', 'time')) + vbtm = ncfile.createVariable('obs_tm', 'f4', ('buoypoints', 'fcycle', 'time')) + vbtp = ncfile.createVariable('obs_tp', 'f4', ('buoypoints', 'fcycle', 'time')) + vbdm = ncfile.createVariable('obs_dm', 'f4', ('buoypoints', 'fcycle', 'time')) + vbdp = ncfile.createVariable('obs_dp', 'f4', ('buoypoints', 'fcycle', 'time')) + vbwind = ncfile.createVariable('obs_wind', 'f4', ('buoypoints', 'fcycle', 'time')) + + else: + ncfile.createDimension('time', bhs.shape[1]) + vt = ncfile.createVariable('time', 'f8', ('time')) + vmhs = ncfile.createVariable('model_hs', 'f4', ('buoypoints', 'time')) + vmtm = ncfile.createVariable('model_tm', 'f4', ('buoypoints', 'time')) + vmtp = ncfile.createVariable('model_tp', 'f4', ('buoypoints', 'time')) + vmdm = ncfile.createVariable('model_dm', 'f4', ('buoypoints', 'time')) + vmdp = ncfile.createVariable('model_dp', 'f4', ('buoypoints', 'time')) + vmwn = ncfile.createVariable('model_wind', 'f4', ('buoypoints', 'time')) + + vbhs = ncfile.createVariable('obs_hs', 'f4', ('buoypoints', 'time')) + vbtm = ncfile.createVariable('obs_tm', 'f4', ('buoypoints', 'time')) + vbtp = ncfile.createVariable('obs_tp', 'f4', ('buoypoints', 'time')) + vbdm = ncfile.createVariable('obs_dm', 'f4', ('buoypoints', 'time')) + vbdp = ncfile.createVariable('obs_dp', 'f4', ('buoypoints', 'time')) + vbwind = ncfile.createVariable('obs_wind', 'f4', ('buoypoints', 'time')) + + if gridinfo: + vpdistcoast = ncfile.createVariable('distcoast', 'f4', ('buoypoints')) + vpdepth = ncfile.createVariable('depth', 'f4', ('buoypoints')) + vponi = ncfile.createVariable('GlobalOceansSeas', 'f4', ('buoypoints')) + vocnames = ncfile.createVariable('names_GlobalOceansSeas', 'S1', ('GlobalOceansSeas')) + vphsmz = ncfile.createVariable('HighSeasMarineZones', 'f4', ('buoypoints')) + vhsmznames = ncfile.createVariable('names_HighSeasMarineZones', 'S1', ('HighSeasMarineZones')) + if cyclonemap: + if forecastds > 0: + vcmap = ncfile.createVariable('cyclone', 'f4', ('buoypoints', 'fcycle', 'time')) + else: + vcmap = ncfile.createVariable('cyclone', 'f4', ('buoypoints', 'time')) + + # Assign units + vlat.units = 'degrees_north' + vlon.units = 'degrees_east' + vt.units = 'seconds since 1970-01-01T00:00:00+00:00' + vmhs.units = 'm' + vbhs.units = 'm' + vmtm.units = 's' + vbtm.units = 's' + vmtp.units = 's' + vbtp.units = 's' + vmdm.units = 'degrees' + vbdm.units = 'degrees' + vmdp.units = 'degrees' + vbdp.units = 'degrees' + vmwn.unit = 'm/s' + vbwind.unit = 'm/s' + + if gridinfo: + vpdepth.units = 'm' + vpdistcoast.units = 'km' + + # Allocate Data + vstname[:] = stname[:] + vlat[:] = lat[:] + vlon[:] = lon[:] + initime_unix = timegm(strptime(initime, '%Y%m%d%H')) + + if forecastds > 0: + vt[:, :] = nmtime[:, :] + vmhs[:, :, :] = nmhs[:, :, :] + vmtm[:, :, :] = nmtm[:, :, :] + vmtp[:, :, :] = nmtp[:, :, :] + vmdm[:, :, :] = nmdm[:, :, :] + vmdp[:, :, :] = nmdp[:, :, :] + vmwn[:, :, :] = nmwn[:, :, :] + vbhs[:, :, :] = nbhs[:, :, :] + vbtm[:, :, :] = nbtm[:, :, :] + vbtp[:, :, :] = nbtp[:, :, :] + vbdm[:, :, :] = nbdm[:, :, :] + vbdp[:, :, :] = nbdp[:, :, :] + vbwind[:, :, :] = nbwind[:, :, :] + + if gridinfo: + vpdistcoast[:] = pdistcoast[:] + vpdepth[:] = pdepth[:] + vponi[:] = poni[:] + vocnames[:] = ocnames[:] + vphsmz[:] = phsmz[:] + vhsmznames[:] = hsmznames[:] + if cyclonemap: + vcinfo[:] = cinfo[:] + if forecastds > 0: + vcmap[:, :, :] = nfcmap[:, :, :] + else: + vcmap[:, :] = fcmap[:, :] + + fcst_hr = np.full_like(mtime, np.nan, dtype='float64') + for i in range(mtime.shape[0]): + if mtime.ndim > 1: + for j in range(mtime.shape[1]): + if not np.isnan(mtime[i, j]): + fcst_hr[i, j] = (mtime[i, j] - initime_unix) / 3600 + else: + if not np.isnan(mtime[i]): + fcst_hr[i] = (mtime[i] - initime_unix) / 3600 + + # Add fcst_hr + vfcst_hr = ncfile.createVariable('fcst_hr', 'f8', ('buoypoints', 'time')) + vfcst_hr.units = 'hours' + vfcst_hr[:] = fcst_hr + + ncfile.close() + print(' ') + print(f'Done. Netcdf ok. New file saved: WW3.{model_name}_Buoy{ftag}_{initime}to{fintime}.nc') + + # File deletion section + files_to_delete = [ + os.path.join(output_dir, '*_contents.txt'), + os.path.join(output_dir, '*_id.txt'), + os.path.join(output_dir, '*.spec_tar') + ] + + for pattern in files_to_delete: + for file in glob.glob(pattern): + try: + os.remove(file) + print(f'Deleted file: {file}') + except OSError as e: + print(f'Error deleting file {file}: {e}') + + # Folder deletion section + try: + shutil.rmtree(extracted_folder) + print(f'Deleted folder: {extracted_folder}') + except OSError as e: + print(f'Error deleting folder {extracted_folder}: {e}') + + print('Temporary files and folder deleted.') diff --git a/ww3tools/wread.py b/ww3tools/wread.py old mode 100755 new mode 100644 index 74ce53c..08c296a --- a/ww3tools/wread.py +++ b/ww3tools/wread.py @@ -60,7 +60,6 @@ import matplotlib import time -import timeit from time import strptime from calendar import timegm import pandas as pd @@ -68,14 +67,15 @@ import netCDF4 as nc import numpy as np from pylab import * -import yaml -import re -import os +import matplotlib.pyplot as plt import sys +import pandas as pd from matplotlib import ticker # import pickle import sys import warnings; warnings.filterwarnings("ignore") +import tarfile +import math def readconfig(fname): @@ -1619,3 +1619,184 @@ def spec_ww3(*args): del mtime,mdate,lat,lon,wnds,wndd,freq,freq1,freq2,dfreq,pwst,dire,d1sp,dspec + + +def spec1_ww3(*args): + ''' + WAVEWATCH III, wave spectrum, netcdf (.nc) or text (.spec) format + Input: file names (list of file names), and station names (list of station names) + Output: list of dictionaries containing: + time(seconds since 1970), time(datetime64), lat, lon; Arrays: freq, dfreq, pwst, d1sp, dire, dspec, wnds, wndd + ''' + + if len(args) < 2: + sys.exit('Two inputs are required: list of file names and list of station names') + + fnames = args[0] + stnames = args[1] + sk = 1 + if len(args) > 2: + sk = int(args[2]) + if len(args) > 3: + sys.exit('Too many inputs') + + results = [] + + for fname in fnames: + for stname in stnames: + try: + with open(fname) as fp: + nt = fp.read().count(stname) + + if nt >= 1: + with open(fname) as fp: + cabc = fp.readline().strip().split() + nf = int(cabc[3]) + nd = int(cabc[4]) + npo = int(cabc[5]) + + freq = np.zeros(nf, 'f') + dire = np.zeros(nd, 'f') + dspec = np.zeros((nt, nf, nd), 'f') + adire = np.zeros(dire.shape) + adspec = np.zeros(dspec.shape) + mtime = np.zeros((nt), 'd') + + k = 0 + # Reading frequencies + for i in range(0, int(np.floor(nf / 8))): + line = fp.readline().strip().split() + for j in range(8): + freq[k] = float(line[j]) + k += 1 + + if (nf % 8) > 0: + line = fp.readline().strip().split() + for i in range(nf % 8): + freq[k] = float(line[i]) + k += 1 + + # Calculate dfreq using geometric progression + dfreq = np.zeros(freq.shape[0], 'f') + alpha = (freq[-1] / freq[-2]) + for i in range(freq.shape[0]): + if i == 0: + dfreq[i] = freq[i] * (np.sqrt(alpha) - 1) + elif i == (freq.shape[0] - 1): + dfreq[i] = freq[i] * (1 - 1 / np.sqrt(alpha)) + else: + dfreq[i] = freq[i] * (np.sqrt(alpha) - 1 / np.sqrt(alpha)) + + k = 0 + # Reading directions + for i in range(0, int(np.floor(nd / 7))): + line = fp.readline().strip().split() + for j in range(7): + dire[k] = float(line[j]) * 180 / np.pi + k += 1 + + if (nd % 7) > 0: + line = fp.readline().strip().split() + for i in range(nd % 7): + dire[k] = float(line[i]) * 180 / np.pi + k += 1 + + auxs = np.zeros((nf * nd), 'f') + wnds = np.zeros((nt), 'f') + wndd = np.zeros((nt), 'f') + + for t in range(nt): + cabc = fp.readline().strip().split() + mtime[t] = np.double(timegm(strptime(cabc[0] + cabc[1][0:2], '%Y%m%d%H'))) + cabc = fp.readline().strip().split() + + if len(cabc) > 8: + lat, lon, depth = float(cabc[2]), float(cabc[3]), float(cabc[4]) + wnds[t], wndd[t] = float(cabc[5]), float(cabc[6]) + elif len(cabc) == 8: + lat_lon_parts = cabc[2].strip("'").split('-') + if len(lat_lon_parts) == 2: + lat, lon, depth = float(lat_lon_parts[0]), -float(lat_lon_parts[1]), float(cabc[3]) + else: + lat = float(cabc[2][:6]) + lon = float(cabc[2][6:]) + depth = float(cabc[3]) + wnds[t], wndd[t] = float(cabc[4]), float(cabc[5]) + elif len(cabc) == 7: + + lat_lon_parts = cabc[1].split() + lat, lon = float(lat_lon_parts[0]), float(lat_lon_parts[1]) + depth = float(cabc[2]) + wnds[t], wndd[t] = float(cabc[3]), float(cabc[4]) + else: + continue + + k = 0 + for i in range(0, int(np.floor((nf * nd) / 7.))): + line = fp.readline().strip().split() + for j in range(7): + auxs[k] = float(line[j]) + k += 1 + + if (nf * nd % 7) > 0: + line = fp.readline().strip().split() + for i in range(nf * nd % 7): + auxs[k] = float(line[i]) + k += 1 + + for ic in range(nf): + for il in range(nd): + dspec[t, ic, il] = auxs[il * nf + ic] + + mdate = pd.to_datetime(mtime, unit='s').strftime('%Y-%m-%dT%H:%M:%S.%f') + freq1 = np.copy(freq) + freq2 = np.copy(freq) + + pwst = np.zeros((dspec.shape[0], nf), 'f') + for t in range(dspec.shape[0]): + for il in range(nf): + pwst[t, il] = sum(dspec[t, il, :] * (2 * np.pi) / nd) + pwst[t, :] *= dfreq + + adspec = np.copy(dspec) + inddire = np.argmin(dire) + for t in range(dspec.shape[0]): + adspec[t, :, 0:nd - (inddire + 1)] = dspec[t, :, (inddire + 1):] + adspec[t, :, nd - (inddire + 1):nd] = dspec[t, :, :(inddire + 1)] + for i in range(nd): + dspec[t, :, i] = adspec[t, :, nd - i - 1] + adspec[t, :, :int(nd / 2)] = dspec[t, :, int(nd / 2):] + adspec[t, :, int(nd / 2):] = dspec[t, :, :int(nd / 2)] + dspec[t, :, :] = adspec[t, :, :] + + dire = np.sort(dire) + + d1sp = np.zeros((dspec.shape[0], nf), 'f') + for t in range(dspec.shape[0]): + for il in range(nf): + a = np.sum(dspec[t, il, :] * np.sin((np.pi * dire) / 180.) / np.sum(dspec[t, il, :])) + b = np.sum(dspec[t, il, :] * np.cos((np.pi * dire) / 180.) / np.sum(dspec[t, il, :])) + aux = np.degrees(np.arctan2(a, b)) + if aux < 0: + aux += 360. + d1sp[t, il] = aux + + m0 = np.sum(pwst, axis=1) + hs = 4 * np.sqrt(m0) + max_index = np.argmax(pwst, axis=1) + f_max = freq[max_index] + tp = 1 / f_max + + result = { + 'time': mtime, 'date': mdate, 'latitude': lat, 'longitude': lon, 'depth': depth, + 'wind_spd': wnds, 'wind_dir': wndd, 'freq': freq, 'freq1': freq1, 'freq2': freq2, + 'deltafreq': dfreq, 'pspec': pwst, 'theta': dire, 'dmspec': d1sp, 'dirspec': dspec, + 'Hs': hs, 'Tp': tp, 'station_name': stname + } + + results.append(result) + except Exception as e: + print(f"Skipping file {fname} for station {stname}: {str(e)}") + continue + + return results