From 5e7a70e47b98dd4427c9a773284419bf8ea42e92 Mon Sep 17 00:00:00 2001 From: cuill Date: Wed, 19 Oct 2022 10:08:17 -0400 Subject: [PATCH 01/31] Add function to get element centroid and area --- pyschism/mesh/base.py | 75 ++++++++++++++++++++++++++++++++++--------- 1 file changed, 60 insertions(+), 15 deletions(-) diff --git a/pyschism/mesh/base.py b/pyschism/mesh/base.py index cf0dc0d7..4c9e4221 100644 --- a/pyschism/mesh/base.py +++ b/pyschism/mesh/base.py @@ -268,6 +268,66 @@ def get_node_ball(self): ine = np.array([np.array(ine[i]) for i in np.arange(NP)], dtype='O') return nne, ine + def compute_centroid(self): + elnode = self.array + NE = self.elements.__len__() + depth = self.nodes.values + + x_centr, y_centr, dp_centr = np.zeros([3, NE]) + + if np.any(elnode.mask): + mask = elnode.mask[:, -1] + #centroid for tris + x_centr[mask] = self.nodes.coords[elnode[mask, :3], 0].mean(axis=1) + y_centr[mask] = self.nodes.coords[elnode[mask, :3], 1].mean(axis=1) + dp_centr[mask] = depth[elnode[mask, :3]].mean(axis=1) + + #centroid for quads + x_centr[~mask] = self.nodes.coords[elnode[~mask, :], 0].mean(axis=1) + y_centr[~mask] = self.nodes.coords[elnode[~mask, :], 1].mean(axis=1) + dp_centr[~mask] = depth[elnode[mask, :]].mean(axis=1) + else: + x_centr = self.nodes.coords[elnode, 0].mean(axis=1) + y_centr = self.nodes.coords[elnode, 1].mean(axis=1) + dp_centr = depth[elnode].mean(axis=1) + + return x_centr, y_centr, dp_centr + + def get_areas(self, crs: Union[CRS, str] = None): + xy = self.nodes.get_xy(crs) + x = xy[:, 0] + y = xy[:, 1] + + elnode = self.array + x1 = x[elnode[:, 0]]; y1 = y[elnode[:, 0]] + x2 = x[elnode[:, 1]]; y2 = y[elnode[:, 1]] + x3 = x[elnode[:, 2]]; y3 = y[elnode[:, 2]] + if np.any(elnode.mask): + x4 = x[elnode[:, 3]]; y4 = y[elnode[:, 3]] + mask = elnode.mask[:, -1] + x4[mask] = x1[mask]; y4[mask] = y1[mask] + else: + x4 = x1; y4 = y1 + area=((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)+(x3-x1)*(y4-y1)-(x4-x1)*(y3-y1))/2 + + return area + + ##self.gdf method is very slow + #if self.nodes.crs.is_geographic: + # elements = [] + # for row in self.gdf.itertuples(): + # aeqd = CRS.from_user_input( + # f"+proj=aeqd +R=6371000 +units=m " + # f"+lat_0={row.geometry.centroid.y} +lon_0={row.geometry.centroid.x}" + # ) + # current_to_aeqd = Transformer.from_crs( + # self.nodes.crs, aeqd, always_xy=True + # ).transform + # elements.append(ops.transform(current_to_aeqd, row.geometry)) + # return [element.area for element in elements] + #else: + # return [row.geometry.area for row in self.gdf.itertuples()] + def get_triangulation_mask(self, element_mask): triangulation_mask = [] @@ -288,21 +348,6 @@ def get_triangulation_mask(self, element_mask): return np.array(triangulation_mask) - def get_areas(self): - if self.nodes.crs.is_geographic: - elements = [] - for row in self.gdf.itertuples(): - aeqd = CRS.from_user_input( - f"+proj=aeqd +R=6371000 +units=m " - f"+lat_0={row.geometry.centroid.y} +lon_0={row.geometry.centroid.x}" - ) - current_to_aeqd = Transformer.from_crs( - self.nodes.crs, aeqd, always_xy=True - ).transform - elements.append(ops.transform(current_to_aeqd, row.geometry)) - return [element.area for element in elements] - else: - return [row.geometry.area for row in self.gdf.itertuples()] @property def array(self): From 437c10a8c931c536b6be8d366f572613f0dbbabf Mon Sep 17 00:00:00 2001 From: cuill Date: Wed, 19 Oct 2022 10:16:56 -0400 Subject: [PATCH 02/31] fixed an error --- pyschism/mesh/base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyschism/mesh/base.py b/pyschism/mesh/base.py index 4c9e4221..ff90cf00 100644 --- a/pyschism/mesh/base.py +++ b/pyschism/mesh/base.py @@ -285,7 +285,7 @@ def compute_centroid(self): #centroid for quads x_centr[~mask] = self.nodes.coords[elnode[~mask, :], 0].mean(axis=1) y_centr[~mask] = self.nodes.coords[elnode[~mask, :], 1].mean(axis=1) - dp_centr[~mask] = depth[elnode[mask, :]].mean(axis=1) + dp_centr[~mask] = depth[elnode[~mask, :]].mean(axis=1) else: x_centr = self.nodes.coords[elnode, 0].mean(axis=1) y_centr = self.nodes.coords[elnode, 1].mean(axis=1) From 18363b642f9ddd63abd31eca5f352cb5babf3719 Mon Sep 17 00:00:00 2001 From: cuill Date: Wed, 19 Oct 2022 11:51:15 -0400 Subject: [PATCH 03/31] added script to generate nc format of source/sink and convert prec to vsource --- examples/combine_sources/README | 40 +++ .../combine_sources/generate_sourcesink_nc.py | 98 +++++++ .../source_sink/sflux2source/SourceSinkIn.py | 250 ++++++++++++++++++ .../source_sink/sflux2source/TimeHistory.py | 238 +++++++++++++++++ 4 files changed, 626 insertions(+) create mode 100644 examples/combine_sources/README create mode 100644 examples/combine_sources/generate_sourcesink_nc.py create mode 100644 pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py create mode 100644 pyschism/forcing/source_sink/sflux2source/TimeHistory.py diff --git a/examples/combine_sources/README b/examples/combine_sources/README new file mode 100644 index 00000000..98af4a8f --- /dev/null +++ b/examples/combine_sources/README @@ -0,0 +1,40 @@ +generate_sourcesink_nc.py +The script will combine two sources (NWM + sflux_prc) and write vsource/msource/vsink to one netcdf file. + +Required inputs: +1. path to the run directory: + rundir = "/sciclone/schism10/lcui01/schism20/ICOGS/ICOGS2D/GA_RUN02" + +2. model start time: + model_start_time_str = ""2018-08-01 00:00:00" + +The script assumes the following files are in the rundir/: + hgrid.gr3 + + source_sink.in # These four + vsource.th + vsink.th + msource.th + + sflux/flux_prc_1.****.nc (precipitation will be interpolated to elements and converted to volume source) + +Finally, the script will write one netcdf file into rundir/Orignial+sflux_source_sink/source.nc, which looks like this: +netcdf source { +dimensions: + nsources = 1283607 ; + nsinks = 34 ; + ntracers = 2 ; + time_msource = 2 ; + time_vsource = 1489 ; + time_vsink = 1489 ; + one = 1 ; +variables: + int source_elem(nsources) ; + float vsource(time_vsource, nsources) ; + float msource(time_msource, ntracers, nsources) ; + int sink_elem(nsinks) ; + float vsink(time_vsink, nsinks) ; + float time_step_vsource(one) ; + float time_step_msource(one) ; + float time_step_vsink(one) ; +} diff --git a/examples/combine_sources/generate_sourcesink_nc.py b/examples/combine_sources/generate_sourcesink_nc.py new file mode 100644 index 00000000..0cce0dea --- /dev/null +++ b/examples/combine_sources/generate_sourcesink_nc.py @@ -0,0 +1,98 @@ +from time import time +import numpy as np +from glob import glob +from scipy import spatial +import xarray as xr +import pandas as pd + + +from pyschism.mesh import Hgrid +from pyschism.forcing.source_sink.sflux2source.SourceSinkIn import source_sink + + +def nearest_neighbour(points_a, points_b): + tree = spatial.cKDTree(points_b) + return tree.query(points_a)[1] + +def Sflux_var_at_points(sflux_dir, sflux_set='air_1', points=None, var_name=None, nt=None): + ''' + extract variable values at specified points from SCHISM-format sflux files + ''' + + sflux_files = sorted(glob(f'{sflux_dir}/sflux_{sflux_set}.*.nc')) + # read the first sflux file + sflux = xr.open_dataset(sflux_files[0]) + if nt is None: + nt = len(sflux['time'].data) + sflux.close() + + idx = nearest_neighbour(np.c_[points[:, 0], points[:, 1]], np.c_[sflux.lon.data.ravel(), sflux.lat.data.ravel()]) + sflux_files = sorted(glob(f'{sflux_dir}/sflux_{sflux_set}.*.nc')) + nt_all = nt * len(sflux_files) + npts = len(idx) + + var_at_pts = np.empty((nt_all, npts), dtype=float) * np.nan + time = np.empty((nt_all), dtype='datetime64[ns]') + + for i, sflux_file in enumerate(sflux_files): + sflux = xr.open_dataset(sflux_file) + var = sflux[var_name].data + var = var.reshape(*var.shape[:-2], -1) + var_at_pts[i*nt:(i+1)*nt, :] = var[:nt, idx] + time[i*nt:(i+1)*nt] = sflux.time.data[:nt] + print(f'processing {sflux_file}') + sflux.close() + + return [time, var_at_pts] + + +if __name__ == "__main__": + rundir = '/sciclone/schism10/lcui01/schism20/ICOGS/ICOGS2D/GA_RUN02' + model_start_time_str = "2018-08-01 00:00:00" + + # read hgrid and get element info + hgrid = Hgrid.open(f'{rundir}/hgrid.gr3', crs='epsg:4326') + lon_ctr, lat_ctr, dp_ctr = hgrid.elements.compute_centroid() + #get element area (crs=26918 to project the coordinate) + area = hgrid.elements.get_areas(26918) + + ## read the original source/sink files + original_source_sink = source_sink(source_dir=rundir, start_time_str=model_start_time_str) + + # make another set of source/sink files based on sflux + # read variable values from sflux + var_at_ele_center = Sflux_var_at_points(f'{rundir}/sflux/', sflux_set='prc_1', points=np.c_[lon_ctr, lat_ctr], var_name='prate') + # read time info from sflux + times = pd.to_datetime(var_at_ele_center[0].astype(str)) + start_time_str = times[0].strftime("%Y-%m-%d %H:%M:%S") + timedeltas = (times-times[0]).total_seconds().astype(int) + nt = times.shape[0] + # calculate new sources for all elements: prate x area + sources_sflux = np.zeros((nt, hgrid.elements.__len__()), dtype=float) + for it in range(nt): # for each time + sources_sflux[it, :] = 1e-3 * var_at_ele_center[1][it] * area # prate is in kg/m^2/s, i.e., 1e-3 m3/m2/s + + ## make new source sink + added_ss = source_sink( + source_dir=None, + source_eles=(np.array(range(hgrid.elements.__len__()))+1).tolist(), + sink_eles=[1], # dummy + start_time_str=start_time_str, + timedeltas=timedeltas, + vsource_data=sources_sflux + ) + added_ss.update_vars() + + ## add original source and sflux source + print('Adding two sources:') + t0 = time() + total_ss = original_source_sink + added_ss + if np.isnan(total_ss.vsource.data.astype(float)).any(): + raise Exception('nan found in sources') + if np.isnan(total_ss.vsink.data.astype(float)).any(): + raise Exception('nan found in sinks') + + print(f'It took {(time()-t0)/60} minutes to add two sources') + + ## write + total_ss.nc_writer(f'{rundir}/Orignial+sflux_source_sink/') diff --git a/pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py b/pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py new file mode 100644 index 00000000..05136be2 --- /dev/null +++ b/pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py @@ -0,0 +1,250 @@ +import os +import copy + +import numpy as np +from scipy import interpolate +import pandas as pd +from netCDF4 import Dataset + +from pyschism.forcing.source_sink.sflux2source.TimeHistory import TimeHistory + +def BinA(A, B): + x = A + y = B + + index = np.argsort(x) + sorted_x = x[index] + sorted_index = np.searchsorted(sorted_x, y) + + yindex = np.take(index, sorted_index, mode="clip") + mask = x[yindex] != y + + result = np.ma.array(yindex, mask=mask) + + return [np.logical_not(np.ma.getmask(result)), np.ma.getmask(result)] + + +class SourceSinkIn(): + """ class for *.prop or other similar formats""" + def __init__(self, filename=None, number_of_groups=2, ele_groups=[]): + self.n_group = number_of_groups # 0: source; 1: sink + if filename is not None: + """Read a usgs data file and initialize from its content """ + self.source_file = filename + self.np_group = [] + self.ip_group = [] + + with open(self.source_file) as fin: + for k in range(0, self.n_group): + self.np_group.append(int(fin.readline().split()[0])) + print("Points in Group " + str(k+1) + ": " + str(self.np_group[k])) + self.ip_group.append(np.empty((self.np_group[k]), dtype=int)) + for i in range(0, self.np_group[k]): + self.ip_group[k][i] = int(fin.readline()) + fin.readline() # empty line between groups + if self.np_group[k] > 0: + print("p first: " + str(self.ip_group[k][0])) + print("p last: " + str(self.ip_group[k][-1])) + else: + self.np_group = [len(x) for x in ele_groups] + self.ip_group = [np.array(x) for x in ele_groups] + + def writer(self, filename=None): + if filename is None: + filename = self.source_file + + with open(filename, 'w') as fout: + for k in range(0, self.n_group): + print("Points in Group " + str(k+1) + ": " + str(self.np_group[k])) + fout.write(f"{self.np_group[k]}\n") + for i in range(0, self.np_group[k]): + fout.write(f"{self.ip_group[k][i]}\n") + fout.write("\n") # empty line + + def add_points(self, added_points=np.array([]), group=None): + self.np_group[group] += len(added_points) + self.ip_group[group] = np.r_[self.ip_group[group], added_points] + + def get_max(self): + """get max prop value""" + + def get_ele(self, pval): + """return ele id with specified prop value""" + + +class source_sink(): + """class for handling all source/sink inputs""" + def __init__(self, source_dir=None, start_time_str='2000-01-01 00:00:00', timedeltas=[0.0, 86400.0*365], + source_eles=[], sink_eles=[], vsource_data=None, vsink_data=None): + + dummy_source_sink = TimeHistory( + file_name=None, start_time_str='2000-01-01 00:00:00', + data_array=np.c_[np.array(timedeltas), np.array(timedeltas)*0.0], columns=['datetime', '1'] + ) + + if source_dir is None: + nt = len(timedeltas) + nsources = len(source_eles) + nsinks = len(sink_eles) + + if vsource_data is None: + vsource_data = np.zeros([nt, nsources]) + if vsink_data is None: + vsink_data = np.zeros([nt, nsinks]) + + if nsources > 0: + self.vsource = TimeHistory(file_name=None, start_time_str=start_time_str, + data_array=np.c_[np.array(timedeltas), vsource_data], + columns=['datetime']+source_eles) + self.msource = TimeHistory(file_name=None, start_time_str=start_time_str, + data_array=np.c_[np.array(timedeltas), -9999*np.ones([nt, nsources]), np.zeros([nt, nsources])], + columns=['datetime']+source_eles+source_eles) + else: + source_eles = [1] # dummy with 0 vsource + self.vsource = copy.deepcopy(dummy_source_sink) + self.msource = TimeHistory(file_name=None, start_time_str=start_time_str, + data_array=np.c_[np.array(timedeltas), -9999*np.ones([nt, 1]), np.zeros([nt, 1])], + columns=['datetime', '1', '1']) + + if nsinks > 0: + self.vsink = TimeHistory(file_name=None, start_time_str=start_time_str, + data_array=np.c_[np.array(timedeltas), vsink_data], + columns=['datetime']+sink_eles) + else: + sink_eles = [1] # dummy with 0 vsource + self.vsink = copy.deepcopy(dummy_source_sink) + + self.source_sink_in = SourceSinkIn(filename=None, number_of_groups=2, ele_groups=[source_eles, sink_eles]) + + else: # read from existing files + self.source_sink_in = SourceSinkIn(f"{source_dir}/source_sink.in") + self.source_dir = source_dir + + if self.source_sink_in.np_group[0] > 0: + print('reading vsource\n') + self.vsource = TimeHistory(f"{source_dir}/vsource.th", start_time_str=start_time_str, columns=['datetime']+self.source_sink_in.ip_group[0].tolist()) + print('reading msource\n') + self.msource = TimeHistory(f"{source_dir}/msource.th", start_time_str=start_time_str, columns=['datetime']+self.source_sink_in.ip_group[0].tolist()+self.source_sink_in.ip_group[0].tolist()) + source_eles = self.source_sink_in.ip_group[0] + else: + self.vsource = copy.deepcopy(dummy_source_sink) + self.msource = TimeHistory(file_name=None, start_time_str=start_time_str, + data_array=np.c_[np.array(timedeltas), -9999*np.ones([nt, 1]), np.zeros([nt, 1])], + columns=['datetime', '1', '1']) + source_eles = [1] + + + if self.source_sink_in.np_group[1] > 0: + print('reading vsink\n') + self.vsink = TimeHistory(f"{source_dir}/vsink.th", start_time_str=start_time_str, columns=['datetime']+self.source_sink_in.ip_group[1].tolist()) + self.vsink.df.columns = self.vsink.df.columns.map(str) + sink_eles = self.source_sink_in.ip_group[1] + else: # dummy vsink + self.vsink = copy.deepcopy(dummy_source_sink) + sink_eles = [1] + + # accomodate dummy sources/sinks + self.source_sink_in = SourceSinkIn(filename=None, number_of_groups=2, ele_groups=[source_eles, sink_eles]) + + self.vsource.df.columns = self.vsource.df.columns.map(str) + self.msource.df.columns = self.msource.df.columns.map(str) + self.vsink.df.columns = self.vsink.df.columns.map(str) + + self.nsource = self.source_sink_in.np_group[0] + self.nsink = self.source_sink_in.np_group[1] + self.source_eles = source_eles + self.sink_eles = sink_eles + self.ntracers = 2 + + def update_vars(self): + self.vsource.df_propagate() + self.msource.df_propagate() + self.vsink.df_propagate() + self.nsource = self.source_sink_in.np_group[0] + self.nsink = self.source_sink_in.np_group[1] + self.source_eles = self.source_sink_in.ip_group[0] + self.sink_eles = self.source_sink_in.ip_group[1] + self.ntracers = 2 + + self.vsource.data = self.vsource.df.values[:, 1:] + + def __add__(self, other): + ''' + Add source_sink B to source_sink A, + retaining A's time stamps + ''' + A = copy.deepcopy(self) + B = other + for i, source_sink in enumerate(['vsource', 'vsink']): # 0: source; 1: sink + [_, new_eles_inds] = BinA(A.source_sink_in.ip_group[i], + B.source_sink_in.ip_group[i]) + [existing_eles_inds, _] = BinA(B.source_sink_in.ip_group[i], + A.source_sink_in.ip_group[i]) + new_eles = B.source_sink_in.ip_group[i][new_eles_inds] + existing_eles = A.source_sink_in.ip_group[i][existing_eles_inds] + + A.source_sink_in.add_points(added_points=new_eles, group=i) + + B_ss = eval(f"B.{source_sink}") + A_ss = eval(f"A.{source_sink}") + + f_interp = interpolate.interp1d(B_ss.time, B_ss.df.iloc[:, 1:].to_numpy().T) + B_df_interp = pd.DataFrame(data=f_interp(A_ss.time).T, columns=B_ss.df.columns[1:]) + + A_ss.df = copy.deepcopy(A_ss.df.join(B_df_interp[new_eles.astype(str)])) + A_ss.df[existing_eles.astype('str')] += B_df_interp[existing_eles.astype('str')] + + A.msource = TimeHistory(file_name=None, + data_array=np.c_[np.array([0.0, 86400*365*100]), + -9999*np.ones([2, len(A.vsource.df.columns)-1]), + np.zeros([2, len(A.vsource.df.columns)-1])], + columns=['datetime']+A.vsource.df.columns[1:].tolist()+A.vsource.df.columns[1:].tolist()) + A_ss.n_time = A_ss.df.shape[0] + A_ss.n_station = A_ss.df.shape[1] - 1 + + A.update_vars() + + return A + + def nc_writer(self, dirname=None): + if dirname is None: + raise Exception("dirname is required.") + os.makedirs(dirname, exist_ok=True) + + #breakpoint() + msource_data = self.msource.data.reshape([self.msource.n_time, self.ntracers, self.nsource]) + fout = Dataset(f'{dirname}/source.nc', 'w', format='NETCDF4') + fout.createDimension('nsources', self.nsource) + fout.createDimension('nsinks', self.nsink) + fout.createDimension('ntracers', self.ntracers) + fout.createDimension('time_msource', self.msource.n_time) + fout.createDimension('time_vsource', self.vsource.n_time) + fout.createDimension('time_vsink', self.vsink.n_time) + fout.createDimension('one', np.atleast_1d(1)) + + fout.createVariable('source_elem', 'i', ('nsources')) + fout['source_elem'][:] = self.source_sink_in.ip_group[0] + + fout.createVariable('vsource', 'f', ('time_vsource', 'nsources')) + fout['vsource'][:,:] = self.vsource.data + + fout.createVariable('msource', 'f', ('time_msource', 'ntracers', 'nsources')) + fout['msource'][:,:,:] = msource_data + + fout.createVariable('sink_elem', 'i', ('nsinks')) + fout['sink_elem'][:] = self.source_sink_in.ip_group[1] + + fout.createVariable('vsink', 'f', ('time_vsink', 'nsinks')) + fout['vsink'][:,:] = self.vsink.data + #C.vars.extend(['time_step_vsource','time_step_msource','time_step_vsink']) + + fout.createVariable('time_step_vsource', 'f', ('one')) + fout['time_step_vsource'][:] = self.vsource.delta_t + + fout.createVariable('time_step_msource', 'f', ('one')) + fout['time_step_msource'][:] = self.msource.delta_t + + fout.createVariable('time_step_vsink', 'f', ('one')) + fout['time_step_vsink'][:] = self.vsink.delta_t + + fout.close() diff --git a/pyschism/forcing/source_sink/sflux2source/TimeHistory.py b/pyschism/forcing/source_sink/sflux2source/TimeHistory.py new file mode 100644 index 00000000..aa9f4f93 --- /dev/null +++ b/pyschism/forcing/source_sink/sflux2source/TimeHistory.py @@ -0,0 +1,238 @@ +from datetime import datetime, timedelta +import pickle +import copy +from numbers import Integral, Real + +from pandas import read_csv, to_datetime, to_numeric +import pandas as pd +from matplotlib import pyplot +import numpy + +def running_mean1(X, N): + Y = X * numpy.nan + if (X.ndim > 1): # multi-dimensional + for i, x in enumerate(X.transpose()): + Y[:, i] = numpy.convolve(x, numpy.ones((N,))/N)[(N-1):] + else: # 1-d array + Y = numpy.convolve(X, numpy.ones((N,))/N)[(N-1):] + + return Y + +def running_mean(X, N): + cumsum = numpy.cumsum(numpy.insert(X, 0, 0)) + return (cumsum[N:] - cumsum[:-N]) / float(N) + + +class TimeHistory(): + + + """Class for manipulating *.th""" + + def __init__(self, file_name=None, start_time_str="2000-01-01 00:00:00", mask_val=None, sec_per_time_unit=1.0, data_array=None, columns=None): + + if mask_val is None: + mask_val = -9999999 + self.mask_val = mask_val + + """initialize """ + self.source_file = file_name + self.start_time_str = start_time_str + self.n_time = -99999999 + self.n_station = -99999999 + self.time = None + self.data = None + self.delta_t = -99999999.9 + self.datetime = None + self.sec_per_time_unit = sec_per_time_unit + + if file_name is None: + self.df = pd.DataFrame(data_array) + else: + self.df = read_csv(file_name, delim_whitespace=True, index_col=False, header=None) + if columns is not None: + self.df.columns = columns + + self.df.rename(columns={0: 'datetime'}, inplace=True) + self.df_propagate() # populate vars from dataframe + + def df_propagate(self): + [self.n_time, self.n_station] = list(self.df.shape) + self.n_station -= 1 # first col is time + if isinstance(self.df['datetime'][0], (Real, Integral)): + self.time = self.df['datetime'].to_numpy() + self.datetime = [datetime.strptime(self.start_time_str, '%Y-%m-%d %H:%M:%S') + + timedelta(seconds=x*self.sec_per_time_unit) for x in self.time] + self.df['datetime'] = self.datetime + elif isinstance(self.df['datetime'][0], datetime): + time_delta = self.df['datetime'] - self.df['datetime'][0] + self.time = numpy.array([x.total_seconds() for x in time_delta]) + self.datetime = self.df['datetime'] + else: + raise Exception('unknown datatype for the first column, neither datetime or float') + + + # mask invalid values + self.data = self.df.iloc[:, 1:].to_numpy(dtype=float) + self.data[abs(self.data-float(self.mask_val)) < 1e-5] = numpy.nan + print("Number of times: " + str(self.n_time)) + print("Number of stations: " + str(self.n_station)) + + self.delta_t = self.time[1] - self.time[0] + if self.delta_t < 50: # probably in days + self.time *= 86400.0 + self.delta_t *= 86400.0 + print("Time unit converted from days to seconds") + else: + print("Time unit: seconds") + + def get_running_mean(self, station_idx, window): + """ sort time series by time average""" + if station_idx is None: + valid_idx = range(0, self.n_station) + else: + valid_idx = station_idx + print(valid_idx) + + self.data = self.df.iloc[:, 1:].to_numpy(dtype=float) + self.data[abs(self.data-float(self.mask_val)) < 1e-5] = numpy.nan + data_rm = running_mean(self.data[:, valid_idx], window) + + return [data_rm] + + def get_time_average(self, station_idx, start_time_str=None, end_time_str=None): + if station_idx == []: + valid_idx = range(0, self.n_station) + else: + valid_idx = station_idx + print(valid_idx) + + if start_time_str is None: + idx1 = 0 + else: + idx1 = self.get_snapshot(start_time_str) + + if end_time_str is None: + idx2 = self.n_time + else: + idx2 = self.get_snapshot(end_time_str) + + self.data = self.df.iloc[:, 1:].to_numpy(dtype=float) + self.data[abs(self.data-float(self.mask_val)) < 1e-5] = numpy.nan + data_mean = numpy.mean(self.data[idx1:idx2, valid_idx], axis=0) + + return data_mean + + def sort_time_average(self, station_idx): + """ sort time series by time average""" + if station_idx == []: + valid_idx = range(0, self.n_station) + else: + valid_idx = station_idx + print(valid_idx) + + self.data = self.df.iloc[:, 1:].to_numpy(dtype=float) + self.data[abs(self.data-float(self.mask_val)) < 1e-5] = numpy.nan + data_mean = numpy.mean(self.data[:, valid_idx], axis=0) + id_sort = numpy.argsort(data_mean) + print("Preview on sorted: ") + print(data_mean[id_sort[-1:-11:-1]]) + + return [id_sort, data_mean] + + def get_snapshot(self, this_time_str=None, this_time_s_1970=None): + """ get a snapshot closest to the input time""" + if this_time_str is None: + raise Exception("needs this_time_str to get snapshot") + this_time_sec = timedelta(datetime.strftime(this_time_str, "%Y-%m-%d %H:%M:%S") + - self.datetime[0]).seconds + + t_idx = numpy.argmin(numpy.abs(self.time - this_time_sec)) + return(t_idx) + + def get_max_station(self): + """ get max among all stations """ + self.data = self.df.iloc[:, 1:].to_numpy(dtype=float) + self.data[abs(self.data-float(self.mask_val)) < 1e-5] = numpy.nan + data_sum = numpy.sum(self.data, axis=0) + max_idx = numpy.argmax(data_sum) + print("Station with the largest time-average value: " + str(max_idx)) + + return max_idx + + def simple_plot(self, idx, **args): + self.data = self.df.iloc[:, 1:].to_numpy(dtype=float) + self.data[abs(self.data-float(self.mask_val)) < 1e-5] = numpy.nan + pyplot.plot(self.datetime, self.data[:, idx], **args) + + def time_unit_conversion(self, unit_conv): + """convert time unit""" + """Col 0 is time: unit_conv from sec to day is 1/86400""" + self.time = self.time * unit_conv + + def data_unit_conversion(self, unit_conv): + """convert data unit""" + self.data = self.df.iloc[:, 1:].to_numpy(dtype=float) + self.data[abs(self.data-float(self.mask_val)) < 1e-5] = numpy.nan + self.data = self.data * unit_conv + + def export_subset(self, station_idx=None, time_idx=None, i_reset_time=False, subset_filename=None): + import os + """extract a subset from the original *.th""" + """by station_idx and time_idx""" + """[] idx means no operation""" + """if i_reset_time == True, then the subset *.th starts from 0 time""" + + if subset_filename is None: + subset_filename = os.path.abspath(self.source_file) + '.subset' + + if station_idx is None: + station_idx = numpy.array(range(self.n_station)) + else: + station_idx = numpy.array(station_idx) + + if time_idx is None: + time_idx = numpy.array(range(self.n_time)) + else: + time_idx = numpy.array(time_idx) + + self.data = self.df.iloc[:, 1:].to_numpy(dtype=float) + self.data[abs(self.data-float(self.mask_val)) < 1e-5] = numpy.nan + + sub_data = copy.deepcopy(self.data) + sub_time = copy.deepcopy(self.time) + if (len(station_idx) > 0): + sub_data = sub_data[:, station_idx] + if (len(time_idx) > 0): + sub_data = sub_data[time_idx, :] + sub_time = self.time[time_idx] + + if i_reset_time: + new_start_time_str = datetime.strftime(datetime.strptime(self.start_time_str, "%Y-%m-%d %H:%M:%S") +timedelta(seconds=sub_time[0]*self.sec_per_time_unit), "%Y-%m-%d %H:%M:%S") + sub_time = sub_time - sub_time[0] + else: + new_start_time_str = self.start_time_str + + with open(subset_filename, 'w') as fout: + for i, _ in enumerate(sub_time): + fout.write(str(sub_time[i]) + " " + + ' '.join(map(str, sub_data[i, :])) + + "\n") + + return TimeHistory(file_name=subset_filename, start_time_str=new_start_time_str, mask_val=self.mask_val, sec_per_time_unit=self.sec_per_time_unit, + data_array=None, columns=None) + + def writer(self, out_file_name): + self.df_propagate() + """ write self """ + with open(out_file_name, 'w') as fout: + for i, _ in enumerate(self.time): + fout.write(str(self.time[i]) + " " + + ' '.join(map(str, self.data[i, :])) + + "\n") + + @classmethod + def timestamp_microsecond(cls, epoch, utc_time): + """ convert date to microseconds """ + time_diff = utc_time - epoch + assert time_diff.resolution == timedelta(microseconds=1) + return (time_diff.days * 86400 + time_diff.seconds) * 10**6 + time_diff.microseconds From fbf346faffd99c094f98fc25d3376cb3dd78ce34 Mon Sep 17 00:00:00 2001 From: cuill Date: Wed, 30 Nov 2022 11:04:07 -0500 Subject: [PATCH 04/31] added dummy sinks when nsink=0 and changed regions in shaprio as optional --- pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py | 2 +- pyschism/mesh/gridgr3.py | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py b/pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py index 05136be2..ae9ca6ac 100644 --- a/pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py +++ b/pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py @@ -176,6 +176,7 @@ def __add__(self, other): A = copy.deepcopy(self) B = other for i, source_sink in enumerate(['vsource', 'vsink']): # 0: source; 1: sink + if i == 1 and len(A.source_sink_in.ip_group[i]) <= 1: continue [_, new_eles_inds] = BinA(A.source_sink_in.ip_group[i], B.source_sink_in.ip_group[i]) [existing_eles_inds, _] = BinA(B.source_sink_in.ip_group[i], @@ -211,7 +212,6 @@ def nc_writer(self, dirname=None): raise Exception("dirname is required.") os.makedirs(dirname, exist_ok=True) - #breakpoint() msource_data = self.msource.data.reshape([self.msource.n_time, self.ntracers, self.nsource]) fout = Dataset(f'{dirname}/source.nc', 'w', format='NETCDF4') fout.createDimension('nsources', self.nsource) diff --git a/pyschism/mesh/gridgr3.py b/pyschism/mesh/gridgr3.py index e0bf28fb..a6fe6ae0 100644 --- a/pyschism/mesh/gridgr3.py +++ b/pyschism/mesh/gridgr3.py @@ -112,7 +112,8 @@ def transform_to_cpp(coords, lonc, latc): class Shapiro(Gr3Field): @classmethod - def slope_filter(cls, hgrid, shapiro_vals1, depths, shapiro_max, threshold_slope, regions, shapiro_vals2, flags, lonc, latc): + def slope_filter(cls, hgrid, shapiro_vals1, depths, shapiro_max, threshold_slope, + regions=None, shapiro_vals2=None, flags=None, lonc=None, latc=None): """ https://github.com/wzhengui/pylibs/blob/1ee35efaa2d52fa682113126d84846ba33318f99/Utility/schism_file.py#L265 """ From 18b2d52edb815e02301462203ca2d1a9ae19ebd9 Mon Sep 17 00:00:00 2001 From: cuill Date: Fri, 2 Dec 2022 13:55:29 -0500 Subject: [PATCH 05/31] fixed negative vsource and added restart option for *.th.nc --- pyschism/forcing/hycom/hycom2schism.py | 290 +++++++++++------- .../source_sink/sflux2source/SourceSinkIn.py | 31 +- 2 files changed, 199 insertions(+), 122 deletions(-) diff --git a/pyschism/forcing/hycom/hycom2schism.py b/pyschism/forcing/hycom/hycom2schism.py index bceec245..1098ea22 100644 --- a/pyschism/forcing/hycom/hycom2schism.py +++ b/pyschism/forcing/hycom/hycom2schism.py @@ -177,7 +177,7 @@ def __init__(self, hgrid, vgrid=None): self.hgrid = hgrid self.vgrid = Vgrid.default() if vgrid is None else vgrid - def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D=True, TS=True, UV=True, adjust2D=False, lats=None, msl_shifts=None): + def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D=True, TS=True, UV=True, restart=False, adjust2D=False, lats=None, msl_shifts=None): outdir = pathlib.Path(outdir) self.start_date = start_date @@ -226,7 +226,7 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D= one=1 #ndt=np.zeros([ntimes]) - if elev2D: + if elev2D and restart == False: #timeseries_el=np.zeros([ntimes,NOP,nComp1]) #create netcdf dst_elev = Dataset(outdir / 'elev2D.th.nc', 'w', format='NETCDF4') @@ -246,8 +246,11 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D= dst_elev.createVariable('time_series', 'f', ('time', 'nOpenBndNodes', 'nLevels', 'nComponents')) #dst_elev['time_series'][:,:,:,:] = timeseries_el + elif restart: + dst_elev = Dataset(outdir / 'elev2D.th.nc', 'a', format='NETCDF4') + time_idx_restart = dst_elev['time'][:].shape[0] - if TS: + if TS and restart == False: #timeseries_s=np.zeros([ntimes,NOP,nvrt,nComp1]) dst_salt = Dataset(outdir / 'SAL_3D.th.nc', 'w', format='NETCDF4') #dimensions @@ -268,7 +271,7 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D= #temp #timeseries_t=np.zeros([ntimes,NOP,nvrt,nComp1]) - dst_temp = Dataset(outdir / 'TEM_3D.th.nc', 'w', format='NETCDF4') + dst_temp = Dataset(outdir / 'TEM_3D.th.nc', 'w', format='NETCDF4') #dimensions dst_temp.createDimension('nOpenBndNodes', NOP) dst_temp.createDimension('one', one) @@ -284,8 +287,11 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D= dst_temp.createVariable('time_series', 'f', ('time', 'nOpenBndNodes', 'nLevels', 'nComponents')) #dst_temp['time_series'][:,:,:,:] = timeseries_t + elif restart: + dst_salt = Dataset(outdir / 'SAL_3D.th.nc', 'a', format='NETCDF4') + dst_temp = Dataset(outdir / 'TEM_3D.th.nc', 'a', format='NETCDF4') - if UV: + if UV and restart == False: #timeseries_uv=np.zeros([ntimes,NOP,nvrt,nComp2]) dst_uv = Dataset(outdir / 'uv3D.th.nc', 'w', format='NETCDF4') #dimensions @@ -304,12 +310,25 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D= dst_uv.createVariable('time_series', 'f', ('time', 'nOpenBndNodes', 'nLevels', 'nComponents')) #dst_uv['time_series'][:,:,:,:] = timeseries_uv + elif restart: + dst_uv = Dataset(outdir / 'uv3D.th.nc', 'a', format='NETCDF4') + logger.info('**** Accessing GOFS data*****') t0=time() - for it, date in enumerate(self.timevector): + if restart == False: + timevector = self.timevector + it0 = 0 + elif restart: + timevector = self.timevector[time_idx_restart+1:] + it0 = time_idx_restart+1 + + for it1, date in enumerate(timevector): + + it = it0 + it1 + database=get_database(date) - logger.info(f'Fetching data for {date} from database {database}') + logger.info(f'Fetching data for {it}, {date} from database {database}') #loop over each open boundary ind1 = 0 @@ -527,7 +546,7 @@ def compute_nudge(lon, lat, nnode, opbd, out): return self.include - def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, rnday): + def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, rnday, restart = False): outdir = pathlib.Path(outdir) @@ -579,13 +598,60 @@ def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, one=1 ntimes=self.rnday+1 - timeseries_s=np.zeros([ntimes,nNode,nvrt,one]) - timeseries_t=np.zeros([ntimes,nNode,nvrt,one]) - ndt=np.zeros([ntimes]) + #timeseries_s=np.zeros([ntimes,nNode,nvrt,one]) + #timeseries_t=np.zeros([ntimes,nNode,nvrt,one]) + #ndt=np.zeros([ntimes]) + if restart: + dst_temp = Dataset(outdir / 'TEM_nu.nc', 'a', format='NETCDF4') + dst_salt = Dataset(outdir / 'SAL_nu.nc', 'a', format='NETCDF4') + time_idx_restart = dst_temp['time'][:].shape[0] + else: + dst_temp = Dataset(outdir / 'TEM_nu.nc', 'w', format='NETCDF4') + #dimensions + dst_temp.createDimension('node', nNode) + dst_temp.createDimension('nLevels', nvrt) + dst_temp.createDimension('one', one) + dst_temp.createDimension('time', None) + #variables + dst_temp.createVariable('time', 'f', ('time',)) + #dst_temp['time'][:] = ndt + + dst_temp.createVariable('map_to_global_node', 'i4', ('node',)) + dst_temp['map_to_global_node'][:] = include+0 + + dst_temp.createVariable('tracer_concentration', 'f', ('time', 'node', 'nLevels', 'one')) + #dst_temp['tracer_concentration'][:,:,:,:] = timeseries_t + + #salinity + dst_salt = Dataset(outdir / 'SAL_nu.nc', 'w', format='NETCDF4') + #dimensions + dst_salt.createDimension('node', nNode) + dst_salt.createDimension('nLevels', nvrt) + dst_salt.createDimension('one', one) + dst_salt.createDimension('time', None) + #variables + dst_salt.createVariable('time', 'f', ('time',)) + #dst_salt['time'][:] = ndt + + dst_salt.createVariable('map_to_global_node', 'i4', ('node',)) + dst_salt['map_to_global_node'][:] = include+0 + + dst_salt.createVariable('tracer_concentration', 'f', ('time', 'node', 'nLevels', 'one')) + #dst_salt['tracer_concentration'][:,:,:,:] = timeseries_s + logger.info('**** Accessing GOFS data*****') + if restart: + timevector = self.timevector[time_idx_restart+1:] + it0 = time_idx_restart + else: + timevector = self.timevector + it0 = 0 + t0=time() - for it, date in enumerate(self.timevector): + for it1, date in enumerate(timevector): + + it = it0 + it1 database=get_database(date) logger.info(f'Fetching data for {date} from database {database}') @@ -632,52 +698,26 @@ def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, logger.info('****Interpolation starts****') - ndt[it]=it + #ndt[it]=it #salt + dst_salt['time'][it] = it salt_int = interp_to_points_3d(dep, y2, x2, bxyz, salt) salt_int = salt_int.reshape(zcor2.shape) - timeseries_s[it,:,:,0]=salt_int + #timeseries_s[it,:,:,0]=salt_int + dst_salt['tracer_concentration'][it,:,:,0] = salt_int #temp + dst_temp['time'][it] = it temp_int = interp_to_points_3d(dep, y2, x2, bxyz, ptemp) temp_int = temp_int.reshape(zcor2.shape) - timeseries_t[it,:,:,0]=temp_int + #timeseries_t[it,:,:,0]=temp_int + dst_temp['tracer_concentration'][it,:,:,0] = temp_int ds.close() - with Dataset(outdir / 'TEM_nu.nc', 'w', format='NETCDF4') as dst: - #dimensions - dst.createDimension('node', nNode) - dst.createDimension('nLevels', nvrt) - dst.createDimension('one', one) - dst.createDimension('time', None) - #variables - dst.createVariable('time', 'f', ('time',)) - dst['time'][:] = ndt - - dst.createVariable('map_to_global_node', 'i4', ('node',)) - dst['map_to_global_node'][:] = include+1 - - dst.createVariable('tracer_concentration', 'f', ('time', 'node', 'nLevels', 'one')) - dst['tracer_concentration'][:,:,:,:] = timeseries_t - - with Dataset(outdir / 'SAL_nu.nc', 'w', format='NETCDF4') as dst: - #dimensions - dst.createDimension('node', nNode) - dst.createDimension('nLevels', nvrt) - dst.createDimension('one', one) - dst.createDimension('time', None) - #variables - dst.createVariable('time', 'f', ('time',)) - dst['time'][:] = ndt - - dst.createVariable('map_to_global_node', 'i4', ('node',)) - dst['map_to_global_node'][:] = include+1 - - dst.createVariable('tracer_concentration', 'f', ('time', 'node', 'nLevels', 'one')) - dst['tracer_concentration'][:,:,:,:] = timeseries_s - - + #dst_temp.close() + #dst_salt.close() + logger.info(f'Writing *_nu.nc takes {time()-t0} seconds') class DownloadHycom: @@ -690,76 +730,96 @@ def __init__(self, hgrid): xmax = xmax + 360. if xmax < 0 else xmax self.bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) - def fetch_data(self, date): + def fetch_data(self, date, fmt='schism'): + ''' + date: datetime.datetime + fmt: 'schism' - for Fortran code; 'hycom' - raw netCDF from HYCOM + ''' database=get_database(date) logger.info(f'Fetching data for {date} from database {database}') time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 = get_idxs(date, database, self.bbox) - url_ssh = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ - f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \ - f'surf_el[{time_idx}][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' - foutname = f'SSH_{date.strftime("%Y%m%d")}.nc' - logger.info(f'filename is {foutname}') - ds = xr.open_dataset(url_ssh) - ds = convert_longitude(ds) - ds = ds.rename_dims({'lon':'xlon'}) - ds = ds.rename_dims({'lat':'ylat'}) - ds = ds.rename_vars({'lat':'ylat'}) - ds = ds.rename_vars({'lon':'xlon'}) - ds.to_netcdf(foutname, 'w', 'NETCDF3_CLASSIC', unlimited_dims='time') - ds.close() - - url_uv = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ - f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \ - f'water_u[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ - f'water_v[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' - - foutname = f'UV_{date.strftime("%Y%m%d")}.nc' - logger.info(f'filename is {foutname}') - ds = xr.open_dataset(url_uv) - ds = convert_longitude(ds) - ds = ds.rename_dims({'lon':'xlon'}) - ds = ds.rename_dims({'lat':'ylat'}) - ds = ds.rename_vars({'lat':'ylat'}) - ds = ds.rename_vars({'lon':'xlon'}) - ds.to_netcdf(foutname, 'w', 'NETCDF3_CLASSIC', unlimited_dims='time') - ds.close() - - url_ts = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ - f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \ - f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ - f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' - - foutname = f'ST_{date.strftime("%Y%m%d")}.nc' - logger.info(f'filename is {foutname}') - - ds = xr.open_dataset(url_ts) - - temp = ds.water_temp.values - salt = ds.salinity.values - dep = ds.depth.values - - #convert in-situ temperature to potential temperature - ptemp = ConvertTemp(salt, temp, dep) - - #drop water_temp variable and add new temperature variable - ds = ds.drop('water_temp') - ds['temperature']=(['time','depth','lat','lon'], ptemp) - ds.temperature.attrs = { - 'long_name': 'Sea water potential temperature', - 'standard_name': 'sea_water_potential_temperature', - 'units': 'degC' - } - - #ds.assign(water_temp2=ptemp) - #ds.assign.attrs = ds.water_temp.attrs - - ds = convert_longitude(ds) - ds = ds.rename_dims({'lon':'xlon'}) - ds = ds.rename_dims({'lat':'ylat'}) - ds = ds.rename_vars({'lat':'ylat'}) - ds = ds.rename_vars({'lon':'xlon'}) - ds.to_netcdf(foutname, 'w', unlimited_dims='time', encoding={'temperature':{'dtype': 'h', '_FillValue': -30000.,'scale_factor': 0.001, 'add_offset': 20., 'missing_value': -30000.}}) - ds.close() + if fmt == 'schism': + url_ssh = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ + f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \ + f'surf_el[{time_idx}][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' + foutname = f'SSH_{date.strftime("%Y%m%d")}.nc' + logger.info(f'filename is {foutname}') + ds = xr.open_dataset(url_ssh) + ds = convert_longitude(ds) + ds = ds.rename_dims({'lon':'xlon'}) + ds = ds.rename_dims({'lat':'ylat'}) + ds = ds.rename_vars({'lat':'ylat'}) + ds = ds.rename_vars({'lon':'xlon'}) + ds.to_netcdf(foutname, 'w', 'NETCDF3_CLASSIC', unlimited_dims='time') + ds.close() + + url_uv = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ + f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \ + f'water_u[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ + f'water_v[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' + + foutname = f'UV_{date.strftime("%Y%m%d")}.nc' + logger.info(f'filename is {foutname}') + ds = xr.open_dataset(url_uv) + ds = convert_longitude(ds) + ds = ds.rename_dims({'lon':'xlon'}) + ds = ds.rename_dims({'lat':'ylat'}) + ds = ds.rename_vars({'lat':'ylat'}) + ds = ds.rename_vars({'lon':'xlon'}) + ds.to_netcdf(foutname, 'w', 'NETCDF3_CLASSIC', unlimited_dims='time') + ds.close() + + url_ts = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ + f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \ + f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ + f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' + + foutname = f'ST_{date.strftime("%Y%m%d")}.nc' + logger.info(f'filename is {foutname}') + + ds = xr.open_dataset(url_ts) + + temp = ds.water_temp.values + salt = ds.salinity.values + dep = ds.depth.values + + #convert in-situ temperature to potential temperature + ptemp = ConvertTemp(salt, temp, dep) + + #drop water_temp variable and add new temperature variable + ds = ds.drop('water_temp') + ds['temperature']=(['time','depth','lat','lon'], ptemp) + ds.temperature.attrs = { + 'long_name': 'Sea water potential temperature', + 'standard_name': 'sea_water_potential_temperature', + 'units': 'degC' + } + + #ds.assign(water_temp2=ptemp) + #ds.assign.attrs = ds.water_temp.attrs + + ds = convert_longitude(ds) + ds = ds.rename_dims({'lon':'xlon'}) + ds = ds.rename_dims({'lat':'ylat'}) + ds = ds.rename_vars({'lat':'ylat'}) + ds = ds.rename_vars({'lon':'xlon'}) + ds.to_netcdf(foutname, 'w', unlimited_dims='time', encoding={'temperature':{'dtype': 'h', '_FillValue': -30000.,'scale_factor': 0.001, 'add_offset': 20., 'missing_value': -30000.}}) + ds.close() + + elif fmt == 'hycom': + url=f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ + f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \ + f'surf_el[{time_idx}][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ + f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ + f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ + f'water_u[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ + f'water_v[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' + + foutname = f'hycom_{date.strftime("%Y%m%d")}.nc' + ds = xr.open_dataset(url) + ds.to_netcdf(foutname, 'w') + + ds.close() diff --git a/pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py b/pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py index ae9ca6ac..cf11331a 100644 --- a/pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py +++ b/pyschism/forcing/source_sink/sflux2source/SourceSinkIn.py @@ -176,6 +176,8 @@ def __add__(self, other): A = copy.deepcopy(self) B = other for i, source_sink in enumerate(['vsource', 'vsink']): # 0: source; 1: sink + print(f'Adding type {i}:') + #if vsink only contains dummy sink, then skip interpolation (todo: remove dummy sinks) if i == 1 and len(A.source_sink_in.ip_group[i]) <= 1: continue [_, new_eles_inds] = BinA(A.source_sink_in.ip_group[i], B.source_sink_in.ip_group[i]) @@ -189,24 +191,39 @@ def __add__(self, other): B_ss = eval(f"B.{source_sink}") A_ss = eval(f"A.{source_sink}") - f_interp = interpolate.interp1d(B_ss.time, B_ss.df.iloc[:, 1:].to_numpy().T) + #use A's time orgin + time_diff = (B_ss.df.iloc[0, 0] - A_ss.df.iloc[0, 0]).total_seconds() + + f_interp = interpolate.interp1d(B_ss.time + time_diff, B_ss.df.iloc[:, 1:].to_numpy().T) B_df_interp = pd.DataFrame(data=f_interp(A_ss.time).T, columns=B_ss.df.columns[1:]) + # clip the values in case of small negative/positive values due to trunction error during interpolation + + if i == 0: # source + B_df_interp = B_df_interp.clip(0.0, None) + else: + B_df_interp = B_df_interp.clip(None, 0.0) + A_ss.df = copy.deepcopy(A_ss.df.join(B_df_interp[new_eles.astype(str)])) A_ss.df[existing_eles.astype('str')] += B_df_interp[existing_eles.astype('str')] - A.msource = TimeHistory(file_name=None, - data_array=np.c_[np.array([0.0, 86400*365*100]), - -9999*np.ones([2, len(A.vsource.df.columns)-1]), - np.zeros([2, len(A.vsource.df.columns)-1])], - columns=['datetime']+A.vsource.df.columns[1:].tolist()+A.vsource.df.columns[1:].tolist()) A_ss.n_time = A_ss.df.shape[0] A_ss.n_station = A_ss.df.shape[1] - 1 - A.update_vars() + A.msource = TimeHistory( + file_name=None, + data_array=np.c_[np.array([0.0, 86400*365*100]), -9999*np.ones([2, len(A.vsource.df.columns)-1]), np.zeros([2, len(A.vsource.df.columns)-1])], + columns=['datetime']+A.vsource.df.columns[1:].tolist()+A.vsource.df.columns[1:].tolist()) + + A.update_vars() return A + def writer(self, dirname=None): + if dirname is None: + raise Exception("dirname is required.") + os.makedirs(dirname, exist_ok=True) + def nc_writer(self, dirname=None): if dirname is None: raise Exception("dirname is required.") From 2c6602aeb66a5d0990d9d45c99a6dcda07706929 Mon Sep 17 00:00:00 2001 From: Jaime R Calzada Date: Thu, 15 Dec 2022 11:29:39 -0500 Subject: [PATCH 06/31] merge improvements + alias name fixes --- pyschism/forcing/bctides/bctides.py | 395 ++++++++-------------------- pyschism/forcing/bctides/iettype.py | 58 +++- pyschism/forcing/bctides/ifltype.py | 37 ++- pyschism/forcing/bctides/isatype.py | 34 ++- pyschism/forcing/bctides/itetype.py | 23 +- pyschism/forcing/bctides/itrtype.py | 18 +- 6 files changed, 253 insertions(+), 312 deletions(-) diff --git a/pyschism/forcing/bctides/bctides.py b/pyschism/forcing/bctides/bctides.py index 49796578..8ee4f027 100644 --- a/pyschism/forcing/bctides/bctides.py +++ b/pyschism/forcing/bctides/bctides.py @@ -1,10 +1,12 @@ from datetime import datetime, timedelta +from functools import cached_property import pathlib from typing import Dict, Union import logging -from pyschism import dates +import numpy as np +from pyschism import dates from pyschism.mesh.vgrid import Vgrid from pyschism.forcing.bctides import iettype, ifltype, isatype, itetype, itrtype, Tides from pyschism.forcing.bctides.elev2d import Elev2D @@ -26,18 +28,12 @@ def __set__(self, obj, val): if val is not None: if isinstance(val, dict): for bnd_id, ibctype in val.items(): - if not isinstance(val, (self.bctype, type(None))): + if not isinstance(ibctype, (self.bctype, type(None))): raise TypeError( - f"Argument {val} must be of type {self.bctype} " + f"Argument {ibctype} for boundary {bnd_id} must be of type {self.bctype} " f" or None, not type {type(ibctype)}." ) - # TODO - raise NotImplementedError("Need to find the column name") - idxs = obj.gdf[ - (obj.gdf["id"] == bnd_id & obj.gdf["id"] == bnd_id) - ].index.values - for idx in idxs: - obj.gdf.at[idx, val] = obj + obj.gdf.at[np.where(obj.gdf["id"] == bnd_id)[0][0], self.bctype.__name__.lower()] = ibctype else: if not isinstance(val, self.bctype): raise TypeError( @@ -65,6 +61,7 @@ class Bctides(metaclass=BctidesMeta): start_date = dates.StartDate() end_date = dates.EndDate() + rnday = dates.RunDays() def __init__( self, @@ -152,6 +149,25 @@ def write( raise IOError("path exists and overwrite is False") with open(bctides, "w") as f: f.write(str(self)) + # write nudge + for bctype, tracer in {"itetype": "TEM", "isatype": "SAL"}.items(): + for boundary in self.gdf.itertuples(): + data_source = getattr(boundary, bctype) + if data_source is not None: + if hasattr(data_source, 'rlmax'): + # I admit this exec is hacky. + # pros: works well, it's simple, we don't need a return value + # cons: might be confusing to read. + # This generates all the nudges and writes the nudge files. + exec( + f"from pyschism.forcing.bctides.nudge import {tracer}_Nudge;" + f"_tracer = output_directory / f'{tracer}_nudge.gr3' if {tracer.lower()}3D is True else {tracer};" + f"_tr={tracer}_Nudge(self, data_source, rlmax=data_source.rlmax, rnu_day=data_source.rnu_day);" + f'logger.info(f"Writing {tracer} nudge to file ' + + r'{_tracer}");' + "_tr.write(_tracer, overwrite=overwrite)" + ) + break ## write nudge #for bctype, tracer in {"itetype": "TEM", "isatype": "SAL"}.items(): # for boundary in self.gdf.itertuples(): @@ -171,6 +187,25 @@ def write( # "_tr.write(_tracer, overwrite=overwrite)" # ) # break + # write nudge + for bctype, tracer in {"itetype": "TEM", "isatype": "SAL"}.items(): + for boundary in self.gdf.itertuples(): + data_source = getattr(boundary, bctype) + if data_source is not None: + if hasattr(data_source, 'rlmax'): + # I admit this exec is hacky. + # pros: works well, it's simple, we don't need a return value + # cons: might be confusing to read. + # This generates all the nudges and writes the nudge files. + exec( + f"from pyschism.forcing.bctides.nudge import {tracer}_Nudge;" + f"_tracer = output_directory / f'{tracer}_nudge.gr3' if {tracer.lower()}3D is True else {tracer};" + f"_tr={tracer}_Nudge(self, data_source, rlmax=data_source.rlmax, rnu_day=data_source.rnu_day);" + f'logger.info(f"Writing {tracer} nudge to file ' + + r'{_tracer}");' + "_tr.write(_tracer, overwrite=overwrite)" + ) + break def write_elev2D(): _elev2D = output_directory / "elev2D.th.nc" if elev2D is True else elev2D @@ -259,8 +294,7 @@ def get_focing_digit(bctype): if bctype is not None: # sensitive to MRO. return str( - getattr( - bctype, f"{bctype.__class__.__bases__[0].__name__.lower()}") + getattr(bctype, f"{bctype.__class__.__bases__[0].__name__.lower()}") ) return "0" @@ -279,16 +313,15 @@ def get_focing_digit(bctype): ) return "\n".join(f) - @property + @cached_property def gdf(self): - if not hasattr(self, "_gdf"): - self._gdf = self.hgrid.boundaries.open.copy() - self._gdf["iettype"] = None - self._gdf["ifltype"] = None - self._gdf["isatype"] = None - self._gdf["itetype"] = None - self._gdf["itrtype"] = None - return self._gdf + gdf = self.hgrid.boundaries.open.copy() + gdf["iettype"] = None + gdf["ifltype"] = None + gdf["isatype"] = None + gdf["itetype"] = None + gdf["itrtype"] = None + return gdf @property def ntip(self): @@ -298,266 +331,68 @@ def ntip(self): def nbfr(self): return len(self.tides.get_active_forcing_constituents()) - @property - def rnday(self): - return self.end_date - self.start_date - - @property + @cached_property def tides(self): + return TidalConstituentCombiner(self.gdf) + + +class TidalConstituentCombiner(Tides): + def __init__(self, gdf): + self.gdf = gdf + afc = self.get_active_forcing_constituents() + apc = self.get_active_potential_constituents() + for constituent in set([*afc, *apc]): + self.use_constituent( + constituent, + forcing=True if constituent in afc else False, + potential=True if constituent in apc else False, + ) - if not hasattr(self, "_tides"): - - class TidalConstituentCombiner(Tides): - def __init__(self, gdf): - self.gdf = gdf - afc = self.get_active_forcing_constituents() - apc = self.get_active_potential_constituents() - for constituent in set([*afc, *apc]): - self.use_constituent( - constituent, - forcing=True if constituent in afc else False, - potential=True if constituent in apc else False, - ) - - def get_active_forcing_constituents(self): - active_constituents = set() - for row in self.gdf.itertuples(): - if row.iettype is not None: - if row.iettype.iettype in [3, 5]: - [ - active_constituents.add(x) - for x in row.iettype.tides.get_active_constituents() - ] - if row.ifltype is not None: - if row.ifltype.ifltype in [3, 5]: - [ - active_constituents.add(x) - for x in row.ifltype.tides.get_active_constituents() - ] - - return list(active_constituents) - - def get_active_potential_constituents(self): - active_constituents = set() - for row in self.gdf.itertuples(): - if row.iettype is not None: - if row.iettype.iettype in [3, 5]: - [ - active_constituents.add(x) - for x in row.iettype.tides.get_active_potential_constituents() - ] - if row.ifltype is not None: - if row.ifltype.ifltype in [3, 5]: - [ - active_constituents.add(x) - for x in row.ifltype.tides.get_active_potential_constituents() - ] - - return list(active_constituents) - - @property - def constituents(self): - if not hasattr(self, "_constituents"): - self._constituents = sorted( - list( - set( - [ - *self.get_active_potential_constituents(), - *self.get_active_forcing_constituents(), - ] - ) - ) - ) - return self._constituents - - self._tides = TidalConstituentCombiner(self.gdf) - - return self._tides - - -class TidesCombiner(Tides): - def __init__(self, bctides): - self.bctides = bctides + def get_active_forcing_constituents(self): + active_constituents = set() + for row in self.gdf.itertuples(): + if row.iettype is not None: + if row.iettype.iettype in [3, 5]: + [ + active_constituents.add(x) + for x in row.iettype.tides.get_active_constituents() + ] + if row.ifltype is not None: + if row.ifltype.ifltype in [3, 5]: + [ + active_constituents.add(x) + for x in row.ifltype.tides.get_active_constituents() + ] + + return list(active_constituents) def get_active_potential_constituents(self): - const = dict() - # for row in self.bctides.gdf.itertuples(): - - # forcing = data['forcing'] - # if isinstance(forcing, Tides): - # for active in forcing.get_active_potential_constituents(): - # const[active] = True - # return tuple(const.keys()) - - # def get_active_forcing_constituents(self): - # # set active tidal forcing constituents - # const = dict() - # for id, data in self._model_domain.open_boundaries: - # forcing = data['forcing'] - # if isinstance(forcing, Tides): - # for active in forcing.get_active_forcing_constituents(): - # const[active] = True - # return tuple(const.keys()) - - -# ----------- draft -# @property -# def tides(self): -# if not hasattr(self, '_tides'): -# # get the first one you can find, since the Tides object is a -# # singleton. -# tides = None -# for boundary in self.hgrid.boundaries.open.itertuples(): -# if boundary.iettype is not None: -# if hasattr(boundary.iettype, "tides"): -# tides = boundary.iettype.tides -# break -# elif boundary.ifltype is not None: -# if hasattr(boundary.ifltype, "tides"): -# tides = boundary.ifltype.tides -# break -# self._tides = tides -# return self._tides - -# @property -# def tracers(self) -> List[Dict[Any, Union[bctides.itrtype.Itrtype, None]]]: -# # if not hasattr(self, '_tracers'): -# # # tracers: List[Dict[Any, Union[itrtype.Itrtype, None]]] = [] -# # boundary_data = {} -# # for boundary in self.hgrid.boundaries.open.itertuples(): -# # itrtypes = boundary.itrtype -# # if itrtypes is None: -# # tracers.append({}) -# # for tracer in boundary.itr -# # tracers.append() -# # tracer.setdefault( - -# # ) - -# # _itrtype = boundary.itrtype -# # return self._tracers -# # TODO: Cheating for now... -# return [] - -# @property -# def ntip(self): -# if self.tides is None: -# return 0 -# return len(self.tides.get_active_potential_constituents()) - -# @property -# def nbfr(self): -# if self.tides is None: -# return 0 -# return self.tides.nbfr - -# @property -# def Z0(self): -# if hasattr(self.tides, '_Z0'): -# return self.tides._Z0 - -# @Z0.setter -# def Z0(self, Z0): -# self.tides.add_Z0(Z0) - -# @property -# def cutoff_depth(self): -# return self._cutoff_depth - -# @cutoff_depth.setter -# def cutoff_depth(self, cutoff_depth: float): -# self._cutoff_depth = float(cutoff_depth) - -# @property -# def subtidal_database(self): -# return self._subtidal_database - -# @subtidal_database.setter -# def subtidal_database(self, subtidal_database: SubTidalDatabase): -# if subtidal_database is not None: -# # self._subtidal_database = Tides(subtidal_database=subtidal_database) -# else: -# self._subtidal_database = None - -# @property -# def elevation(self): -# return self._elevation - -# @elevation.setter -# def elevation(self, elevation): -# if elevation is not None: -# assert isinstance(elevation, iettype.Iettype) -# self._elevation = elevation - -# @property -# def velocity(self): -# return self._velocity - -# @velocity.setter -# def velocity(self, velocity): -# if velocity is not None: -# assert isinstance(velocity, ifltype.Ifltype) -# self._velocity = velocity - -# @property -# def temperature(self): -# return self._temperature - -# @temperature.setter -# def temperature(self, temperature: Union[itetype.Itetype, None]): -# if temperature is not None: -# assert isinstance(temperature, itetype.Itetype) -# self._temperature = temperature - -# @property -# def salinity(self): -# return self._salinity - -# @salinity.setter -# def salinity(self, salinity: Union[isatype.Isatype, None]): -# if salinity is not None: -# assert isinstance(salinity, isatype.Isatype) -# self._salinity = salinity - - -# class HgridDescriptor: - -# def __set__(self, obj, val: Hgrid): -# if not isinstance(val, Hgrid): -# raise TypeError( -# f'Argument hgrid must be of type {Hgrid}, not type ' -# f'{type(val)}.') -# obj.__dict__['hgrid'] = val - -# def __get__(self, obj, val): -# return obj.__dict__['hgrid'] - - -# class StartDateDescriptor: - -# def __set__(self, obj, val: datetime): -# if not isinstance(val, datetime): -# raise TypeError( -# f'Argument start_date must be of type {datetime}, ' -# f'not type {type(val)}.') -# if datetime_is_naive(val): -# val = pytz.timezone('UTC').localize(val) -# obj.__dict__['start_date'] = val - -# def __get__(self, obj, val): -# return obj.__dict__['start_date'] - - -# class RndayDescriptor: - -# def __set__(self, obj, val: Union[int, float, timedelta]): -# if not isinstance(val, (int, float, timedelta)): -# raise TypeError( -# f'Argument rnday must be of type {int}, {float} or ' -# f'{timedelta}, not type {type(val)}.') -# if not isinstance(val, timedelta): -# val = timedelta(days=val) -# obj.__dict__['rnday'] = val - -# def __get__(self, obj, val) -> timedelta: -# return obj.__dict__['rnday'] + active_constituents = set() + for row in self.gdf.itertuples(): + if row.iettype is not None: + if row.iettype.iettype in [3, 5]: + [ + active_constituents.add(x) + for x in row.iettype.tides.get_active_potential_constituents() + ] + if row.ifltype is not None: + if row.ifltype.ifltype in [3, 5]: + [ + active_constituents.add(x) + for x in row.ifltype.tides.get_active_potential_constituents() + ] + + return list(active_constituents) + + @cached_property + def constituents(self): + return sorted( + list( + set( + [ + *self.get_active_potential_constituents(), + *self.get_active_forcing_constituents(), + ] + ) + ) + ) diff --git a/pyschism/forcing/bctides/iettype.py b/pyschism/forcing/bctides/iettype.py index 1b70f071..86f3a85b 100644 --- a/pyschism/forcing/bctides/iettype.py +++ b/pyschism/forcing/bctides/iettype.py @@ -40,15 +40,34 @@ def iettype(self) -> int: class TidalElevation(Iettype): - def __init__(self, constituents="all", database="hamtide"): + def __init__(self, constituents='all', database="hamtide"): self.tides = Tides(tidal_database=database, constituents=constituents) + + def add_constituent( + self, + name: str, + amplitude, + angular_frequency, + phase=0., + **kwargs, + ): + self.tides.add_constituent( + name=name, + angular_frequency=angular_frequency, + elevation_amplitude=amplitude, + elevation_phase=phase, + velocity_amplitude=self.tides._amplitudes['velocity'].get(name, 0.), + velocity_phase=self.tides._phases['velocity'].get(name, 0.), + **kwargs + ) def get_boundary_string(self, hgrid, boundary, global_constituents: List[str] = None): f = [] + crs = 'epsg:4326' if hgrid.crs is not None else None if global_constituents is None: for constituent in self.tides.get_active_forcing_constituents(): f.append(f"{constituent}") - vertices = hgrid.get_xy(crs="EPSG:4326")[boundary.indexes, :] + vertices = hgrid.get_xy(crs=crs)[boundary.indexes, :] amp, phase = self.tides.get_elevation(constituent, vertices) for i in range(len(boundary.indexes)): f.append(f"{amp[i]:.8e} {phase[i]:.8e}") @@ -60,7 +79,7 @@ def get_boundary_string(self, hgrid, boundary, global_constituents: List[str] = for i in range(len(boundary.indexes)): f.append(f"{0:.8e} {0:.8e}") else: - vertices = hgrid.get_xy(crs="EPSG:4326")[boundary.indexes, :] + vertices = hgrid.get_xy(crs=crs)[boundary.indexes, :] amp, phase = self.tides.get_elevation(constituent, vertices) for i in range(len(boundary.indexes)): f.append(f"{amp[i]:.8e} {phase[i]:.8e}") @@ -72,6 +91,11 @@ def iettype(self): return 3 +class HarmonicElevation(TidalElevation): + + def __init__(self): + self.tides = Tides(constituents=None) + class SpatiallyVaryingTimeHistoryElevation(Iettype): class BaroclinicDatabases(Enum): #RTOFS = hycom.RTOFS @@ -139,9 +163,25 @@ def iettype(self): return -1 -Iettype1 = UniformTimeHistoryElevation -Iettype2 = ConstantElevation -Iettype3 = TidalElevation -Iettype4 = SpatiallyVaryingTimeHistoryElevation -Iettype5 = TidalAndSpatiallyVaryingElevationTimeHistory -Iettype_1 = ZeroElevation +class Iettype1(UniformTimeHistoryElevation): + pass + + +class Iettype2(ConstantElevation): + pass + + +class Iettype3(TidalElevation): + pass + + +class Iettype4(SpatiallyVaryingTimeHistoryElevation): + pass + + +class Iettype5(TidalAndSpatiallyVaryingElevationTimeHistory): + pass + + +class Iettype_1(ZeroElevation): + pass diff --git a/pyschism/forcing/bctides/ifltype.py b/pyschism/forcing/bctides/ifltype.py index 46e99e15..46c5583a 100644 --- a/pyschism/forcing/bctides/ifltype.py +++ b/pyschism/forcing/bctides/ifltype.py @@ -25,13 +25,15 @@ def ifltype(self): class ConstantVelocity(Ifltype): def __init__(self, value): - raise NotImplementedError(f"{self.__class__.__name__}") self.value = value @property def ifltype(self) -> int: return 2 + def get_boundary_string(self, *args, **kwargs): + return f'{self.value}' + class TidalVelocity(Ifltype): def __init__(self, constituents="all", database="hamtide"): @@ -151,10 +153,29 @@ def ifltype(self): return -5 -Ifltype1 = UniformTimeHistoryVelocity -Ifltype2 = ConstantVelocity -Ifltype3 = TidalVelocity -Ifltype4 = SpatiallyVaryingTimeHistoryVelocity -Ifltype5 = TidalAndSpatiallyVaryingVelocityTimeHistory -Ifltype_1 = ZeroVelocity -Flather = ZeroVelocity +class Ifltype1(UniformTimeHistoryVelocity): + pass + + +class Ifltype2(ConstantVelocity): + pass + + +class Ifltype3(TidalVelocity): + pass + + +class Ifltype4(SpatiallyVaryingTimeHistoryVelocity): + pass + + +class Ifltype5(TidalAndSpatiallyVaryingVelocityTimeHistory): + pass + + +class Ifltype_1(ZeroVelocity): + pass + + +class Flather(ZeroVelocity): + pass diff --git a/pyschism/forcing/bctides/isatype.py b/pyschism/forcing/bctides/isatype.py index 4c6a09ba..9f5cb54f 100644 --- a/pyschism/forcing/bctides/isatype.py +++ b/pyschism/forcing/bctides/isatype.py @@ -1,5 +1,6 @@ from abc import abstractmethod from enum import Enum +from matplotlib.pyplot import new_figure_manager from pyschism.forcing.bctides.bctypes import Bctype from pyschism.forcing import hycom @@ -23,14 +24,25 @@ def isatype(self) -> int: class ConstantSalinity(Isatype): - def __init__(self, value): - raise NotImplementedError(f"{self.__class__.__name__}") + def __init__(self, value: float, nudging_factor: float): self.value = value + if not (nudging_factor >= 0) and (nudging_factor <= 1): + raise ValueError(f'Argument `nudging_factor` must be get 0 and let 1, but got {nudging_factor}.') + self.nudging_factor = nudging_factor + super().__init__() + @property def isatype(self) -> int: return 2 + def get_boundary_string(self, *args, **kwargs) -> str: + boundary_string = [ + f'{self.value:0.6f}', + f'{self.nudging_factor:0.6f}', + ] + return '\n'.join(boundary_string) + class SalinityInitialConditions(Isatype): def __init__(self): @@ -74,7 +86,17 @@ def isatype(self): return 4 -Isatype1 = UniformTimeHistorySalinity -Isatype2 = ConstantSalinity -Isatype3 = SalinityInitialConditions -Isatype4 = SpatiallyVaryingTimeHistorySalinity +class Isatype1(UniformTimeHistorySalinity): + pass + + +class Isatype2(ConstantSalinity): + pass + + +class Isatype3(SalinityInitialConditions): + pass + + +class Isatype4(SpatiallyVaryingTimeHistorySalinity): + pass diff --git a/pyschism/forcing/bctides/itetype.py b/pyschism/forcing/bctides/itetype.py index 227c8c73..dd51a8e6 100644 --- a/pyschism/forcing/bctides/itetype.py +++ b/pyschism/forcing/bctides/itetype.py @@ -86,8 +86,21 @@ def itetype(self): return 4 -Itetype0 = Itetype -Itetype1 = UniformTimeHistoryTemperature -Itetype2 = ConstantTemperature -Itetype3 = TemperatureInitialConditions -Itetype4 = SpatiallyVaryingTimeHistoryTemperature +class Itetype0(Itetype): + pass + + +class Itetype1(UniformTimeHistoryTemperature): + pass + + +class Itetype2(ConstantTemperature): + pass + + +class Itetype3(TemperatureInitialConditions): + pass + + +class Itetype4(SpatiallyVaryingTimeHistoryTemperature): + pass diff --git a/pyschism/forcing/bctides/itrtype.py b/pyschism/forcing/bctides/itrtype.py index 2307fa41..f15df641 100644 --- a/pyschism/forcing/bctides/itrtype.py +++ b/pyschism/forcing/bctides/itrtype.py @@ -53,7 +53,17 @@ def itrtype(self): return 4 -Itrtype1 = UniformTimeHistoryTracer -Itrtype2 = ConstantTracer -Itrtype3 = TracerInitialConditions -Itrtype4 = SpatiallyVaryingTimeHistoryTracer +class Itrtype1(UniformTimeHistoryTracer): + pass + + +class Itrtype2(ConstantTracer): + pass + + +class Itrtype3(TracerInitialConditions): + pass + + +class Itrtype4(SpatiallyVaryingTimeHistoryTracer): + pass From caee36a349a761445873ea2f324d2b23bc86fc2a Mon Sep 17 00:00:00 2001 From: cuill Date: Thu, 15 Dec 2022 11:37:41 -0500 Subject: [PATCH 07/31] Added param to specify ocean bounday --- pyschism/forcing/bctides/bctides.py | 8 ++++---- pyschism/forcing/hycom/hycom2schism.py | 25 +++++++++++++++---------- 2 files changed, 19 insertions(+), 14 deletions(-) diff --git a/pyschism/forcing/bctides/bctides.py b/pyschism/forcing/bctides/bctides.py index 49796578..50b0da6b 100644 --- a/pyschism/forcing/bctides/bctides.py +++ b/pyschism/forcing/bctides/bctides.py @@ -132,10 +132,10 @@ def write( start_date: datetime = None, end_date: Union[datetime, timedelta] = None, bctides: Union[bool, str] = True, - elev2D: Union[bool, str] = True, - uv3D: Union[bool, str] = True, - tem3D: Union[bool, str] = True, - sal3D: Union[bool, str] = True, + elev2D: Union[bool, str] = False, + uv3D: Union[bool, str] = False, + tem3D: Union[bool, str] = False, + sal3D: Union[bool, str] = False, overwrite: bool = False, parallel_download=False, progress_bar=True, diff --git a/pyschism/forcing/hycom/hycom2schism.py b/pyschism/forcing/hycom/hycom2schism.py index 1098ea22..d335f6a3 100644 --- a/pyschism/forcing/hycom/hycom2schism.py +++ b/pyschism/forcing/hycom/hycom2schism.py @@ -177,7 +177,7 @@ def __init__(self, hgrid, vgrid=None): self.hgrid = hgrid self.vgrid = Vgrid.default() if vgrid is None else vgrid - def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D=True, TS=True, UV=True, restart=False, adjust2D=False, lats=None, msl_shifts=None): + def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, ocean_bnd_ids = [0], elev2D=True, TS=True, UV=True, restart=False, adjust2D=False, lats=None, msl_shifts=None): outdir = pathlib.Path(outdir) self.start_date = start_date @@ -190,8 +190,10 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D= #Get open boundary gdf=self.hgrid.boundaries.open.copy() opbd=[] - for boundary in gdf.itertuples(): - opbd.extend(list(boundary.indexes)) + #for boundary in gdf.itertuples(): + # opbd.extend(list(boundary.indexes)) + for ibnd in ocean_bnd_ids: + opbd.extend(list(gdf.iloc[ibnd].indexes)) blon = self.hgrid.coords[opbd,0] blat = self.hgrid.coords[opbd,1] #logger.info(f'blon min {np.min(blon)}, max {np.max(blon)}') @@ -328,14 +330,15 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D= it = it0 + it1 database=get_database(date) - logger.info(f'Fetching data for {it}, {date} from database {database}') + logger.info(f'Fetching data for {date} from database {database}') #loop over each open boundary ind1 = 0 ind2 = 0 - for boundary in gdf.itertuples(): - - opbd = list(boundary.indexes) + #for boundary in gdf.itertuples(): + for ibnd in ocean_bnd_ids: + #opbd = list(boundary.indexes) + opbd = list(gdf.iloc[ibnd].indexes) ind1 = ind2 ind2 = ind1 + len(opbd) #logger.info(f'ind1 = {ind1}, ind2 = {ind2}') @@ -452,9 +455,10 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, elev2D= class Nudge: - def __init__(self): + def __init__(self, ocean_bnd_ids=[0]): self.include = None + self.ocean_bnd_ids = ocean_bnd_ids def gen_nudge(self, outdir: Union[str, os.PathLike], hgrid, rlmax = 1.5, rnu_day=0.25): @@ -491,8 +495,9 @@ def compute_nudge(lon, lat, nnode, opbd, out): #Get open boundary gdf=hgrid.boundaries.open.copy() opbd=[] - for boundary in gdf.itertuples(): - opbd.extend(list(boundary.indexes)) + #for boundary in gdf.itertuples(): + for ibnd in self.ocean_bnd_ids: + opbd.extend(list(gdf.iloc[ibnd].indexes)) opbd = np.array(opbd) elnode=hgrid.elements.array From d815c9ea6e03b4aab0614a986d6b6316385e02fb Mon Sep 17 00:00:00 2001 From: Jaime R Calzada Date: Thu, 15 Dec 2022 11:58:19 -0500 Subject: [PATCH 08/31] update tides.py --- pyschism/forcing/bctides/tides.py | 85 ++++++++++++++++++++++++++++++- 1 file changed, 83 insertions(+), 2 deletions(-) diff --git a/pyschism/forcing/bctides/tides.py b/pyschism/forcing/bctides/tides.py index 02a6aed7..9744ed16 100644 --- a/pyschism/forcing/bctides/tides.py +++ b/pyschism/forcing/bctides/tides.py @@ -47,6 +47,19 @@ class Tides: major_constituents = MAJOR_CONSTITUENTS minor_constituents = MINOR_CONSTITUENTS constituents = ALL_CONSTITUENTS + _nodal_factors = {} + _earth_equilibrium_arguments = {} + _tidal_species_types = {} + _amplitude_constants = {} + _amplitudes = { + 'elevation': {}, + 'velocity': {}, + } + + _phases = { + 'elevation': {}, + 'velocity': {}, + } def __init__( self, @@ -88,6 +101,9 @@ def __init__( for constituent in constituents: self.use_constituent(constituent) + + + def __iter__(self): for constituent in self.active_constituents: yield (constituent, self.get_tidal_constituent(constituent)) @@ -109,14 +125,43 @@ def __call__(self, start_date: datetime, rnday: Union[int, datetime], self.get_greenwich_factor(start_date, rnday, constituent)) # FACE* # noqa:E501 def get_elevation(self, constituent, vertices): + if constituent.lower() == 'z0': return np.full((vertices.shape[0],), float(self._Z0)), \ np.full((vertices.shape[0],), 0.) + + if constituent in self._amplitudes['elevation']: + amps = self._amplitudes['elevation'][constituent] + phases = self._phases['elevation'][constituent] + if isinstance(amps, (float, int)): + amps = np.full((vertices.shape[0],), float(amps)) + if isinstance(phases, (float, int)): + phases = np.full((vertices.shape[0],), float(phases)) + return amps, phases + + # if vertices is None: + # raise ValueError( + # 'Argument vertices must not be None for tidal constituent ' + # f'{constituent}.') return self.tidal_database.get_elevation(constituent, vertices) def get_velocity(self, constituent, vertices): if constituent.lower() == 'z0': return tuple(4*[np.full((vertices.shape[0],), 0.)]) + + if constituent in self._amplitudes['velocity']: + amps = self._amplitudes['velocity'][constituent] + phases = self._phases['velocity'][constituent] + if isinstance(amps, (float, int)): + amps = np.full((vertices.shape[0],), float(amps)) + if isinstance(phases, (float, int)): + phases = np.full((vertices.shape[0],), float(phases)) + return amps, phases + + # if vertices is None: + # raise ValueError( + # 'Argument vertices must not be None for tidal constituent ' + # f'{constituent}.') return self.tidal_database.get_velocity(constituent, vertices) def use_all(self, potential=True, forcing=True): @@ -148,6 +193,35 @@ def drop_constituent(self, constituent): raise ValueError("Argument constituent must be one of " f"{self.active_constituents}") self._active_constituents.pop(constituent) + + def add_constituent( + self, + name: str, + angular_frequency, + nodal_factor = 1., + earth_equilibrium_argument: float = 0., + amplitude_constant: float = None, + tidal_species_type: int = None, + elevation_amplitude: dict = None, + elevation_phase: dict = 0., + velocity_amplitude: dict = None, + velocity_phase = 0., + ): + # self._amplitudes['elevation'].update({name: amplitude}) + self.orbital_frequencies.update({name: angular_frequency}) + # self. + self._nodal_factors.update({name: nodal_factor}) + self._earth_equilibrium_arguments.update({name: earth_equilibrium_argument}) + self._tidal_species_types.update({name: tidal_species_type}) + self._amplitude_constants.update({name: amplitude_constant}) + self._amplitudes['elevation'].update({name: elevation_amplitude}) + self._amplitudes['velocity'].update({name: velocity_amplitude}) + self._phases['elevation'].update({name: elevation_phase}) + self._phases['velocity'].update({name: velocity_phase}) + self._active_constituents[name] = { + "potential": False, + "forcing": True, + } def add_Z0(self, value): self._active_constituents['Z0'] = { @@ -176,6 +250,7 @@ def get_tidal_species_type(self, constituent): return self.tidal_species_type[constituent] def get_orbital_frequency(self, constituent): + # if constituent in self. return self.orbital_frequencies[constituent] def get_initial_conditions(self, constituent, vertices): @@ -191,8 +266,10 @@ def set_Z0(self, Z0): def _manage_dates(f: Callable): def decorator(self, start_date, rnday, constituent): val = f(self, start_date, rnday, constituent) - del(self.start_date_utc) - del(self.end_date_utc) + if hasattr(self, 'start_date_utc'): + del(self.start_date_utc) + if hasattr(self, 'end_date_utc'): + del(self.end_date_utc) return val return decorator @@ -200,6 +277,8 @@ def decorator(self, start_date, rnday, constituent): def get_nodal_factor(self, start_date: datetime, rnday: Union[float, timedelta], constituent: str): + if constituent not in ALL_CONSTITUENTS: + return self._nodal_factors[constituent] if start_date.tzinfo is not None and \ start_date.tzinfo.utcoffset(start_date) is not None: self.start_date_utc = start_date.astimezone(timezone(timedelta(0))) @@ -300,6 +379,8 @@ def decorator(self, start_date, rnday, constituent): def get_greenwich_factor(self, start_date: datetime, rnday: Union[float, timedelta], constituent: str): + if constituent in self._earth_equilibrium_arguments: + return self._earth_equilibrium_arguments[constituent] if start_date.tzinfo is not None and \ start_date.tzinfo.utcoffset(start_date) is not None: self.start_date_utc = start_date.astimezone(timezone(timedelta(0))) From 6bf85de204d313f25cd3ad22d5916ff9c5ca0d27 Mon Sep 17 00:00:00 2001 From: Jaime R Calzada Date: Thu, 15 Dec 2022 13:28:28 -0500 Subject: [PATCH 09/31] multiboundary bctides 2D gen test --- pyschism/dates.py | 34 +++++---- pyschism/forcing/bctides/bctides.py | 107 +++++++++++++--------------- pyschism/forcing/bctides/bctypes.py | 7 +- pyschism/forcing/bctides/iettype.py | 5 ++ pyschism/forcing/bctides/ifltype.py | 5 ++ pyschism/forcing/bctides/isatype.py | 6 +- pyschism/forcing/bctides/itetype.py | 4 ++ pyschism/forcing/bctides/itrtype.py | 4 ++ pyschism/forcing/bctides/tides.py | 5 +- 9 files changed, 102 insertions(+), 75 deletions(-) diff --git a/pyschism/dates.py b/pyschism/dates.py index 7f18bc87..9f30ecea 100644 --- a/pyschism/dates.py +++ b/pyschism/dates.py @@ -11,14 +11,14 @@ def __init__(self): self.start_date = None def __set__(self, obj, start_date: datetime): - start_date = nearest_cycle() if start_date is None else \ - localize_datetime(start_date).astimezone(pytz.utc) - - if obj.end_date is not None: - if start_date > obj.end_date: - raise ValueError( - 'start_date is greater than end_date. ' - 'Try clearing end_date first.') + if start_date is not None: + start_date = localize_datetime(start_date).astimezone(pytz.utc) + + if obj.end_date is not None: + if start_date > obj.end_date: + raise ValueError( + 'start_date is greater than end_date. ' + 'Try clearing end_date first.') self.start_date = start_date def __get__(self, obj, val) -> datetime: @@ -37,9 +37,13 @@ def __set__(self, obj, run_days: Union[timedelta, float]): if isinstance(run_days, (int, float)): run_days = timedelta(days=float(run_days)) + + elif isinstance(run_days, datetime): + run_days = localize_datetime(run_days) - obj.start_date + elif not isinstance(run_days, timedelta): raise TypeError( - f'Argument run_days must be {timedelta} or float.') + f'Argument run_days must be {timedelta} or float, not {run_days}.') if obj.start_date is not None: obj.end_date = obj.start_date + run_days @@ -60,15 +64,21 @@ def __init__(self): def __set__(self, obj, val: Union[float, timedelta, datetime]): + if isinstance(val, (int, float)): + val = timedelta(days=float(val)) + if isinstance(val, datetime): val = localize_datetime(val).astimezone(pytz.utc) - elif not isinstance(val, timedelta): + elif not isinstance(val, timedelta) and obj.start_date is not None: val = obj.start_date + timedelta(days=float(val)) - elif isinstance(val, timedelta): + elif isinstance(val, timedelta) and obj.start_date is not None: val = obj.start_date + val + else: + val = None + self.end_date = val def __get__(self, obj, val) -> datetime: @@ -109,7 +119,7 @@ def nearest_zulu(input_datetime=None, method='floor'): def nearest_cycle(input_datetime=None, period=6, method='floor'): - + # Bug: method='ceil' doesn´t work correctly for cycle 18 if input_datetime is None: input_datetime = localize_datetime(datetime.utcnow()) assert method in ['floor', 'ceil'] diff --git a/pyschism/forcing/bctides/bctides.py b/pyschism/forcing/bctides/bctides.py index 8ee4f027..c17f222a 100644 --- a/pyschism/forcing/bctides/bctides.py +++ b/pyschism/forcing/bctides/bctides.py @@ -61,7 +61,7 @@ class Bctides(metaclass=BctidesMeta): start_date = dates.StartDate() end_date = dates.EndDate() - rnday = dates.RunDays() + run_days = dates.RunDays() def __init__( self, @@ -140,7 +140,7 @@ def write( if start_date is not None: self.start_date = start_date if end_date is not None: - self.end_date = end_date + self.run_days = end_date # self.tidal_database.write(path, ) output_directory = pathlib.Path(output_directory) output_directory.mkdir(exist_ok=overwrite, parents=True) @@ -154,57 +154,18 @@ def write( for boundary in self.gdf.itertuples(): data_source = getattr(boundary, bctype) if data_source is not None: + import importlib if hasattr(data_source, 'rlmax'): - # I admit this exec is hacky. - # pros: works well, it's simple, we don't need a return value - # cons: might be confusing to read. # This generates all the nudges and writes the nudge files. - exec( - f"from pyschism.forcing.bctides.nudge import {tracer}_Nudge;" - f"_tracer = output_directory / f'{tracer}_nudge.gr3' if {tracer.lower()}3D is True else {tracer};" - f"_tr={tracer}_Nudge(self, data_source, rlmax=data_source.rlmax, rnu_day=data_source.rnu_day);" - f'logger.info(f"Writing {tracer} nudge to file ' - + r'{_tracer}");' - "_tr.write(_tracer, overwrite=overwrite)" - ) - break - ## write nudge - #for bctype, tracer in {"itetype": "TEM", "isatype": "SAL"}.items(): - # for boundary in self.gdf.itertuples(): - # data_source = getattr(boundary, bctype) - # if data_source is not None: - - # # I admit this exec is hacky. - # # pros: works well, it's simple, we don't need a return value - # # cons: might be confusing to read. - # # This generates all the nudges and writes the nudge files. - # exec( - # f"from pyschism.forcing.bctides.nudge import {tracer}_Nudge;" - # f"_tracer = output_directory / f'{tracer}_nudge.gr3' if {tracer.lower()}3D is True else {tracer};" - # f"_tr={tracer}_Nudge(self, data_source, rlmax=data_source.rlmax, rnu_day=data_source.rnu_day);" - # f'logger.info(f"Writing {tracer} nudge to file ' - # + r'{_tracer}");' - # "_tr.write(_tracer, overwrite=overwrite)" - # ) - # break - # write nudge - for bctype, tracer in {"itetype": "TEM", "isatype": "SAL"}.items(): - for boundary in self.gdf.itertuples(): - data_source = getattr(boundary, bctype) - if data_source is not None: - if hasattr(data_source, 'rlmax'): - # I admit this exec is hacky. - # pros: works well, it's simple, we don't need a return value - # cons: might be confusing to read. - # This generates all the nudges and writes the nudge files. - exec( - f"from pyschism.forcing.bctides.nudge import {tracer}_Nudge;" - f"_tracer = output_directory / f'{tracer}_nudge.gr3' if {tracer.lower()}3D is True else {tracer};" - f"_tr={tracer}_Nudge(self, data_source, rlmax=data_source.rlmax, rnu_day=data_source.rnu_day);" - f'logger.info(f"Writing {tracer} nudge to file ' - + r'{_tracer}");' - "_tr.write(_tracer, overwrite=overwrite)" - ) + nudgemod = importlib.import_module('pyschism.forcing.bctides.nudge') + nudgeclass = getattr(nudgemod, f'{tracer}_Nudge') + _tracerfile = locals()[f'{tracer.lower()}3D'] + if _tracerfile is False: + continue + elif _tracerfile is True: + _tracerfile = output_directory / f'{tracer}_nudge.gr3' + nudgeclass(self, data_source, rlmax=data_source.rlmax, rnu_day=data_source.rnu_day + ).write(_tracerfile, overwrite=overwrite) break def write_elev2D(): @@ -292,17 +253,13 @@ def get_forcing_string(self, boundary, global_constituents): def get_focing_digit(bctype): if bctype is not None: - # sensitive to MRO. - return str( - getattr(bctype, f"{bctype.__class__.__bases__[0].__name__.lower()}") - ) + return bctype.forcing_digit return "0" line = [ f"{len(boundary.indexes)}", - *[digit for digit in map(get_focing_digit, bctypes)], + *[str(digit) for digit in map(get_focing_digit, bctypes)], ] - f = [" ".join(line)] for bctype in bctypes: if bctype is not None: @@ -313,6 +270,10 @@ def get_focing_digit(bctype): ) return "\n".join(f) + @property + def rnday(self): + return self.run_days + @cached_property def gdf(self): gdf = self.hgrid.boundaries.open.copy() @@ -396,3 +357,35 @@ def constituents(self): ) ) ) + + +def test(): + from datetime import datetime + import logging + + from pyschism.mesh import Hgrid + from pyschism.forcing.bctides import Bctides, iettype, ifltype, isatype, itetype + + # setup logging + logging.basicConfig( + format="[%(asctime)s] %(name)s %(levelname)s: %(message)s", + force=True, + ) + logging.getLogger("pyschism").setLevel(logging.DEBUG) + + startdate = datetime(2018, 8, 17) + print(startdate) + rnday = 61 + hgrid = Hgrid.open("./hgrid.gr3", crs="epsg:4326") + + # Bctides + iet3 = iettype.Iettype3(constituents='major', database='tpxo') + iet4 = iettype.Iettype4() + iet5 = iettype.Iettype5(iettype3=iet3, iettype4=iet4) + ifl3 = ifltype.Ifltype3(constituents='major', database='tpxo') + ifl4 = ifltype.Ifltype4() + ifl5 = ifltype.Ifltype5(ifltype3=ifl3, ifltype4=ifl4) + isa3 = isatype.Isatype4() + ite3 = itetype.Itetype4() + bctides = Bctides(hgrid, iettype={'1': iet5}, ifltype={'1': ifl5}, isatype=isa3, itetype=ite3) + bctides.write('./', startdate, rnday, bctides=True, elev2D=False, uv3D=True, tem3D=True, sal3D=False, overwrite=True) diff --git a/pyschism/forcing/bctides/bctypes.py b/pyschism/forcing/bctides/bctypes.py index e8b86a6b..c9ecc601 100644 --- a/pyschism/forcing/bctides/bctypes.py +++ b/pyschism/forcing/bctides/bctypes.py @@ -1,4 +1,4 @@ -from abc import abstractmethod +from abc import abstractmethod, abstractproperty from pyschism.forcing.base import ModelForcing @@ -9,5 +9,10 @@ class BoundaryForcing(ModelForcing): def get_boundary_string(self, boundary) -> str: pass + @abstractproperty + def forcing_digit(self): + """ bctype int value """ + raise NotImplementedError + Bctype = BoundaryForcing diff --git a/pyschism/forcing/bctides/iettype.py b/pyschism/forcing/bctides/iettype.py index 86f3a85b..42598ff9 100644 --- a/pyschism/forcing/bctides/iettype.py +++ b/pyschism/forcing/bctides/iettype.py @@ -8,11 +8,16 @@ class Iettype(Bctype): + @property @abstractmethod def iettype(self) -> int: pass + @property + def forcing_digit(self): + return self.iettype + class UniformTimeHistoryElevation(Iettype): def __init__(self, time_history: List[List[float]], boundaries): diff --git a/pyschism/forcing/bctides/ifltype.py b/pyschism/forcing/bctides/ifltype.py index 46c5583a..cf708039 100644 --- a/pyschism/forcing/bctides/ifltype.py +++ b/pyschism/forcing/bctides/ifltype.py @@ -7,11 +7,16 @@ class Ifltype(Bctype): + @property @abstractmethod def ifltype(self) -> int: """Returns integer representig SCHISM ifltype code for bctides.in""" + @property + def forcing_digit(self): + return self.ifltype + class UniformTimeHistoryVelocity(Ifltype): def __init__(self, time_history): diff --git a/pyschism/forcing/bctides/isatype.py b/pyschism/forcing/bctides/isatype.py index 9f5cb54f..3ddfcb6d 100644 --- a/pyschism/forcing/bctides/isatype.py +++ b/pyschism/forcing/bctides/isatype.py @@ -1,17 +1,21 @@ from abc import abstractmethod from enum import Enum -from matplotlib.pyplot import new_figure_manager from pyschism.forcing.bctides.bctypes import Bctype from pyschism.forcing import hycom class Isatype(Bctype): + @property @abstractmethod def isatype(self) -> int: pass + @property + def forcing_digit(self): + return self.isatype + class UniformTimeHistorySalinity(Isatype): def __init__(self, time_history): diff --git a/pyschism/forcing/bctides/itetype.py b/pyschism/forcing/bctides/itetype.py index dd51a8e6..293f4815 100644 --- a/pyschism/forcing/bctides/itetype.py +++ b/pyschism/forcing/bctides/itetype.py @@ -16,6 +16,10 @@ class Itetype(Bctype): def itetype(self) -> int: pass + @property + def forcing_digit(self): + return self.itetype + class UniformTimeHistoryTemperature(Itetype): diff --git a/pyschism/forcing/bctides/itrtype.py b/pyschism/forcing/bctides/itrtype.py index f15df641..f8fbdaf2 100644 --- a/pyschism/forcing/bctides/itrtype.py +++ b/pyschism/forcing/bctides/itrtype.py @@ -7,6 +7,10 @@ class Itrtype(Bctype): def itrtype(self) -> int: '''Returns integer representig SCHISM itrtype code for bctides.in''' + @property + def forcing_digit(self): + return self.itrtype + def get_boundary_string(self, hgrid, boundary): return '' diff --git a/pyschism/forcing/bctides/tides.py b/pyschism/forcing/bctides/tides.py index 9744ed16..65503133 100644 --- a/pyschism/forcing/bctides/tides.py +++ b/pyschism/forcing/bctides/tides.py @@ -55,7 +55,7 @@ class Tides: 'elevation': {}, 'velocity': {}, } - + _phases = { 'elevation': {}, 'velocity': {}, @@ -101,9 +101,6 @@ def __init__( for constituent in constituents: self.use_constituent(constituent) - - - def __iter__(self): for constituent in self.active_constituents: yield (constituent, self.get_tidal_constituent(constituent)) From a1c048eeef585d5632abade67c8b8c0fa3202e26 Mon Sep 17 00:00:00 2001 From: Jaime R Calzada Date: Thu, 15 Dec 2022 14:49:51 -0500 Subject: [PATCH 10/31] add implementation for constant temp forcing (itetype2) --- pyschism/forcing/bctides/bctides.py | 27 ++++++++++++++++++++++----- pyschism/forcing/bctides/itetype.py | 18 +++++++++++++++--- 2 files changed, 37 insertions(+), 8 deletions(-) diff --git a/pyschism/forcing/bctides/bctides.py b/pyschism/forcing/bctides/bctides.py index c17f222a..87ed634c 100644 --- a/pyschism/forcing/bctides/bctides.py +++ b/pyschism/forcing/bctides/bctides.py @@ -295,7 +295,7 @@ def nbfr(self): @cached_property def tides(self): return TidalConstituentCombiner(self.gdf) - + class TidalConstituentCombiner(Tides): def __init__(self, gdf): @@ -359,7 +359,7 @@ def constituents(self): ) -def test(): +def ad_hoc_test(): from datetime import datetime import logging @@ -386,6 +386,23 @@ def test(): ifl4 = ifltype.Ifltype4() ifl5 = ifltype.Ifltype5(ifltype3=ifl3, ifltype4=ifl4) isa3 = isatype.Isatype4() - ite3 = itetype.Itetype4() - bctides = Bctides(hgrid, iettype={'1': iet5}, ifltype={'1': ifl5}, isatype=isa3, itetype=ite3) - bctides.write('./', startdate, rnday, bctides=True, elev2D=False, uv3D=True, tem3D=True, sal3D=False, overwrite=True) + # ite3 = itetype.Itetype4() + bctides = Bctides(hgrid, iettype={'1': iet5}, ifltype={'1': ifl5}, + isatype=isa3, + itetype={'1': itetype.Itetype2(10, 1)} + ) + bctides.write( + './', + startdate, + rnday, + bctides=True, + elev2D=False, + uv3D=False, + tem3D=False, + sal3D=False, + overwrite=True + ) + + +if __name__ == "__main__": + ad_hoc_test() diff --git a/pyschism/forcing/bctides/itetype.py b/pyschism/forcing/bctides/itetype.py index 293f4815..c8add4df 100644 --- a/pyschism/forcing/bctides/itetype.py +++ b/pyschism/forcing/bctides/itetype.py @@ -34,14 +34,26 @@ def itetype(self) -> int: class ConstantTemperature(Itetype): - # def __init__(self, value, tobc: float = 1.): - # self.value = value - # super().__init__(tobc) + def __init__(self, value: float, nudging_factor: float): + self.value = value + if not (nudging_factor >= 0) and (nudging_factor <= 1): + raise ValueError( + 'Argument `nudging_factor` must be >= 0 and <= 1,' + f'but got {nudging_factor}.') + self.nudging_factor = nudging_factor + super().__init__() @property def itetype(self) -> int: return 2 + def get_boundary_string(self, *args, **kwargs) -> str: + boundary_string = [ + f'{self.value:0.6f}', + f'{self.nudging_factor:0.6f}', + ] + return '\n'.join(boundary_string) + class TemperatureInitialConditions(Itetype): From d29fa43bbf8798ac5cbee9442d9e233e6de2440b Mon Sep 17 00:00:00 2001 From: Jaime Calzada Date: Mon, 19 Dec 2022 12:22:56 -0500 Subject: [PATCH 11/31] read param defaults from static file --- pyschism/param/core.py | 5 +- pyschism/param/opt.py | 6 +- pyschism/param/param.nml | 1104 +++++++++++++++++++++++++++++++++ pyschism/param/param.py | 15 +- pyschism/param/schism_init.py | 12 +- 5 files changed, 1126 insertions(+), 16 deletions(-) create mode 100644 pyschism/param/param.nml diff --git a/pyschism/param/core.py b/pyschism/param/core.py index 1d97276e..1c682d5e 100644 --- a/pyschism/param/core.py +++ b/pyschism/param/core.py @@ -7,7 +7,6 @@ import f90nml from pyschism.enums import Stratification -from pyschism.param.schism_init import GitParamTemplate # class IbcType(Enum): @@ -202,7 +201,7 @@ def template(self): @template.setter def template(self, template: Union[str, os.PathLike, None]): if template is True: - template = GitParamTemplate().path + template = pathlib.Path(__file__) / 'param.nml' elif template is False or template is None: template = None else: @@ -212,5 +211,5 @@ def template(self, template: Union[str, os.PathLike, None]): @property def defaults(self): if self.template is None: - return GitParamTemplate().core + return f90nml.read(pathlib.Path(__file__).parent / 'param.nml')["core"] return f90nml.read(self.template)['core'] diff --git a/pyschism/param/opt.py b/pyschism/param/opt.py index cb2e4d6e..b6bcf365 100644 --- a/pyschism/param/opt.py +++ b/pyschism/param/opt.py @@ -1,7 +1,7 @@ from datetime import datetime, timedelta import logging -# import pathlib +import pathlib from typing import Union import f90nml @@ -10,9 +10,7 @@ from pyschism import dates from pyschism.mesh.fgrid import NchiType -from pyschism.param.schism_init import GitParamTemplate -PARAM_DEFAULTS = f90nml.read(GitParamTemplate().path)["opt"] logger = logging.getLogger(__name__) @@ -42,6 +40,8 @@ class OptMeta(type): def __new__(meta, name, bases, attrs): + + PARAM_DEFAULTS = f90nml.read(pathlib.Path(__file__).parent / 'param.nml')["opt"] for key, value in PARAM_DEFAULTS.items(): if key not in attrs: if isinstance(value, list): diff --git a/pyschism/param/param.nml b/pyschism/param/param.nml new file mode 100644 index 00000000..b9453459 --- /dev/null +++ b/pyschism/param/param.nml @@ -0,0 +1,1104 @@ +!parameter inputs via namelist convention. +!(1) Use ' ' (single quotes) for chars; +!(2) integer values are fine for real vars/arrays; +!(3) if multiple entries for a parameter are found, the last one wins - please avoid this +!(4) array inputs follow column major (like FORTRAN) and can spill to multiple lines. +! Values can be separated by commas or spaces. +!(5) space allowed before/after '=' + +&CORE +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Core (mandatory) parameters; no defaults +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Pre-processing option. Useful for checking grid errors etc. Need to use 1 +! core only for compute (plus necessary scribe cores). Under scribe I/O, the +! code (the scribe part) will hang but the outputs will be there. Just kill +! the job. + ipre = 0 !Pre-processor flag (1: on; 0: off) + +! Baroclinic/barotropic option. If ibc=0 (baroclinic model), ibtp is not used. + ibc = 0 !Baroclinic option + ibtp = 1 + + rnday = 30 !total run time in days + dt = 100. !Time step in sec + +! Grid for WWM (USE_WWM) + msc2 = 24 !same as msc in .nml ... for consitency check between SCHISM and WWM + mdc2 = 30 !same as mdc in .nml + +! Define # of tracers in tracer modules (if enabled) + ntracer_gen = 2 !user defined module (USE_GEN) + ntracer_age = 4 !age calculation (USE_AGE). Must be =2*N where N is # of age tracers + sed_class = 5 !SED3D (USE_SED) + eco_class = 27 !EcoSim (USE_ECO): must be between [25,60] + +! Global output controls + nspool = 36 !output step spool + ihfskip = 864 !stack spool; every ihfskip steps will be put into 1_*, 2_*, etc... +/ + +&OPT +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Optional parameters. The values shown below are default unless otherwise noted +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +!----------------------------------------------------------------------- +! Optional offline partitioning using NO_PARMETIS. If on, need partition.prop (like global_to_local.prop). +!----------------------------------------------------------------------- + + +!----------------------------------------------------------------------- +! 2nd pre-proc flag. If ipre2/=0, code will output some diagnotic outputs and stop. +! These outputs are: drag coefficients (Cdp) +!----------------------------------------------------------------------- + ipre2 = 0 + +!----------------------------------------------------------------------- +! Option to only solve tracer transport (and bypass most b-tropic solver) +! Usage: turn the flag on ('1' or '2') and turn on your tracer modules and make sure +! (1) your inputs are consistent with the original hydro-only run (with additional tracers +! of course); (2) hydro-only run results are in hydro_out/schout*.nc, which must +! have 'hvel_side' (not normal hvel), 'elev', 'diffusivity', 'temp_elem', 'salt_elem' (for +! new scribe outputs, use corresponding files/var names); (3) dt above is +! multiple of _output_ step used in the original hydro-only run +! (as found in in hydro_out/schout*.nc); e.g. dt = 1 hour. (4). When itransport_only=2, +! additional variables ('totalSuspendedLoad','sedBedStress') are needed. +! Hotstart should work also, but you'd probably not use an aggressively large dt especially +! when air-sea exchange is involved. +!----------------------------------------------------------------------- + itransport_only = 0 + +!----------------------------------------------------------------------- +! Option to add self-attracting and loading tide (SAL) into tidal potential +! (usually for basin-scale applications). +! If iloadtide=0, no SAL. +! If iloadtide=1, needs inputs: loadtide_[FREQ].gr3, +! where [FREQ] are freq names (shared with tidal potential, in upper cases) +! and the _two_ 'depths' inside are amplitude (m) and phases (degrees behind GMT), +! interpolated from global tide model (e.g. FES2014). In this option, SAL is +! lumped into tidal potential so it shares some parameters with tidal potential +! in bctides.in (cut-off depth, frequencies). +! If iloadtide=2 or 3, use a simple scaling for gravity approach (in this option, +! SAL is applied everywhere and does not share parameters with tidal potential). +! For iloadtide=2, a const scaling (1-loadtide_coef) is used; for iloadtide=3, the scaling is +! dependent on depth (Stepanov & Hughes 2004) with max of loadtide_coef. +!----------------------------------------------------------------------- + iloadtide = 0 + loadtide_coef = 0.1 !only used if iloadtide=2,3 (for '3', the default should be 0.12) + +! Starting time + start_year = 2000 !int + start_month = 1 !int + start_day = 1 !int + start_hour = 0 !double + utc_start = 8 !double + +!----------------------------------------------------------------------- +! Coordinate option: 1: Cartesian; 2: lon/lat (hgrid.gr3=hgrid.ll in this case, +! and orientation of element is outward of earth) +! Notes for lon/lat: make sure hgrid.ll and grid in sflux are consistent in +! longitude range! +!----------------------------------------------------------------------- + ics = 1 !Coordinate option + +!----------------------------------------------------------------------- +! Hotstart option. 0: cold start; 1: hotstart with time reset to 0; 2: +! continue from the step in hotstart.nc +!----------------------------------------------------------------------- + ihot = 0 + +!----------------------------------------------------------------------- +! Equation of State type used +! ieos_type=0: UNESCO 1980 (nonlinear); =1: linear function of T ONLY, i.e. +! \rho=eos_b+eos_a*T, where eos_a<=0 in kg/m^3/C +!----------------------------------------------------------------------- + ieos_type = 0 + ieos_pres = 0 !used only if ieos_type=0. 0: without pressure effects + eos_a = -0.1 !needed if ieos_type=1; should be <=0 + eos_b = 1001. !needed if ieos_type=1 + +!----------------------------------------------------------------------- +! Main ramp option +!----------------------------------------------------------------------- +! nramp = 1 !ramp-up option (1: on; 0: off) + dramp = 1. !ramp-up period in days for b.c. etc (no ramp-up if <=0) + +! nrampbc = 0 !ramp-up flag for baroclinic force + drampbc = 0. !ramp-up period in days for baroclinic force + +!----------------------------------------------------------------------- +! Method for momentum advection. 0: ELM; 1: upwind (not quite working yet) +!----------------------------------------------------------------------- + iupwind_mom = 0 + +!----------------------------------------------------------------------- +! Methods for computing velocity at nodes. +! If indvel=0, conformal linear shape function is used; if indvel=1, averaging method is used. +! For indvel=0, a stabilization method is needed (see below). +!----------------------------------------------------------------------- + indvel = 0 + +!----------------------------------------------------------------------- +! 2 stabilization methods, mostly for indvel=0. +! (1) Horizontal viscosity option. ihorcon=0: no viscosity is used; =1: Lapacian; +! =2: bi-harmonic. If ihorcon=1, horizontal viscosity _coefficient_ (<=1/8, related +! to diffusion number) is given in hvis_coef0, and the diffusion # +! is problem dependent; [0.001-1/8] seems to work well. +! If ihorcon=2, diffusion number is given by hvis_coef0 (<=0.025). +! If indvel=1, no horizontal viscosity is needed. +! (2) Shapiro filter (see below) +! +! For non-eddying regime applications (nearshore, estuary, river), an easiest option is: +! indvel=0, ishapiro=1 (shapiro0=0.5), ihorcon=inter_mom=0. +! For applications that include eddying regime, refer to the manual. +!----------------------------------------------------------------------- + ihorcon = 0 + hvis_coef0 = 0.025 !const. diffusion # if ihorcon/=0; <=0.025 for ihorcon=2, <=0.125 for ihorcon=1 +! cdh = 0.01 !needed only if ihorcon/=0; land friction coefficient - not active yet + +!----------------------------------------------------------------------- +! 2nd stabilization method via Shapiro filter. This should normally be used +! if indvel=0. To transition between eddying/non-eddying regimes, use +! indvel=0, ihorcon/=0, and ishapiro=-1 (requiring shapiro.gr3) or 2 (Smagorinsky-like filter). +!----------------------------------------------------------------------- + ishapiro = 1 !options + niter_shap = 1 !needed if ishapiro/=0: # of iterations with Shapiro filter + !shapiro0: Shapiro filter strength, needed only if ishapiro=1 or 2 + !If ishapiro=1, shapiro0 is the filter strength (max is 0.5). + !If ishapiro=2, shapiro0 is the coefficient in tanh(). Experiences so far suggest ~1.e3 + !If ishapiro=-1, this is not used as the filter strength is read in from shapiro.gr3 + shapiro0 = 0.5 + +!----------------------------------------------------------------------- +! Implicitness factor (0.5 Stokes drift advection (xy), Coriolis + fwvor_advz_stokes = 1 ! --> Stokes drift advection (z) , Coriolis + fwvor_gradpress = 1 ! --> Pressure term + fwvor_breaking = 1 ! --> Wave breaking + fwvor_streaming = 1 ! --> Wave streaming (works with iwbl /= 0) + fwvor_wveg = 0 ! --> Wave dissipation by vegetation acceleration term + fwvor_wveg_NL = 0 ! --> Non linear intrawave vegetation force (see Dean and Bender, 2006 or van Rooijen et al., 2016 for details) + cur_wwm = 0 ! Coupling current in WWM + ! 0: surface layer current + ! 1: depth-averaged current + ! 2: computed according to Kirby and Chen (1989) + wafo_obcramp = 0 ! Ramp on wave forces at open boundary (1: on / 0: off) + ! --> this option requires the input file wafo_ramp.gr3 + ! which defines the ramp value (between 0 and 1) + ! at nodes over the whole domain + +!----------------------------------------------------------------------- +! Bed deformation option (0: off; 1: vertical deformation only; 2: 3D bed deformation). +! If imm=1, bdef.gr3 is needed; if imm=2, user needs to update depth info etc +! in the code (not working for ics=2 yet). +!----------------------------------------------------------------------- + imm = 0 + ibdef = 10 !needed if imm=1; # of steps used in deformation + +!----------------------------------------------------------------------- +! Reference latitude for beta-plane approximation when ncor=1 (not used if ics=2) +!----------------------------------------------------------------------- + slam0 = -124 !lon - not really used + sfea0 = 45 !lat + +!----------------------------------------------------------------------- +! Option to deal with under resolution near steep slopes in deeper depths +! 0: use h[12,_bcc below; /=0: use hw_* below +!----------------------------------------------------------------------- + iunder_deep = 0 + +!----------------------------------------------------------------------- +! Baroclinicity calculation in off/nearshore with iunder_deep=ibc=0. +! The 'below-bottom' gradient +! is zeroed out if h>=h2_bcc (i.e. like Z) or uses const extrap +! (i.e. like terrain-following) if h<=h1_bcc(h1_bcc + +!----------------------------------------------------------------------- +! Hannah-Wright-like ratio & depth used to account for under-resolution +! in a b-clinic model. Used only if ibc=0 and iunder_deep/=0. +! The b-clinic force at prism centers is calculated with a reconstruction +! method in horizontal that has a stencil of an element and its adjacent elements. +! If the depths change is too much between the elem and its adjacent elem +! the under-resolution occurs (with steep bottom slope) and b-clinic force +! needs to be zeroed out below the (higher) bottom, specifically, if +! max(2 depths)>=hw_depth and abs(diff(2 depths))>=hw_ratio*max(2 depths). +!----------------------------------------------------------------------- + hw_depth = 1.e6 !threshold depth in [m] + hw_ratio = 0.5 !ratio + +!----------------------------------------------------------------------- +! Hydraulic model option. If ihydraulics/=0, hydraulics.in +! is required. This option cannot be used with non-hydrostatic model. +!----------------------------------------------------------------------- + ihydraulics = 0 + +!----------------------------------------------------------------------- +! Point sources/sinks option (0: no; 1: ASCII inputs; -1: netcdf). +! If =1, needs source_sink.in (list of elements), +! vsource,th, vsink.th, and msource.th (the source/sink values must be single precision). +! If =-1, all info is in source.nc +! and each type of volume/mass source/sink can have its own time step and +! # of records. +! Source/sinks can be specified at an elem more +! than once, and the code will accumulate the volumes; for mass conc, +! values are applied at _net_ source elem (no summation for conc). + +! meth_sink: options to treat sinks @ dry elem: no speacial treatment if meth_sink=0; +! zero out sink if the elem is dry if meth_sink=1 +!----------------------------------------------------------------------- + if_source = 0 +! nramp_ss = 1 !needed if if_source=1; ramp-up flag for source/sinks + dramp_ss = 2 !needed if if_source/=0; ramp-up period in days for source/sinks (no ramp-up if <=0) + meth_sink = 1 !options to treat sinks @ dry elem + +! Specify vertical level to inject source concentration for each tracer _module_. +! Code will extrapolate below bottom/above surface unless level 0 is specified. In +! the latter case, the incoming tracer mass will be distributed across all vertical layers - +! this approach works best in shallow waters. + lev_tr_source(1) = -9 !T + lev_tr_source(2) = -9 !S + lev_tr_source(3) = -9 !GEN + lev_tr_source(4) = -9 !AGE + lev_tr_source(5) = -9 !SED3D + lev_tr_source(6) = -9 !EcoSim + lev_tr_source(7) = -9 !ICM + lev_tr_source(8) = -9 !CoSINE + lev_tr_source(9) = -9 !Feco + lev_tr_source(10) = -9 !TIMOR + lev_tr_source(11) = -9 !FABM + lev_tr_source(12) = -9 !DVD + +!----------------------------------------------------------------------- +! Horizontal diffusivity option. if ihdif=1, horizontal diffusivity is given in hdif.gr3 +!----------------------------------------------------------------------- + ihdif = 0 + +!----------------------------------------------------------------------- +! Bottom friction. +! nchi=0: drag coefficients specified in drag.gr3; nchi=-1: Manning's +! formulation (even for 3D prisms) with n specified in manning.gr3. +! nchi=1: bottom roughness (in meters) specified in rough.gr3 (and in this case, negative +! or 0 depths in rough.gr3 indicate time-independent Cd, not roughness!). +! Cd is calculated using the log law, when dzb>=dzb_min; when dzb=h_tvd and the flag in +! tvd.prop = 1 for the elem; otherwise upwind is used for efficiency. +! itr_met=3 (horizontal TVD) or 4 (horizontal WENO): implicit TVD in the vertical dimension. +! Also if itr_met==3 and h_tvd>=1.e5, some parts of the code are bypassed for efficiency. +! Controls for WENO are not yet in place. +!----------------------------------------------------------------------- + itr_met = 3 + h_tvd = 5. !cut-off depth (m) + !If itr_met=3 or 4, need the following 2 tolerances of convergence. The convergence + !is achieved when sqrt[\sum_i(T_i^s+1-T_i^s)^2]<=eps1_tvd_imp*sqrt[\sum_i(T_i^s)^2]+eps2_tvd_imp + eps1_tvd_imp = 1.e-4 !suggested value is 1.e-4, but for large suspended load, need to use a smaller value (e.g. 1.e-9) + eps2_tvd_imp = 1.e-14 + + !Optional hybridized ELM transport for efficiency + ielm_transport = 0 !1: turn on + max_subcyc = 10 !used only if ielm_transport/=0. Max # of subcycling per time step in transport allowed + + !if itr_met = 4, the following parameters are needed + !if itr_met=4 and ipre=1, diagnostic outputs are generated for weno accuracy and stencil quality, + ! see subroutine weno_diag in src/Hydro/misc_subs.F90 for details + ip_weno = 2 !order of accuracy: 0- upwind; 1- linear polynomial, 2nd order; 2- quadratic polynomial, 3rd order + courant_weno=0.5 !Courant number for weno transport + nquad = 2 !number of quad points on each side, nquad= 1 or 2 + ntd_weno = 1 !order of temporal discretization: (1) Euler (default); (3): 3rd-order Runge-Kutta (only for benchmarking) + epsilon1 = 1.e-15 !coefficient for 2nd order weno smoother (larger values are more prone to numerical dispersion) + epsilon2 = 1.e-10 !1st coefficient for 3rd order weno smoother (larger values are more prone to numerical dispersion + !, 1.e-10 should be fairly safe, recommended values: 1.e-8 ~ 1.e-6 (from the real applications so far) + i_prtnftl_weno = 0 !option for writing nonfatal errors on invalid temp. or salinity for density: (0) off; (1) on. + + !Inactive at the moment: + epsilon3 = 1.e-25 !2nd coefficient for 3rd order weno smoother (inactive at the moment) + !Elad filter has not been implemented yet; preliminary tests showed it might not be necessary + ielad_weno = 0 !ielad, if ielad=1, use ELAD method to suppress dispersion (inactive at the moment) + small_elad = 1.e-4 !small (inactive at the moment) + +!----------------------------------------------------------------------- +! Atmos. option. nws=3 is reserved for coupling with atmospheric model. +! If nws=0, no atmos. forcing is applied. If nws=1, atmos. +! variables are read in from wind.th. If nws=2, atmos. variables are +! read in from sflux_ files. +! If nws=4, ascii format is used for wind and atmos. pressure at each node (see source code). +! If nws=-1 (requires USE_PAHM), use Holland parametric wind model (barotropic only with wind and atmos. pressure). +! In this case, the Holland model is called every step so wtiminc is not used. An extra +! input is needed: hurricane-track.dat +! +! Stress calculation: +! If nws=2, ihconsv=1 and iwind_form=0, the stress is calculated from heat exchange +! routine; +! Otherwise if iwind_form=-1, the stress is calculated from Pond & Pichard formulation; +! if iwind_form=1, Hwang (2018) formulation (Cd tapers off at high wind). +! If WWM is enabled and icou_elfe_wwm>0 and iwind_form=-2 or -3, stress is overwritten by WWM: +! If iwind_form=-2, stress=rho_air*ufric^2; scaled by rho_water +! If iwind_form=-3, the stress is calculated according to the param. of Donelan et al. (1993) based on the wave age. +! In all cases, if USE_ICE the stress in ice-covered portion is calculated by ICE routine. +!----------------------------------------------------------------------- + nws = 0 + wtiminc = 150. !time step for atmos. forcing. Default: same as dt +! nrampwind = 1 !ramp-up option for atmos. forcing + drampwind = 1. !ramp-up period in days for wind (no ramp-up if <=0) + iwindoff = 0 !needed only if nws/=0; '1': needs windfactor.gr3 + iwind_form = 1 !needed if nws/=0 + !If IMPOSE_NET_FLUX is on and nws=2, read in net _surface_ heat flux as var 'dlwrf' + !(Downward Long Wave) in sflux_rad (solar radiation is still used separately), + !and if PREC_EVAP is on, also read in net P-E as 'prate' (Surface Precipitation Rate) in sflux_prc. + +!----------------------------------------------------------------------- +! Heat and salt exchange using Zeng's (default) or Fairall (USE_BULK_FAIRALL) scheme. +! isconsv=1 needs ihconsv=1; ihconsv=1 needs nws=2. +! If isconsv=1, need to compile with precip/evap module turned on. +! If at least one of ihconsv and isconsv is set to 1, +! locally turning off air-sea exchange in shallow waters (<0.25 m) is recommended. +! Options for locally turning off heat/salt exchange are specified by +! i_hmin_airsea_ex, i_hmin_salt_ex respectively: +! 0: not turned off locally; +! 1: locally turned off, based on grid depth, i.e., +! off when the local grid depth (z in hgrid) < hmin_airsea_ex or hmin_salt_ex; +! 2 (recommended): locally turned off, based on local total water depth, i.e., +! off when the local total water depth < hmin_airsea_ex or hmin_salt_ex, +! hmin_airsea_ex, hmin_salt_ex: +! 0.2 m is recommended for both "1" and "2" +!----------------------------------------------------------------------- + ihconsv = 0 !heat exchange option + isconsv = 0 !evaporation/precipitation model + i_hmin_airsea_ex = 2 ! no effect if ihconsv=0 + hmin_airsea_ex = 0.2 ![m], no effect if ihconsv=0 + i_hmin_salt_ex = 2 ! no effect if isconsv=0 + hmin_salt_ex = 0.2 ![m], no effect if isconsv=0 + iprecip_off_bnd = 0 !if /=0, precip will be turned off near land bnd + +!----------------------------------------------------------------------- +! Turbulence closure. +!----------------------------------------------------------------------- + itur = 3 !Default: 0 + dfv0 = 1.e-2 !needed if itur=0 + dfh0 = 1.e-4 !needed if itur=0 + mid = 'KL' !needed if itur=3,5. Use KE if itur=5 + stab = 'KC' !needed if itur=3 or 5. Use 'GA' if turb_met='MY'; otherwise use 'KC'. + xlsc0 = 0.1 !needed if itur=3 or 5. Scale for surface & bottom mixing length (>0) + +!----------------------------------------------------------------------- +! Sponge layer for elevation and vel. +! If inu_elev=0, no relaxation is applied to elev. +! If inu_elev=1, relax. constants (in 1/sec, i.e. relax*dt<=1) are specified in elev_nudge.gr3 +! and applied to eta=0 (thus a depth=0 means no relaxation). +! Similarly for inu_uv (with input uv_nudge.gr3) +!----------------------------------------------------------------------- + inu_elev = 0 + inu_uv = 0 + +!----------------------------------------------------------------------- +! Nudging options for tracers. If inu_[MOD]=0, no nudging is used. If inu_tr=1, +! nudge to initial condition according to relaxation constants specified. +! If inu_tr=2, nudge to values in [MOD]_nu.nc (with step 'step_nu_tr'). +! The relaxation constants = [horizontal relax (specified in [MOD]_nudge.gr3) + or x +! vertical relax] times dt, where vertical relax is a linear function of +! vnh[1,2] and vnf[1,2], and [MOD] are tracer model names. 'nu_sum_mult' decides +! '+' or 'x' in the calculation of final relax. +! Code will ignore junk values (<=-99) inside [MOD]_nu.nc, so 1 way to avoid +! nudging for a tracer is to set its nudged values to -9999 +!----------------------------------------------------------------------- + inu_tr(1) = 0 !T + inu_tr(2) = 0 !S + inu_tr(3) = 0 !GEN + inu_tr(4) = 0 !Age + inu_tr(5) = 0 !SED3D + inu_tr(6) = 0 !EcoSim + inu_tr(7) = 0 !ICM + inu_tr(8) = 0 !CoSINE + inu_tr(9) = 0 !FIB + inu_tr(10) = 0 !TIMOR + inu_tr(11) = 0 !FABM + inu_tr(12) = 0 !DVD (must=0) + + nu_sum_mult=1 !1: final relax is sum of horizontal&vertical; 2: product + vnh1 = 400 !vertical nudging depth 1 + vnf1 = 0. !vertical relax \in [0,1] + vnh2 = 500 !vertical nudging depth 2 (must >vnh1) + vnf2 = 0. !vertical relax + + step_nu_tr = 86400. !time step [sec] in all [MOD]_nu.nc (for inu_[MOD]=2) + +!----------------------------------------------------------------------- +! Cut-off depth for cubic spline interpolation near bottom when computing horizontal gradients +! e.g. using hgrad_nodes() (radiation stress, and gradients of qnon and qhat in non-hydro model). +! If depth > h_bcc1 ('deep'), +! a min. (e.g. max bottom z-cor for the element) is imposed in the spline and so a more +! conservative method is used without extrapolation beyond bottom; +! otherwise constant extrapolation below bottom is used. +!----------------------------------------------------------------------- + h_bcc1 = 100. !h_bcc1 + +!----------------------------------------------------------------------- +! Dimensioning parameters for inter-subdomain btrack. +! If error occurs like 'bktrk_subs: overflow' or 'MAIN: nbtrk > mxnbt' +! gradually increasing these will solve the problem +!----------------------------------------------------------------------- + s1_mxnbt = 0.5 + s2_mxnbt = 3.5 + +!----------------------------------------------------------------------- +! Flag for harmonic analysis for elevation. If used , need to turn on USE_HA +! in Makefile, and input harm.in. Otherwise set it to 0. Hotstart ihot=2 is not working with HA +! Outputs are harme_* and use combine_outHA to combine. +!----------------------------------------------------------------------- + iharind = 0 + +!----------------------------------------------------------------------- +! Conservation check option. If iflux=1, some fluxes are computed +! in regions specified in fluxflag.prop (regional number from -1 to an arbitrary integer). +! in output flux.out, positive means flux from region n to region n-1 (n>=1) +! output file flux.dat need to be saved before continuing hotstart +!----------------------------------------------------------------------- + iflux = 0 + iflux_out_format = 0 !0: basic flux output; 1: more detailed outputs + +!----------------------------------------------------------------------- +! Test flags for debugging. These flags should be turned off normally. +!----------------------------------------------------------------------- +! Williamson test #5 (zonal flow over an isolated mount); if +! on, ics must =2 +!----------------------------------------------------------------------- + izonal5 = 0 !"0" - no test; otherwise on + +!----------------------------------------------------------------------- +! Rotating Gausshill test with stratified T,S (1: on; 0: off) +! Surface T,S read in from *.ic; code generates stratification +!----------------------------------------------------------------------- + ibtrack_test = 0 + +!----------------------------------------------------------------------- +! Rouse profile test (1: on; 0: off) +! If on, must turn on USE_TIMOR +!----------------------------------------------------------------------- + irouse_test = 0 + +!----------------------------------------------------------------------- +! Flag to choose FIB model for bacteria decay (used with USE_FIB) +! flag_fib = 1 - Constant decay rate (/day) in .gr3 format +! (kkfib_1.gr3 and kkfib_2.gr3) +! flag_fib = 2 - Decay rate computed from Canteras et al., 1995 +! flag_fib = 3 - Decay rate computed from Servais et al., 2007 +!---------------------------------------------------------------------- + flag_fib = 1 + +!---------------------------------------------------------------------- +! Marsh model parameters (only if USE_MARSH is on) +!---------------------------------------------------------------------- + slr_rate = 120. !sea-level rise rate in mm/yr, times morphological acceleration if used + +!---------------------------------------------------------------------- +! Vegetation model +! If isav=1, need 4 extra inputs: (1) sav_D.gr3 (depth is stem diameter in meters); +! (2) sav_N.gr3 (depth is # of stems per m^2); +! (3) sav_h.gr3 (height of canopy in meters); +! (4) sav_cd.gr3 (drag coefficient). +! If one of these depths=0 at a node, the code will set all to 0. +! If USE_MARSH is on and isav=1, all .gr3 must have constant depths! +!---------------------------------------------------------------------- + isav = 0 !on/off flag + +!---------------------------------------------------------------------- +! Coupling step with ICE module. +!---------------------------------------------------------------------- + nstep_ice = 1 !call ice module every nstep_ice steps of SCHISM + + +!---------------------------------------------------------------------- +! Specify level #'s if age module is invoked (USE_AGE), for 1st half of tracers only +!---------------------------------------------------------------------- + level_age = 9, -999 !default: -999 (all levels) + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! Physical constants +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +!----------------------------------------------------------------------- +! Earth's radii at pole and equator (to define an ellipsoid) +!----------------------------------------------------------------------- + rearth_pole = 6378206.4 + rearth_eq = 6378206.4 + +!----------------------------------------------------------------------- +! Specific heat of water (C_p) in J/kg/K +!----------------------------------------------------------------------- + shw = 4184.d0 + +!----------------------------------------------------------------------- +! Reference water density for Boussinesq approximation +!----------------------------------------------------------------------- + rho0 = 1000.d0 !kg/m^3 + +!----------------------------------------------------------------------- +! Fraction of vertical flux closure adjustment applied at surface, then subtracted +! from all vertical fluxes. This is currently done for T,S only +! 0.0 <= vclose_surf_frac < 1.0 +! 1: fully from surface (i.e. no correction as before); 0: fully from bottom +!----------------------------------------------------------------------- + vclose_surf_frac=1.0 + +!----------------------------------------------------------------------- +! Option to enforce strict mass conservation for each tracer model (only works with itr_met=3,4) +! At moment the scheme has not accounted for bottom 'leaking' (e.g. in SED), +! so iadjust_mass_consv0(5) must =0 +!----------------------------------------------------------------------- + iadjust_mass_consv0(1)=0 !T + iadjust_mass_consv0(2)=0 !S + iadjust_mass_consv0(3)=0 !GEN + iadjust_mass_consv0(4)=0 !AGE + iadjust_mass_consv0(5)=0 !SED3D (code won't allow non-0 for this module) + iadjust_mass_consv0(6)=0 !EcoSim + iadjust_mass_consv0(7)=0 !ICM + iadjust_mass_consv0(8)=0 !CoSiNE + iadjust_mass_consv0(9)=0 !Feco + iadjust_mass_consv0(10)=0 !TIMOR + iadjust_mass_consv0(11)=0 !FABM + iadjust_mass_consv0(12)=0 !DVD (must=0) + +! For ICM, impose mass conservation for depths larger than a threshold by considering prism +! volume change from step n to n+1. rinflation_icm is the max ratio btw H^{n+1} and H^n allowed. + h_massconsv = 2. ![m] + rinflation_icm = 1.e-3 + +/ + +&SCHOUT +!----------------------------------------------------------------------- +! Output section - all optional. Values shown are default unless otherwise stated, +! and default for most global outputs is off +!----------------------------------------------------------------------- + +!----------------------------------------------------------------------- +! Main switch to control netcdf. If =0, SCHISM won't output nc files +! at all (useful for other programs like ESMF to output) +!----------------------------------------------------------------------- + nc_out = 1 + +!----------------------------------------------------------------------- +! UGRID option for _3D_ outputs under scribed IO (out2d*.nc always has meta +! data info). If iof_ugrid/=0, 3D outputs will also have UGRID metadata (at +! the expense of file size). +!----------------------------------------------------------------------- + iof_ugrid = 0 + +!----------------------------------------------------------------------- +! Option for hotstart outputs +!----------------------------------------------------------------------- + nhot = 0 !1: output *_hotstart every 'hotout_write' steps + nhot_write = 8640 !must be a multiple of ihfskip if nhot=1 + +!----------------------------------------------------------------------- +! Station output option. If iout_sta/=0, need output skip (nspool_sta) and +! a station.in. If ics=2, the cordinates in station.in must be in lon., lat, +! and z (positive upward; not used for 2D variables). +!----------------------------------------------------------------------- + iout_sta = 0 + nspool_sta = 10 !needed if iout_sta/=0; mod(nhot_write,nspool_sta) must=0 + +!----------------------------------------------------------------------- +! Global output options +! The variable names that appear in nc output are shown in {} +!----------------------------------------------------------------------- + iof_hydro(1) = 1 !0: off; 1: on - elev. [m] {elevation} 2D + iof_hydro(2) = 0 !air pressure [Pa] {airPressure} 2D + iof_hydro(3) = 0 !air temperature [C] {airTemperature} 2D + iof_hydro(4) = 0 !Specific humidity [-] {specificHumidity} 2D + iof_hydro(5) = 0 !Net downward solar (shortwave) radiation after albedo [W/m/m] {solarRadiation} 2D + iof_hydro(6) = 0 !sensible flux (positive upward) [W/m/m] {sensibleHeat} 2D + iof_hydro(7) = 0 !latent heat flux (positive upward) [W/m/m] {latentHeat} 2D + iof_hydro(8) = 0 !upward longwave radiation (positive upward) [W/m/m] {upwardLongwave} 2D + iof_hydro(9) = 0 !downward longwave radiation (positive downward) [W/m/m] {downwardLongwave} 2D + iof_hydro(10) = 0 !total flux=-flsu-fllu-(radu-radd) [W/m/m] {totalHeat} 2D + iof_hydro(11) = 0 !evaporation rate [kg/m/m/s] {evaporationRate} 2D + iof_hydro(12) = 0 !precipitation rate [kg/m/m/s] {precipitationRate} 2D + iof_hydro(13) = 0 !Bottom stress vector [kg/m/s^2(Pa)] {bottomStressX,Y} 2D vector + iof_hydro(14) = 0 !wind velocity vector [m/s] {windSpeedX,Y} 2D vector + iof_hydro(15) = 0 !wind stress vector [m^2/s/s] {windStressX,Y} 2D vector + iof_hydro(16) = 0 !depth-averaged vel vector [m/s] {depthAverageVelX,Y} 2D vector + iof_hydro(17) = 0 !vertical velocity [m/s] {verticalVelocity} 3D + iof_hydro(18) = 0 !water temperature [C] {temperature} 3D + iof_hydro(19) = 0 !water salinity [PSU] {salinity} 3D + iof_hydro(20) = 0 !water density [kg/m^3] {waterDensity} 3D + iof_hydro(21) = 0 !vertical eddy diffusivity [m^2/s] {diffusivity} 3D + iof_hydro(22) = 0 !vertical eddy viscosity [m^2/s] {viscosity} 3D + iof_hydro(23) = 0 !turbulent kinetic energy {turbulentKineticEner} 3D + iof_hydro(24) = 0 !turbulent mixing length [m] {mixingLength} 3D +! iof_hydro(25) = 1 !z-coord {zCoordinates} 3D - this flag should be on for visIT etc + iof_hydro(26) = 1 !horizontal vel vector [m/s] {horizontalVelX,Y} 3D vector + iof_hydro(27) = 0 !horizontal vel vector defined @side [m/s] {horizontalSideVelX,Y} 3D vector + iof_hydro(28) = 0 !vertical vel. @elem [m/s] {verticalVelAtElement} 3D + iof_hydro(29) = 0 !T @prism centers [C] {temperatureAtElement} 3D + iof_hydro(30) = 0 !S @prism centers [PSU] {salinityAtElement} 3D + iof_hydro(31) = 0 !Barotropic pressure gradient force vector (m.s-2) @side centers {pressure_gradient} 2D vector + +!----------------------------------------------------------------------- +! Outputs from optional modules. Only uncomment these if the USE_* is on +!----------------------------------------------------------------------- +!----------------------------------------------------------------------- +! Outputs from DVD (USE_DVD must be on in Makefile) +!----------------------------------------------------------------------- +! iof_dvd(1) = 1 !num mixing for S (PSU^2/s) {DVD_1} 3D elem + +!----------------------------------------------------------------------- +! Outputs from WWM (USE_WWM must be on in Makefile) +!----------------------------------------------------------------------- +! iof_wwm(1) = 0 !sig. height (m) {sigWaveHeight} 2D +! iof_wwm(2) = 0 !Mean average period (sec) - TM01 {meanWavePeriod} 2D +! iof_wwm(3) = 0 !Zero down crossing period for comparison with buoy (s) - TM02 {zeroDowncrossPeriod} 2D +! iof_wwm(4) = 0 !Average period of wave runup/overtopping - TM10 {TM10} 2D +! iof_wwm(5) = 0 !Mean wave number (1/m) {meanWaveNumber} 2D +! iof_wwm(6) = 0 !Mean wave length (m) {meanWaveLength} 2D +! iof_wwm(7) = 0 !Mean average energy transport direction (degr) - MWD in NDBC? {meanWaveDirection} 2D +! iof_wwm(8) = 0 !Mean directional spreading (degr) {meanDirSpreading} 2D +! iof_wwm(9) = 0 !Discrete peak period (sec) - Tp {peakPeriod} 2D +! iof_wwm(10) = 0 !Continuous peak period based on higher order moments (sec) {continuousPeakPeriod} 2D +! iof_wwm(11) = 0 !Peak phase vel. (m/s) {peakPhaseVel} 2D +! iof_wwm(12) = 0 !Peak n-factor {peakNFactor} 2D +! iof_wwm(13) = 0 !Peak group vel. (m/s) {peakGroupVel} 2D +! iof_wwm(14) = 0 !Peak wave number {peakWaveNumber} 2D +! iof_wwm(15) = 0 !Peak wave length {peakWaveLength} 2D +! iof_wwm(16) = 0 !Peak (dominant) direction (degr) {dominantDirection} 2D +! iof_wwm(17) = 0 !Peak directional spreading {peakSpreading} 2D +! iof_wwm(18) = 0 !Discrete peak direction (radian?) {discretePeakDirectio} 2D +! iof_wwm(19) = 0 !Orbital vel. (m/s) {orbitalVelocity} 2D +! iof_wwm(20) = 0 !RMS Orbital vel. (m/s) {rmsOrbitalVelocity} 2D +! iof_wwm(21) = 0 !Bottom excursion period (sec?) {bottomExcursionPerio} 2D +! iof_wwm(22) = 0 !Bottom wave period (sec) {bottomWavePeriod} 2D +! iof_wwm(23) = 0 !Uresell number based on peak period {UresellNumber} 2D +! iof_wwm(24) = 0 !Friction velocity (m/s?) {frictionalVelocity} 2D +! iof_wwm(25) = 0 !Charnock coefficient {CharnockCoeff} 2D +! iof_wwm(26) = 0 !Rougness length {rougnessLength} 2D + +! iof_wwm(27) = 0 !Roller energy dissipation rate (W/m²) @nodes {Drol} 2D +! iof_wwm(28) = 0 !Total wave energy dissipation rate by depth-induced breaking (W/m²) @nodes {wave_sbrtot} 2D +! iof_wwm(29) = 0 !Total wave energy dissipation rate by bottom friction (W/m²) @nodes {wave_sbftot} 2D +! iof_wwm(30) = 0 !Total wave energy dissipation rate by whitecapping (W/m²) @nodes {wave_sdstot} 2D +! iof_wwm(31) = 0 !Total wave energy dissipation rate by vegetation (W/m²) @nodes {wave_svegtot} 2D +! iof_wwm(32) = 0 !Total wave energy input rate from atmospheric forcing (W/m²) @nodes {wave_sintot} 2D +! iof_wwm(33) = 0 !WWM_energy vector {waveEnergyDirX,Y} 2D vector + +! iof_wwm(34) = 0 !Vertical Stokes velocity (m.s-1) @sides and whole levels {stokes_wvel} 3D +! iof_wwm(35) = 0 !Wave force vector (m.s-2) computed by wwm @side centers and whole levels {waveForceX,Y} 3D vector + +! iof_wwm(36) = 0 !Horizontal Stokes velocity (m.s-1) @nodes and whole levels {stokes_hvel} 3D vector +! iof_wwm(37) = 0 !Roller contribution to horizontal Stokes velocity (m.s-1) @nodes and whole levels {roller_stokes_hvel} 3D vector + +!----------------------------------------------------------------------- +! Tracer module outputs. In most cases, actual # of outputs depends on # of tracers used +!----------------------------------------------------------------------- +! Outputs for user-defined tracer module (USE_GEN) +!----------------------------------------------------------------------- +! iof_gen(1) = 0 !1st tracer {GEN_1} 3D +! iof_gen(2) = 0 !2nd tracer {GEN_2} 3D + +!----------------------------------------------------------------------- +! Outputs for (age). Indices from "1" to "ntracer_age/2"; [days] +!----------------------------------------------------------------------- +! iof_age(1) = 0 !{AGE_1} 3D +! iof_age(2) = 0 !{AGE_2} 3D + +!----------------------------------------------------------------------- +! Specific outputs in SED3D (USE_SED must be on in Makefile; +! otherwise these are not needed) +!----------------------------------------------------------------------- +! iof_sed(1) = 0 ! total bed thickness @elem (m) {sedBedThickness} 2D +! iof_sed(2) = 0 ! total bed age over all layers @elem (sec) {sedBedAge} 2D +! iof_sed(3) = 0 ! Sediment transport roughness length @elem (m) (sedTransportRough) {z0st} 2D +! iof_sed(4) = 0 !current-ripples roughness length @elem (m) (sedRoughCurrentRippl) {z0cr} 2D +! iof_sed(5) = 0 !sand-waves roughness length (m) @elem (z0sw_elem) {sedRoughSandWave} 2D +! iof_sed(6) = 0 !wave-ripples roughness length @elem (m) (z0wr_elem) {sedRoughWaveRipple} 2D + +! iof_sed(7) = 0 !bottom depth _change_ from init. condition (m) {sedDepthChange} 2D +! iof_sed(8) = 0 ! Bed median grain size in the active layer (mm) {sedD50} 2D +! iof_sed(9) = 0 ! Bottom shear stress (Pa) {sedBedStress} 2D +! iof_sed(10) = 0 ! Bottom roughness lenghth (mm) {sedBedRoughness} 2D +! iof_sed(11) = 0 ! sediment porosity in the top layer (-) {sedPorocity} 2D +! iof_sed(12) = 0 ! erosion flux for suspended transport (kg/m/m/s) {sedErosionalFlux} 2D +! iof_sed(13) = 0 ! deposition flux for suspended transport (kg/m/m/s) {sedDepositionalFlux} 2D +! iof_sed(14) = 0 ! Bedload transport rate vector due to wave acceleration (kg/m/s) {sedBedloadTransportX,Y} 2D vector + +! Example of using 2 classes +! iof_sed(15) = 0 !Bedload transport rate _vector_ (kg.m-1.s-1) for 1st tracer (one output per class) {sedBedload[X,Y]_1} 2D vector +! iof_sed(16) = 0 !Bedload transport of 2nd class {sedBedFraction_[X,Y]_2} 2D vector +! iof_sed(17) = 0 !Bed fraction 1st tracer (one output per class) [-] {sedBedFraction_1} 2D +! iof_sed(18) = 0 !Bed fraction of 2nd class {sedBedFraction_2} 2D + +! iof_sed(19) = 0 !conc. of 1st class (one output need by each class) [g/L] {sedConcentration_1} 3D +! iof_sed(20) = 0 !conc. of 2nd class {sedConcentration_2} 3D + +! iof_sed(21) = 0 !total suspended concentration (g/L) {totalSuspendedLoad} 3D + +!----------------------------------------------------------------------- +! EcoSim outputs +!----------------------------------------------------------------------- +! iof_eco(1) = 0 {ECO_1} 3D + +!----------------------------------------------------------------------- +! ICM outputs +!----------------------------------------------------------------------- +! !core Module +! iof_icm_core(1) = 1 !PB1 +! iof_icm_core(2) = 1 !PB2 +! iof_icm_core(3) = 1 !PB3 +! iof_icm_core(4) = 1 !RPOC +! iof_icm_core(5) = 1 !LPOC +! iof_icm_core(6) = 1 !DOC +! iof_icm_core(7) = 1 !RPON +! iof_icm_core(8) = 1 !LPON +! iof_icm_core(9) = 1 !DON +! iof_icm_core(10) = 1 !NH4 +! iof_icm_core(11) = 1 !NO3 +! iof_icm_core(12) = 1 !RPOP +! iof_icm_core(13) = 1 !LPOP +! iof_icm_core(14) = 1 !DOP +! iof_icm_core(15) = 1 !PO4 +! iof_icm_core(16) = 1 !COD +! iof_icm_core(17) = 1 !DOX +! +! !silica module +! iof_icm_silica(1) = 1 !SU +! iof_icm_silica(2) = 1 !SA +! +! !zooplankton module +! iof_icm_zb(1) = 1 !ZB1 +! iof_icm_zb(2) = 1 !ZB2 +! +! !pH model +! iof_icm_ph(1) = 1 !TIC +! iof_icm_ph(2) = 1 !ALK +! iof_icm_ph(3) = 1 !CA +! iof_icm_ph(4) = 1 !CACO3 +! +! !CBP model +! iof_icm_cbp(1) = 1 !SRPOC +! iof_icm_cbp(2) = 1 !SRPON +! iof_icm_cbp(3) = 1 !SRPOP +! iof_icm_cbp(4) = 1 !PIP +! +! !SAV model +! iof_icm_sav(1) = 1 !stleaf: total leaf biomass @elem [gC/m^2] +! iof_icm_sav(2) = 1 !ststem: total stem biomass @elem [gC/m^2] +! iof_icm_sav(3) = 1 !stroot: total root biomass @elem [gC/m^2] +! iof_icm_sav(4) = 1 !sht: canopy height @elem [m] +! +! !VEG model +! iof_icm_veg(1) = 1 !vtleaf1: leaf biomass group 1 [gC/m^2] +! iof_icm_veg(2) = 1 !vtleaf2: leaf biomass group 2 [gC/m^2] +! iof_icm_veg(3) = 1 !vtleaf3: leaf biomass group 3 [gC/m^2] +! iof_icm_veg(4) = 1 !vtstem1: stem biomass group 1 [gC/m^2] +! iof_icm_veg(5) = 1 !vtstem2: stem biomass group 2 [gC/m^2] +! iof_icm_veg(6) = 1 !vtstem3: stem biomass group 3 [gC/m^2] +! iof_icm_veg(7) = 1 !vtroot1: root biomass group 1 [gC/m^2] +! iof_icm_veg(8) = 1 !vtroot2: root biomass group 2 [gC/m^2] +! iof_icm_veg(9) = 1 !vtroot3: root biomass group 3 [gC/m^2] +! iof_icm_veg(10) = 1 !vht1: canopy height group 1 [m] +! iof_icm_veg(11) = 1 !vht2: canopy height group 2 [m] +! iof_icm_veg(12) = 1 !vht3: canopy height group 3 [m] +! +! !SFM model +! iof_icm_sed(1) = 1 !bPOC1 (g.m-3) +! iof_icm_sed(2) = 1 !bPOC2 (g.m-3) +! iof_icm_sed(3) = 1 !bPOC3 (g.m-3) +! iof_icm_sed(4) = 1 !bPON1 (g.m-3) +! iof_icm_sed(5) = 1 !bPON2 (g.m-3) +! iof_icm_sed(6) = 1 !bPON3 (g.m-3) +! iof_icm_sed(7) = 1 !bPOP1 (g.m-3) +! iof_icm_sed(8) = 1 !bPOP2 (g.m-3) +! iof_icm_sed(9) = 1 !bPOP3 (g.m-3) +! iof_icm_sed(10) = 1 !bNH4 (g.m-3) +! iof_icm_sed(11) = 1 !bNO3 (g.m-3) +! iof_icm_sed(12) = 1 !bPO4 (g.m-3) +! iof_icm_sed(13) = 1 !bH2S (g.m-3) +! iof_icm_sed(14) = 1 !bCH4 (g.m-3) +! iof_icm_sed(15) = 1 !bPOS (g.m-3) +! iof_icm_sed(16) = 1 !bSA (g.m-3) +! iof_icm_sed(17) = 1 !bstc: surface transfer coeff. (m/day) +! iof_icm_sed(18) = 1 !bSTR: benthic stress (day) +! iof_icm_sed(19) = 1 !bThp: consective days of hypoxia (day) +! iof_icm_sed(20) = 1 !bTox: consective days of oxic condition after hypoxia (day) +! iof_icm_sed(21) = 1 !SOD (g.m-2.day-1) +! iof_icm_sed(22) = 1 !JNH4 (g.m-2.day-1) +! iof_icm_sed(23) = 1 !JNO3 (g.m-2.day-1) +! iof_icm_sed(24) = 1 !JPO4 (g.m-2.day-1) +! iof_icm_sed(25) = 1 !JSA (g.m-2.day-1) +! iof_icm_sed(26) = 1 !JCOD (g.m-2.day-1) +! +! !Benthic Algae model +! iof_icm_ba(1) = 1 !BA (g[C].m-2) +! +! !ICM Debug Outputs (need coding, for developers) +! iof_icm_dbg(1) = 1 !2D ICM debug variables +! iof_icm_dbg(2) = 1 !3D ICM debug variables + +!----------------------------------------------------------------------- +! CoSINE outputs: all 3D +!----------------------------------------------------------------------- +! iof_cos(1) = 0 !NO3 (uM) +! iof_cos(2) = 0 !SiO4(uM) +! iof_cos(3) = 0 !NH4 (uM) +! iof_cos(4) = 0 !S1 (uM) +! iof_cos(5) = 0 !S2 (uM) +! iof_cos(6) = 0 !Z1 (uM) +! iof_cos(7) = 0 !Z2 (uM) +! iof_cos(8) = 0 !DN (uM) +! iof_cos(9) = 0 !DSi (uM) +! iof_cos(10) = 0 !PO4 (uM) +! iof_cos(11) = 0 !DOX (uM) +! iof_cos(12) = 0 !CO2 (uM) +! iof_cos(13) = 0 !ALK (uM) + +!----------------------------------------------------------------------- +! Fecal indicating bacteria module +!----------------------------------------------------------------------- +! iof_fib(1) = 0 ! {FIB_1} 3D + +!----------------------------------------------------------------------- +! Specific outputs in SED2D (USE_SED2D must be on in Makefile; +! otherwise these are not needed) - not implemented yet +!----------------------------------------------------------------------- +! iof_sed2d(1) = 0 !bottom depth _change_ from init. condition (m) {SED2D_depth_change} +! iof_sed2d(2) = 0 !drag coefficient used in transport formulae SED2D_Cd{} +! iof_sed2d(3) = 0 !Courant number (b.qtot.dt / h.dx) {SED2D_cflsed} +! iof_sed2d(4) = 0 !Top layer d50 (m) {SED2D_d50} +! iof_sed2d(5) = 0 !total transport rate vector (kg/m/s) {SED2D_total_transport} +! iof_sed2d(6) = 0 !suspended tranport rate vector (kg/m/s) {SED2D_susp_load} +! iof_sed2d(7) = 0 !bedload transport rate vector (kg/m/s) {SED2D_bed_load} +! iof_sed2d(8) = 0 !time averaged total transport rate vector (kg/m/s) {SED2D_average_transport} +! iof_sed2d(9) = 0 !bottom slope vector (m/m); negative uphill {SED2D_bottom_slope} +! iof_sed2d(10) = 0 !Total roughness length @elem (m) (z0eq) {z0eq} +! iof_sed2d(11) = 0 !current-ripples roughness length @elem (m) (z0cr) {z0cr} +! iof_sed2d(12) = 0 !sand-waves roughness length @elem (m) (z0sw) {z0sw} +! iof_sed2d(13) = 0 !wave-ripples roughness length @elem (m) (z0wr) {z0wr} + +!----------------------------------------------------------------------- +! marsh flags (USE_MARSH on) +!----------------------------------------------------------------------- +! iof_marsh(1) = 0 ! {marshFlag} 2D elem + +!----------------------------------------------------------------------- +! Ice module outputs (if USE_ICE is on) +!----------------------------------------------------------------------- +! iof_ice(1)= 0 !strain rate @ elem [1/sec] {iceStrainRate} 2D +! iof_ice(2)= 0 !ice advective velcoity vector [m/s] {iceVelocityX,Y} 2D vector +! iof_ice(3)= 0 !net heat flux to ocean (>0 warm up SST) [W/m/m] {iceNetHeatFlux} 2D +! iof_ice(4)= 0 !net fresh water flux to ocean (>0 freshens up SSS) [kg/s/m/m] {iceFreshwaterFlux} 2D +! iof_ice(5)= 0 !ice temperature [C] at air-ice interface {iceTopTemperature} 2D + +! iof_ice(6)= 0 !ice volume [m] {iceTracer_1} 2D +! iof_ice(7)= 0 !ice concentration [-] {iceTracer_2} 2D +! iof_ice(8)= 0 !snow volume [m] {iceTracer_3} 2D + +!----------------------------------------------------------------------- +! Analysis module outputs (USE_ANALYSIS) +!----------------------------------------------------------------------- +! iof_ana(1) = 0 !min time step at each element over all subcycles in horizontal transport solver [s] {minTransportTimeStep} 2D +! iof_ana(2) = 0 !x-component of \nabla air_pres /\rho_0 [m/s/s] {airPressureGradientX} 2D +! iof_ana(3) = 0 !y-component of \nabla air_pres /\rho_0 [m/s/s] {airPressureGradientY} 2D +! iof_ana(4) = 0 !\alpha*g*\nabla \Psi [m/s/s] (gradient of tidal potential) {tidePotentialGradX} 2D +! iof_ana(5) = 0 !\alpha*g*\nabla \Psi [m/s/s] {tidePotentialGradY} 2D +! iof_ana(6) = 0 !\nabla \cdot (\mu \nabla u) [m/s/s] (horizontal momentum mixing) {horzontalViscosityX} 3D side +! iof_ana(7) = 0 !\nabla \cdot (\mu \nabla v) [m/s/s] {horzontalViscosityY} 3D side +! iof_ana(8) = 0 !-g/rho0* \int_z^\eta dr_dx dz [m/s/s] (b-clinic gradient) {baroclinicForceX} 3D side +! iof_ana(9) = 0 !-g/rho0* \int_z^\eta dr_dy dz [m/s/s] {baroclinicForceY} 3D side +! iof_ana(10) = 0 !d (\nu du/dz)/dz [m/s/s] - no vegetation effects (vertical momentum mixing) {verticalViscosityX} 3D side +! iof_ana(11) = 0 !d (\nu dv/dz)/dz [m/s/s] - no vegetation effects {verticalViscosityY} 3D side +! iof_ana(12) = 0 !(u \cdot \nabla) u [m/s/s] (momentum advection) {mommentumAdvectionX} 3D side +! iof_ana(13) = 0 !(u \cdot \nabla) u [m/s/s] {mommentumAdvectionY} 3D side +! iof_ana(14) = 0 !gradient Richardson number [-] {gradientRichardson} 3D +/ diff --git a/pyschism/param/param.py b/pyschism/param/param.py index 87873893..e56b3176 100644 --- a/pyschism/param/param.py +++ b/pyschism/param/param.py @@ -1,7 +1,8 @@ from datetime import timedelta +from typing import Union import logging import pathlib -from typing import Union +import tempfile import f90nml @@ -10,11 +11,9 @@ from pyschism.param.core import CORE from pyschism.param.opt import OPT from pyschism.param.schout import SCHOUT - +from pyschism.param import schism_init # from pyschism.stations import Stations -# PARAM_TEMPLATE = pathlib.Path(__file__).parent / 'param.nml.template' -from pyschism.param.schism_init import GitParamTemplate logger = logging.getLogger(__name__) @@ -48,8 +47,12 @@ def write(self, path, overwrite=False, use_template=False): if path.is_file() and not overwrite: raise IOError(f"File {path} exists and overwrite=False") if use_template: - PARAM_TEMPLATE = GitParamTemplate().path if use_template is True else use_template - f90nml.patch(PARAM_TEMPLATE, self.to_dict(), path) + PARAM_TEMPLATE = pathlib.Path(__file__).parent / 'param.nml' if use_template is True else use_template + schism_param_sample = schism_init.read_schism_param_sample_patched(PARAM_TEMPLATE) + tmpfile = tempfile.NamedTemporaryFile() + with open(tmpfile.name, 'w') as fh: + fh.write(schism_param_sample) + f90nml.patch(tmpfile.name, self.to_dict(), path) else: with open(path, "w") as f: f.write(str(self)) diff --git a/pyschism/param/schism_init.py b/pyschism/param/schism_init.py index fee5ea09..d8da0abb 100644 --- a/pyschism/param/schism_init.py +++ b/pyschism/param/schism_init.py @@ -111,7 +111,8 @@ def schout(self): return SchoutParser(self.schism_init_f90) -def patch_sample_param_to_avoid_f90nml_bug(schism_param_sample): +def read_schism_param_sample_patched(schism_param_file): + schism_param_sample = pathlib.Path(schism_param_file).read_text().decode('utf-8').split('\n') pattern = r"^[! _A-Za-z0-9]+\([*\d]*(\d+)[^\d]*\)[ ]*=" var_collection = [] for i, line in enumerate(schism_param_sample): @@ -132,7 +133,7 @@ def patch_sample_param_to_avoid_f90nml_bug(schism_param_sample): else: schism_param_sample.insert(i + 1, f"{var}({varcount}) = 0") break - + return '\n'.join(schism_param_sample) class Singleton(type): def __init__(cls, name, bases, dict): @@ -149,8 +150,11 @@ class GitParamTemplate(metaclass=Singleton): def __init__(self, branch="master"): url = f"https://raw.githubusercontent.com/schism-dev/schism/{branch}/sample_inputs/param.nml" response = urllib.request.urlopen(url) - schism_param_sample = response.read().decode("utf-8").split("\n") - patch_sample_param_to_avoid_f90nml_bug(schism_param_sample) + tmpfile = tempfile.NamedTemporaryFile() + with open(tmpfile.name) as fh: + fh.write(response.read().decode("utf-8")) + # schism_param_sample = response.read().decode("utf-8").split("\n") + schism_param_sample = read_schism_param_sample_patched(tmpfile.name) tmpfile1 = tempfile.NamedTemporaryFile() with open(tmpfile1.name, "w") as f: f.write("\n".join(schism_param_sample)) From c5d1a32819e478f0643e6ce497c61f12f0553883 Mon Sep 17 00:00:00 2001 From: cuill Date: Fri, 6 Jan 2023 11:18:14 -0500 Subject: [PATCH 12/31] Added parameters for nudge distance/strength --- pyschism/forcing/hycom/hycom2schism.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/pyschism/forcing/hycom/hycom2schism.py b/pyschism/forcing/hycom/hycom2schism.py index d335f6a3..3d9c0c2a 100644 --- a/pyschism/forcing/hycom/hycom2schism.py +++ b/pyschism/forcing/hycom/hycom2schism.py @@ -523,7 +523,7 @@ def compute_nudge(lon, lat, nnode, opbd, out): #logger.info(f'len of nudge idxs is {len(idxs)}') logger.info(f'It took {time() -t0} sencods to calcuate nudge coefficient') - nudge = [f"{rlmax}, {rnu_day}"] + nudge = [f"rlmax={rlmax}, rnu_day={rnu_day}"] nudge.extend("\n") nudge.append(f"{NE} {NP}") nudge.extend("\n") @@ -551,7 +551,7 @@ def compute_nudge(lon, lat, nnode, opbd, out): return self.include - def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, rnday, restart = False): + def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, rnday, restart = False, rlmax = None, rnu_day=None): outdir = pathlib.Path(outdir) @@ -565,8 +565,12 @@ def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, vd=Vgrid.open(vgrid) sigma=vd.sigma + #define nudge zone and strength + rlmax = 1.5 if rlmax is None else rlmax + rnu_day = 0.25 if rnu_day is None else rnu_day + logger.info(f'Max relax distance is {rlmax} degree, max relax strengh is {rnu_day} days.') #Get the index for nudge - include = self.gen_nudge(outdir,hgrid) + include = self.gen_nudge(outdir,hgrid, rlmax = rlmax, rnu_day=rnu_day) #get coords of SCHISM loni=hgrid.nodes.coords[:,0] From 43cd9e69dbdb5ea965f315ddcc39b8426219c29b Mon Sep 17 00:00:00 2001 From: cuill Date: Mon, 23 Jan 2023 13:49:16 -0500 Subject: [PATCH 13/31] fixed a bug for time_idx_restart --- pyschism/forcing/hycom/hycom2schism.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/pyschism/forcing/hycom/hycom2schism.py b/pyschism/forcing/hycom/hycom2schism.py index 3d9c0c2a..b5dc374a 100644 --- a/pyschism/forcing/hycom/hycom2schism.py +++ b/pyschism/forcing/hycom/hycom2schism.py @@ -248,7 +248,7 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, ocean_b dst_elev.createVariable('time_series', 'f', ('time', 'nOpenBndNodes', 'nLevels', 'nComponents')) #dst_elev['time_series'][:,:,:,:] = timeseries_el - elif restart: + elif elev2D and restart: dst_elev = Dataset(outdir / 'elev2D.th.nc', 'a', format='NETCDF4') time_idx_restart = dst_elev['time'][:].shape[0] @@ -289,9 +289,10 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, ocean_b dst_temp.createVariable('time_series', 'f', ('time', 'nOpenBndNodes', 'nLevels', 'nComponents')) #dst_temp['time_series'][:,:,:,:] = timeseries_t - elif restart: + elif TS and restart: dst_salt = Dataset(outdir / 'SAL_3D.th.nc', 'a', format='NETCDF4') dst_temp = Dataset(outdir / 'TEM_3D.th.nc', 'a', format='NETCDF4') + time_idx_restart = dst_salt['time'][:].shape[0] if UV and restart == False: #timeseries_uv=np.zeros([ntimes,NOP,nvrt,nComp2]) @@ -312,8 +313,9 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, ocean_b dst_uv.createVariable('time_series', 'f', ('time', 'nOpenBndNodes', 'nLevels', 'nComponents')) #dst_uv['time_series'][:,:,:,:] = timeseries_uv - elif restart: + elif UV and restart: dst_uv = Dataset(outdir / 'uv3D.th.nc', 'a', format='NETCDF4') + time_idx_restart = dst_uv['time'][:].shape[0] logger.info('**** Accessing GOFS data*****') t0=time() @@ -652,7 +654,7 @@ def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, logger.info('**** Accessing GOFS data*****') if restart: timevector = self.timevector[time_idx_restart+1:] - it0 = time_idx_restart + it0 = time_idx_restart+1 else: timevector = self.timevector it0 = 0 @@ -786,7 +788,7 @@ def fetch_data(self, date, fmt='schism'): f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' - foutname = f'ST_{date.strftime("%Y%m%d")}.nc' + foutname = f'TS_{date.strftime("%Y%m%d")}.nc' logger.info(f'filename is {foutname}') ds = xr.open_dataset(url_ts) From 0a05af91580d21fb7ea5f2179aef689e1a0a534e Mon Sep 17 00:00:00 2001 From: cuill Date: Fri, 3 Feb 2023 16:30:33 -0500 Subject: [PATCH 14/31] relove issues when boundary crosses the prime meridian --- pyschism/forcing/hycom/hycom2schism.py | 100 ++++++++++++++----------- 1 file changed, 57 insertions(+), 43 deletions(-) diff --git a/pyschism/forcing/hycom/hycom2schism.py b/pyschism/forcing/hycom/hycom2schism.py index b5dc374a..6662700b 100644 --- a/pyschism/forcing/hycom/hycom2schism.py +++ b/pyschism/forcing/hycom/hycom2schism.py @@ -58,10 +58,10 @@ def get_database(date, Bbox=None): logger.info('No data for {date}') return database -def get_idxs(date, database, bbox): +def get_idxs(date, database, bbox, blonc, blatc): - if date.strftime("%Y-%m-%d") >= datetime.now().strftime("%Y-%m-%d"): - date2 = datetime.now() - timedelta(days=1) + if date >= datetime.utcnow(): + date2 = datetime.utcnow() - timedelta(days=1) baseurl = f'https://tds.hycom.org/thredds/dodsC/{database}/FMRC/runs/GLBy0.08_930_FMRC_RUN_{date2.strftime("%Y-%m-%dT12:00:00Z")}?depth[0:1:-1],lat[0:1:-1],lon[0:1:-1],time[0:1:-1]' else: baseurl=f'https://tds.hycom.org/thredds/dodsC/{database}?lat[0:1:-1],lon[0:1:-1],time[0:1:-1],depth[0:1:-1]' @@ -73,8 +73,20 @@ def get_idxs(date, database, bbox): lon=ds['lon'][:] lat=ds['lat'][:] dep=ds['depth'][:] - lat_idxs=np.where((lat>=bbox.ymin-2.0)&(lat<=bbox.ymax+2.0))[0] - lon_idxs=np.where((lon>=bbox.xmin-2.0) & (lon<=bbox.xmax+2.0))[0] + + #check if hycom's lon is the same range as schism's + if not (bbox.xmin >= lon.min() and bbox.xmax <= lon.max()): + if lon.min() >= 0: + logger.info(f'Convert HYCOM longitude from [0, 360) to [-180, 180):') + idxs = lon>=180 + lon[idxs] = lon[idxs]-360 + elif lon.min() <= 0: + logger.info(f'Convert HYCOM longitude from [-180, 180) to [0, 360):') + idxs = lon<=0 + lon[idxs] = lon[idxs]+360 + + lat_idxs=np.where((lat>=bbox.ymin-0.5)&(lat<=bbox.ymax+0.5))[0] + lon_idxs=np.where((lon>=bbox.xmin-0.5) & (lon<=bbox.xmax+0.5))[0] lon=lon[lon_idxs] lat=lat[lat_idxs] #logger.info(lon_idxs) @@ -86,14 +98,11 @@ def get_idxs(date, database, bbox): lat_idx2=lat_idxs[-1].item() #logger.info(f'lat_idx1 is {lat_idx1}, lat_idx2 is {lat_idx2}') - for ilon in np.arange(len(lon)): - if lon[ilon] > 180: - lon[ilon] = lon[ilon]-360. #lonc=(np.max(lon)+np.min(lon))/2.0 - #logger.info(f'lonc is {lonc}') - #latc=(np.max(lat)+np.min(lat))/2.0 - #logger.info(f'latc is {latc}') - x2, y2=transform_ll_to_cpp(lon, lat) + #logger.info(f'lonc is {blonc}') + ##latc=(np.max(lat)+np.min(lat))/2.0 + #logger.info(f'latc is {blatc}') + x2, y2=transform_ll_to_cpp(lon, lat, blonc, blatc) idxs=np.where( date == times)[0] #check if time_idx is empty @@ -117,10 +126,6 @@ def get_idxs(date, database, bbox): return time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 def transform_ll_to_cpp(lon, lat, lonc=-77.07, latc=24.0): - #lonc=(np.max(lon)+np.min(lon))/2.0 - #logger.info(f'lonc is {lonc}') - #latc=(np.max(lat)+np.min(lat))/2.0 - #logger.info(f'latc is {latc}') longitude=lon/180*np.pi latitude=lat/180*np.pi radius=6378206.4 @@ -135,6 +140,12 @@ def interp_to_points_3d(dep, y2, x2, bxyz, val): idxs = np.where(abs(val) > 10000) val[idxs] = float('nan') + if not np.all(x2[:-1] <= x2[1:]): + logger.info('x2 is not in stricitly ascending order! Sorting x2 and val') + idxs = np.argsort(x2) + x2 = x2[idxs] + val = val[:, :, idxs] + val_fd = sp.interpolate.RegularGridInterpolator((dep,y2,x2),np.squeeze(val),'linear', bounds_error=False, fill_value = float('nan')) val_int = val_fd(bxyz) idxs = np.isnan(val_int) @@ -150,6 +161,12 @@ def interp_to_points_2d(y2, x2, bxy, val): idxs = np.where(abs(val) > 10000) val[idxs] = float('nan') + if not np.all(x2[:-1] <= x2[1:]): + logger.info('x2 is not in stricitly ascending order! Sorting x2 and val') + idxs = np.argsort(x2) + x2 = x2[idxs] + val = val[:, idxs] + val_fd = sp.interpolate.RegularGridInterpolator((y2,x2),np.squeeze(val),'linear', bounds_error=False, fill_value = float('nan')) val_int = val_fd(bxy) idxs = np.isnan(val_int) @@ -346,7 +363,10 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, ocean_b #logger.info(f'ind1 = {ind1}, ind2 = {ind2}') blon = self.hgrid.coords[opbd,0] blat = self.hgrid.coords[opbd,1] - xi,yi = transform_ll_to_cpp(blon, blat) + blonc = blon.mean() + blatc = blat.mean() + #logger.info(f'blonc = {blon.mean()}, blatc = {blat.mean()}') + xi,yi = transform_ll_to_cpp(blon, blat, blonc, blatc) bxy = np.c_[yi, xi] if TS or UV: @@ -361,21 +381,22 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, ocean_b xmin, xmax = np.min(blon), np.max(blon) ymin, ymax = np.min(blat), np.max(blat) + bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) - if date.strftime("%Y-%m-%d") >= datetime(2017, 2, 1).strftime("%Y-%m-%d") and \ - date.strftime("%Y-%m-%d") < datetime(2017, 6, 1).strftime("%Y-%m-%d") or \ - date.strftime("%Y-%m-%d") >= datetime(2017, 10, 1).strftime("%Y-%m-%d"): - xmin = xmin + 360. if xmin < 0 else xmin - xmax = xmax + 360. if xmax < 0 else xmax - bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) - else: - bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) + #if date >= datetime(2017, 2, 1) and \ + # date < datetime(2017, 6, 1) or \ + # date >= datetime(2017, 10, 1): + # xmin = xmin + 360. if xmin < 0 else xmin + # xmax = xmax + 360. if xmax < 0 else xmax + # bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) + #else: + # bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) #logger.info(f'xmin is {xmin}, xmax is {xmax}') - time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 = get_idxs(date, database, bbox) + time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 = get_idxs(date, database, bbox, blonc, blatc) - if date.strftime("%Y-%m-%d") >= datetime.now().strftime("%Y-%m-%d"): - date2 = datetime.now() - timedelta(days=1) + if date >= datetime.utcnow(): + date2 = datetime.utcnow() - timedelta(days=1) url = f'https://tds.hycom.org/thredds/dodsC/{database}/FMRC/runs/GLBy0.08_930_FMRC_RUN_' + \ f'{date2.strftime("%Y-%m-%dT12:00:00Z")}?depth[0:1:-1],lat[{lat_idx1}:1:{lat_idx2}],' + \ f'lon[{lon_idx1}:1:{lon_idx2}],time[{time_idx}],' + \ @@ -398,7 +419,7 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, ocean_b ds=Dataset(url) dep=ds['depth'][:] - logger.info('****Interpolation starts****') + logger.info(f'****Interpolation starts for boundary {ibnd}****') #ndt[it]=it*24*3600. @@ -589,7 +610,9 @@ def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, #Get open nudge array nlon = hgrid.coords[include, 0] nlat = hgrid.coords[include, 1] - xi,yi = transform_ll_to_cpp(nlon, nlat) + nlonc = nlon.mean() + nlatc = nlat.mean() + xi,yi = transform_ll_to_cpp(nlon, nlat, nlonc, nlatc) bxy = np.c_[yi, xi] zcor2=zcor[include,:] @@ -669,22 +692,13 @@ def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, xmin, xmax = np.min(nlon), np.max(nlon) ymin, ymax = np.min(nlat), np.max(nlat) - - if date.strftime("%Y-%m-%d") >= datetime(2017, 2, 1).strftime("%Y-%m-%d") and \ - date.strftime("%Y-%m-%d") < datetime(2017, 6, 1).strftime("%Y-%m-%d") or \ - date.strftime("%Y-%m-%d") >= datetime(2017, 10, 1).strftime("%Y-%m-%d"): - logger.info('Convert xmin and xmax') - xmin = xmin + 360. if xmin < 0 else xmin - xmax = xmax + 360. if xmax < 0 else xmax - bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) - else: - bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) + bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) #logger.info(f'xmin is {xmin}, xmax is {xmax}') - time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 = get_idxs(date, database, bbox) + time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 = get_idxs(date, database, bbox, nlonc, nlatc) - if date.strftime("%Y-%m-%d") >= datetime.now().strftime("%Y-%m-%d"): - date2 = datetime.now() - timedelta(days=1) + if date >= datetime.utcnow(): + date2 = datetime.utcnow() - timedelta(days=1) url = f'https://tds.hycom.org/thredds/dodsC/{database}/FMRC/runs/GLBy0.08_930_FMRC_RUN_' + \ f'{date2.strftime("%Y-%m-%dT12:00:00Z")}?depth[0:1:-1],lat[{lat_idx1}:1:{lat_idx2}],' + \ f'lon[{lon_idx1}:1:{lon_idx2}],time[{time_idx}],' + \ From 52dec166264328c5d6e48e23ad13bf866933043e Mon Sep 17 00:00:00 2001 From: cuill Date: Wed, 8 Feb 2023 10:46:02 -0500 Subject: [PATCH 15/31] added tmpdir option when size of the downloaded file exceeds sys tmp's limit --- pyschism/forcing/nws/nws2/era5.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/pyschism/forcing/nws/nws2/era5.py b/pyschism/forcing/nws/nws2/era5.py index 7bd77217..5da23da8 100644 --- a/pyschism/forcing/nws/nws2/era5.py +++ b/pyschism/forcing/nws/nws2/era5.py @@ -23,13 +23,15 @@ class ERA5DataInventory: - def __init__(self, start_date=None, rnday: Union[float, timedelta] = 4, bbox=None): + def __init__(self, start_date=None, rnday: Union[float, timedelta] = 4, bbox=None, tmpdir=None): self.start_date = start_date self.rnday = rnday self.end_date = self.start_date+timedelta(self.rnday + 1) self.client=cdsapi.Client() self._bbox = bbox + if tmpdir is not None: + self.tmpdir = tmpdir r = self.client.retrieve( 'reanalysis-era5-single-levels', @@ -60,6 +62,10 @@ def tmpdir(self): self._tmpdir = tempfile.TemporaryDirectory() return pathlib.Path(self._tmpdir.name) + @tmpdir.setter + def tmpdir(self, value): + self._tmpdir = value + @property def files(self): return sorted(list(self.tmpdir.glob('**/era5_*.nc'))) @@ -274,6 +280,7 @@ def write( bbox: Bbox = None, overwrite: bool=False, output_interval: int = 1, + tmpdir = None, ): self.start_date=start_date self.rnday=rnday @@ -294,6 +301,7 @@ def write( self.start_date, self.rnday, bbox, + tmpdir = tmpdir ) logger.info('Finished downloading ERA5') From 18f701f8e0e3dd8298aad78893e8154935b65a64 Mon Sep 17 00:00:00 2001 From: cuill Date: Wed, 8 Feb 2023 10:48:34 -0500 Subject: [PATCH 16/31] refined download hycom data --- pyschism/forcing/hycom/hycom2schism.py | 194 ++++++++++++------------- 1 file changed, 96 insertions(+), 98 deletions(-) diff --git a/pyschism/forcing/hycom/hycom2schism.py b/pyschism/forcing/hycom/hycom2schism.py index 6662700b..8408dff6 100644 --- a/pyschism/forcing/hycom/hycom2schism.py +++ b/pyschism/forcing/hycom/hycom2schism.py @@ -55,10 +55,10 @@ def get_database(date, Bbox=None): elif date >= datetime(1994, 1, 1) and date < datetime(2016, 1, 1): database = f'GLBv0.08/expt_53.X/data/{date.year}' else: - logger.info('No data for {date}') + raise ValueError(f'No data fro {date}!') return database -def get_idxs(date, database, bbox, blonc, blatc): +def get_idxs(date, database, bbox, lonc=None, latc=None): if date >= datetime.utcnow(): date2 = datetime.utcnow() - timedelta(days=1) @@ -98,11 +98,13 @@ def get_idxs(date, database, bbox, blonc, blatc): lat_idx2=lat_idxs[-1].item() #logger.info(f'lat_idx1 is {lat_idx1}, lat_idx2 is {lat_idx2}') - #lonc=(np.max(lon)+np.min(lon))/2.0 - #logger.info(f'lonc is {blonc}') - ##latc=(np.max(lat)+np.min(lat))/2.0 - #logger.info(f'latc is {blatc}') - x2, y2=transform_ll_to_cpp(lon, lat, blonc, blatc) + if lonc is None: + lonc = lon.mean() + logger.info(f'lonc is {lonc}') + if latc is None: + latc = lat.mean() + logger.info(f'latc is {latc}') + x2, y2=transform_ll_to_cpp(lon, lat, lonc, latc) idxs=np.where( date == times)[0] #check if time_idx is empty @@ -393,7 +395,7 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, ocean_b # bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) #logger.info(f'xmin is {xmin}, xmax is {xmax}') - time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 = get_idxs(date, database, bbox, blonc, blatc) + time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 = get_idxs(date, database, bbox, lonc=blonc, latc=blatc) if date >= datetime.utcnow(): date2 = datetime.utcnow() - timedelta(days=1) @@ -695,7 +697,7 @@ def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) #logger.info(f'xmin is {xmin}, xmax is {xmax}') - time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 = get_idxs(date, database, bbox, nlonc, nlatc) + time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 = get_idxs(date, database, bbox, lonc=nlonc, latc=nlatc) if date >= datetime.utcnow(): date2 = datetime.utcnow() - timedelta(days=1) @@ -747,104 +749,100 @@ def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, class DownloadHycom: - def __init__(self, hgrid): + def __init__(self, hgrid=None, bbox=None): - xmin, xmax = hgrid.coords[:, 0].min(), hgrid.coords[:, 0].max() - ymin, ymax = hgrid.coords[:, 1].min(), hgrid.coords[:, 1].max() - xmin = xmin + 360. if xmin < 0 else xmin - xmax = xmax + 360. if xmax < 0 else xmax - self.bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) + if hgrid is None and bbox is None: + raise ValueError('Either hgrid or bbox must be provided!') - def fetch_data(self, date, fmt='schism'): + if hgrid is not None: + xmin, xmax = hgrid.coords[:, 0].min(), hgrid.coords[:, 0].max() + ymin, ymax = hgrid.coords[:, 1].min(), hgrid.coords[:, 1].max() + self.bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) + elif bbox is not None: + self.bbox = bbox + + def fetch_data(self, start_date, rnday=1, fmt='schism', bnd=False, nudge=False, outdir=None): ''' - date: datetime.datetime + start_date: datetime.datetime + rnday: integer fmt: 'schism' - for Fortran code; 'hycom' - raw netCDF from HYCOM + bnd: file names are SSH_*.nc, TS_*.nc, UV_*.nc used in gen_hot_3Dth_from_hycom.f90 + nudge: file name is TS_*.nc used in gen_nudge_from_hycom.f90 + outdir: directory for output files ''' + timevector = np.arange( + start_date, + start_date + timedelta(days=rnday+1), + timedelta(days=1) + ).astype(datetime) - database=get_database(date) - logger.info(f'Fetching data for {date} from database {database}') - - time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 = get_idxs(date, database, self.bbox) - - if fmt == 'schism': - url_ssh = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ - f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \ - f'surf_el[{time_idx}][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' - foutname = f'SSH_{date.strftime("%Y%m%d")}.nc' - logger.info(f'filename is {foutname}') - ds = xr.open_dataset(url_ssh) - ds = convert_longitude(ds) - ds = ds.rename_dims({'lon':'xlon'}) - ds = ds.rename_dims({'lat':'ylat'}) - ds = ds.rename_vars({'lat':'ylat'}) - ds = ds.rename_vars({'lon':'xlon'}) - ds.to_netcdf(foutname, 'w', 'NETCDF3_CLASSIC', unlimited_dims='time') - ds.close() - - url_uv = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ - f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \ - f'water_u[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ - f'water_v[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' - - foutname = f'UV_{date.strftime("%Y%m%d")}.nc' - logger.info(f'filename is {foutname}') - ds = xr.open_dataset(url_uv) - ds = convert_longitude(ds) - ds = ds.rename_dims({'lon':'xlon'}) - ds = ds.rename_dims({'lat':'ylat'}) - ds = ds.rename_vars({'lat':'ylat'}) - ds = ds.rename_vars({'lon':'xlon'}) - ds.to_netcdf(foutname, 'w', 'NETCDF3_CLASSIC', unlimited_dims='time') - ds.close() - - url_ts = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ - f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \ - f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ - f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' + for i, date in enumerate(timevector): + database=get_database(date) + logger.info(f'Fetching data for {date} from database {database}') - foutname = f'TS_{date.strftime("%Y%m%d")}.nc' - logger.info(f'filename is {foutname}') + time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 = get_idxs(date, database, self.bbox) - ds = xr.open_dataset(url_ts) - - temp = ds.water_temp.values - salt = ds.salinity.values - dep = ds.depth.values - - #convert in-situ temperature to potential temperature - ptemp = ConvertTemp(salt, temp, dep) + if fmt == 'schism': + url_ssh = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ + f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}:8:{time_idx}],' + \ + f'surf_el[{time_idx}][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ + f'water_u[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ + f'water_v[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ + f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ + f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' - #drop water_temp variable and add new temperature variable - ds = ds.drop('water_temp') - ds['temperature']=(['time','depth','lat','lon'], ptemp) - ds.temperature.attrs = { - 'long_name': 'Sea water potential temperature', - 'standard_name': 'sea_water_potential_temperature', - 'units': 'degC' - } - - #ds.assign(water_temp2=ptemp) - #ds.assign.attrs = ds.water_temp.attrs - - ds = convert_longitude(ds) - ds = ds.rename_dims({'lon':'xlon'}) - ds = ds.rename_dims({'lat':'ylat'}) - ds = ds.rename_vars({'lat':'ylat'}) - ds = ds.rename_vars({'lon':'xlon'}) - ds.to_netcdf(foutname, 'w', unlimited_dims='time', encoding={'temperature':{'dtype': 'h', '_FillValue': -30000.,'scale_factor': 0.001, 'add_offset': 20., 'missing_value': -30000.}}) - ds.close() + foutname = f'hycom_{date.strftime("%Y%m%d")}.nc' + #foutname = f'TS_{i+1}.nc' + logger.info(f'filename is {foutname}') + ds = xr.open_dataset(url_ssh) + + #convert in-situ temperature to potential temperature + temp = ds.water_temp.values + salt = ds.salinity.values + dep = ds.depth.values + + ptemp = ConvertTemp(salt, temp, dep) + #drop water_temp variable and add new temperature variable + ds = ds.drop('water_temp') + ds['temperature']=(['time','depth','lat','lon'], ptemp) + ds.temperature.attrs = { + 'long_name': 'Sea water potential temperature', + 'standard_name': 'sea_water_potential_temperature', + 'units': 'degC' + } + + ds = convert_longitude(ds) + ds = ds.rename_dims({'lon':'xlon'}) + ds = ds.rename_dims({'lat':'ylat'}) + ds = ds.rename_vars({'lat':'ylat'}) + ds = ds.rename_vars({'lon':'xlon'}) + ds.to_netcdf(foutname, 'w', unlimited_dims='time', encoding={'temperature':{'dtype': 'h', '_FillValue': -30000.,'scale_factor': 0.001, 'add_offset': 20., 'missing_value': -30000.}}) + ds.close() - elif fmt == 'hycom': - url=f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ - f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \ - f'surf_el[{time_idx}][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ - f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ - f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ - f'water_u[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ - f'water_v[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' + #link output file to Fortran required files + dir_path = os.path.abspath(outdir) + logger.info(f'current dir is {dir_path}') + src = f'{dir_path}/{foutname}' + if bnd: + names = ['SSH', 'TS', 'UV'] + for name in names: + dst = f'{dir_path}/{name}_{i+1}.nc' + os.symlink(src, dst) + elif nudge: + dst = f'{dir_path}/TS_{i+1}.nc' + os.symlink(src, dst) + + elif fmt == 'hycom': + url=f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ + f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \ + f'surf_el[{time_idx}][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ + f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ + f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ + f'water_u[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ + f'water_v[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' - foutname = f'hycom_{date.strftime("%Y%m%d")}.nc' - ds = xr.open_dataset(url) - ds.to_netcdf(foutname, 'w') + foutname = f'hycom_{date.strftime("%Y%m%d")}.nc' + ds = xr.open_dataset(url) + ds.to_netcdf(foutname, 'w') - ds.close() + ds.close() From 5e44d2dc23a65899acc74b3943417bef7d734d93 Mon Sep 17 00:00:00 2001 From: cuill Date: Wed, 8 Feb 2023 11:34:53 -0500 Subject: [PATCH 17/31] Check if lon needs to be converted --- pyschism/forcing/hycom/hycom2schism.py | 41 +++++++++++--------------- 1 file changed, 17 insertions(+), 24 deletions(-) diff --git a/pyschism/forcing/hycom/hycom2schism.py b/pyschism/forcing/hycom/hycom2schism.py index 8408dff6..eb4cfd1d 100644 --- a/pyschism/forcing/hycom/hycom2schism.py +++ b/pyschism/forcing/hycom/hycom2schism.py @@ -75,7 +75,9 @@ def get_idxs(date, database, bbox, lonc=None, latc=None): dep=ds['depth'][:] #check if hycom's lon is the same range as schism's + same = True if not (bbox.xmin >= lon.min() and bbox.xmax <= lon.max()): + same = False if lon.min() >= 0: logger.info(f'Convert HYCOM longitude from [0, 360) to [-180, 180):') idxs = lon>=180 @@ -100,10 +102,10 @@ def get_idxs(date, database, bbox, lonc=None, latc=None): if lonc is None: lonc = lon.mean() - logger.info(f'lonc is {lonc}') + #logger.info(f'lonc is {lonc}') if latc is None: latc = lat.mean() - logger.info(f'latc is {latc}') + #logger.info(f'latc is {latc}') x2, y2=transform_ll_to_cpp(lon, lat, lonc, latc) idxs=np.where( date == times)[0] @@ -125,7 +127,7 @@ def get_idxs(date, database, bbox, lonc=None, latc=None): ds.close() - return time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 + return time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2, same def transform_ll_to_cpp(lon, lat, lonc=-77.07, latc=24.0): longitude=lon/180*np.pi @@ -385,17 +387,7 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, ocean_b ymin, ymax = np.min(blat), np.max(blat) bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) - #if date >= datetime(2017, 2, 1) and \ - # date < datetime(2017, 6, 1) or \ - # date >= datetime(2017, 10, 1): - # xmin = xmin + 360. if xmin < 0 else xmin - # xmax = xmax + 360. if xmax < 0 else xmax - # bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) - #else: - # bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) - #logger.info(f'xmin is {xmin}, xmax is {xmax}') - - time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 = get_idxs(date, database, bbox, lonc=blonc, latc=blatc) + time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2, _ = get_idxs(date, database, bbox, lonc=blonc, latc=blatc) if date >= datetime.utcnow(): date2 = datetime.utcnow() - timedelta(days=1) @@ -607,7 +599,6 @@ def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, #compute zcor zcor = depth[:,None]*sigma nvrt=zcor.shape[1] - #logger.info(f'zcor at node 1098677 is {zcor[1098676,:]}') #Get open nudge array nlon = hgrid.coords[include, 0] @@ -619,9 +610,7 @@ def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, zcor2=zcor[include,:] idxs=np.where(zcor2 > 5000) - #logger.info(idxs) zcor2[idxs]=5000.0-1.0e-6 - #logger.info(f'zcor2 at node 200 is {zcor2[199,:]}') #construct schism grid x2i=np.tile(xi,[nvrt,1]).T @@ -697,7 +686,7 @@ def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) #logger.info(f'xmin is {xmin}, xmax is {xmax}') - time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 = get_idxs(date, database, bbox, lonc=nlonc, latc=nlatc) + time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2, _ = get_idxs(date, database, bbox, lonc=nlonc, latc=nlatc) if date >= datetime.utcnow(): date2 = datetime.utcnow() - timedelta(days=1) @@ -770,17 +759,18 @@ def fetch_data(self, start_date, rnday=1, fmt='schism', bnd=False, nudge=False, nudge: file name is TS_*.nc used in gen_nudge_from_hycom.f90 outdir: directory for output files ''' - timevector = np.arange( - start_date, - start_date + timedelta(days=rnday+1), - timedelta(days=1) + if rnday == 1: + timevector = [start_date] + else: + timevector = np.arange( + start_date, start_date + timedelta(days=rnday+1), timedelta(days=1) ).astype(datetime) for i, date in enumerate(timevector): database=get_database(date) logger.info(f'Fetching data for {date} from database {database}') - time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2 = get_idxs(date, database, self.bbox) + time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2, isLonSame = get_idxs(date, database, self.bbox) if fmt == 'schism': url_ssh = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ @@ -811,7 +801,10 @@ def fetch_data(self, start_date, rnday=1, fmt='schism', bnd=False, nudge=False, 'units': 'degC' } - ds = convert_longitude(ds) + if not isLonSame: + logger.info('Lon is not the same!') + ds = convert_longitude(ds) + ds = ds.rename_dims({'lon':'xlon'}) ds = ds.rename_dims({'lat':'ylat'}) ds = ds.rename_vars({'lat':'ylat'}) From 5756a047e60fffb65b40977902351f898ad05d0f Mon Sep 17 00:00:00 2001 From: cuill Date: Thu, 2 Mar 2023 11:24:19 -0500 Subject: [PATCH 18/31] fixed a bug in nudge when the last node sits on the bnd --- pyschism/forcing/hycom/hycom2schism.py | 11 +++++++---- pyschism/mesh/base.py | 6 ++++++ 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/pyschism/forcing/hycom/hycom2schism.py b/pyschism/forcing/hycom/hycom2schism.py index eb4cfd1d..cfb8d74a 100644 --- a/pyschism/forcing/hycom/hycom2schism.py +++ b/pyschism/forcing/hycom/hycom2schism.py @@ -528,11 +528,14 @@ def compute_nudge(lon, lat, nnode, opbd, out): idxs=np.where(out > 0)[0] idxs_nudge[idxs]=1 + #expand nudging marker to neighbor nodes - idxs=np.where(np.max(out[elnode], axis=1) > 0)[0] - fp=elnode[idxs,-1] < 0 - idxs_nudge[elnode[idxs[fp],:3]]=1 - idxs_nudge[elnode[idxs[~fp],:]]=1 + i34 = hgrid.elements.i34 + fp = i34==3 + idxs=np.where(np.max(out[elnode[fp, 0:3]], axis=1) > 0)[0] + idxs_nudge[elnode[fp,0:3][idxs,:]]=1 + idxs=np.where(np.max(out[elnode[~fp, :]], axis=1) > 0)[0] + idxs_nudge[elnode[~fp,:][idxs,:]]=1 #idxs_nudge=np.delete(idxs_nudge, np.where(idxs_nudge == -99)) idxs=np.where(idxs_nudge == 1)[0] diff --git a/pyschism/mesh/base.py b/pyschism/mesh/base.py index ff90cf00..1aa72db1 100644 --- a/pyschism/mesh/base.py +++ b/pyschism/mesh/base.py @@ -361,6 +361,12 @@ def array(self): self._array = array return self._array + @property + def i34(self): + if not hasattr(self, "_i34"): + self._i34 = np.sum(~self.array.mask, axis=1) + return self._i34 + @property def triangles(self): if not hasattr(self, "_triangles"): From c4237e01a315a64800469d83653ca555e220b65b Mon Sep 17 00:00:00 2001 From: cuill Date: Thu, 2 Mar 2023 15:31:55 -0500 Subject: [PATCH 19/31] check if hgrid has quads --- pyschism/mesh/base.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pyschism/mesh/base.py b/pyschism/mesh/base.py index 1aa72db1..ee1f9681 100644 --- a/pyschism/mesh/base.py +++ b/pyschism/mesh/base.py @@ -364,7 +364,10 @@ def array(self): @property def i34(self): if not hasattr(self, "_i34"): - self._i34 = np.sum(~self.array.mask, axis=1) + if np.any(self.array.mask): + self._i34 = np.sum(~self.array.mask, axis=1) + else: + self._i34 = np.full(self.array.shape[0], 3) return self._i34 @property From 26d87046d61c8cd55b03783b67193add25306d6e Mon Sep 17 00:00:00 2001 From: cuill Date: Thu, 2 Mar 2023 15:32:38 -0500 Subject: [PATCH 20/31] fixed a bug in mapping global nodes --- pyschism/forcing/hycom/hycom2schism.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyschism/forcing/hycom/hycom2schism.py b/pyschism/forcing/hycom/hycom2schism.py index cfb8d74a..ebf07ff0 100644 --- a/pyschism/forcing/hycom/hycom2schism.py +++ b/pyschism/forcing/hycom/hycom2schism.py @@ -645,7 +645,7 @@ def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, #dst_temp['time'][:] = ndt dst_temp.createVariable('map_to_global_node', 'i4', ('node',)) - dst_temp['map_to_global_node'][:] = include+0 + dst_temp['map_to_global_node'][:] = include+1 dst_temp.createVariable('tracer_concentration', 'f', ('time', 'node', 'nLevels', 'one')) #dst_temp['tracer_concentration'][:,:,:,:] = timeseries_t @@ -662,7 +662,7 @@ def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, #dst_salt['time'][:] = ndt dst_salt.createVariable('map_to_global_node', 'i4', ('node',)) - dst_salt['map_to_global_node'][:] = include+0 + dst_salt['map_to_global_node'][:] = include+1 dst_salt.createVariable('tracer_concentration', 'f', ('time', 'node', 'nLevels', 'one')) #dst_salt['tracer_concentration'][:,:,:,:] = timeseries_s From 0850401e52f13aa3ee64edc863c7251c81b4efa3 Mon Sep 17 00:00:00 2001 From: cuill Date: Fri, 3 Mar 2023 16:53:19 -0500 Subject: [PATCH 21/31] Added support for FES2014 --- pyschism/forcing/bctides/fes2014.py | 155 ++++++++++++++++++++++++++++ pyschism/forcing/bctides/tides.py | 4 +- 2 files changed, 158 insertions(+), 1 deletion(-) create mode 100644 pyschism/forcing/bctides/fes2014.py diff --git a/pyschism/forcing/bctides/fes2014.py b/pyschism/forcing/bctides/fes2014.py new file mode 100644 index 00000000..aad2ae73 --- /dev/null +++ b/pyschism/forcing/bctides/fes2014.py @@ -0,0 +1,155 @@ +import logging +import os +import pathlib +import glob + +import appdirs +from netCDF4 import Dataset +import numpy as np +from scipy.interpolate import griddata +from scipy.interpolate.fitpack2 import RectBivariateSpline + +from pyschism.forcing.bctides.base import TidalDataProvider + + +logger = logging.getLogger(__name__) + +FES2014_TIDES_EXTRA = 'ocean_tide_extrapolated' +FES2014_EASTWARD_VEL = 'eastward_velocity' +FES2014_NORTHWARD_VEL = 'northward_velocity' + + +def raise_missing_file(fpath, fname): + raise FileNotFoundError('\n'.join([ + f'No FES2014 file found at "{fpath}".', + 'New users will need to register and download a copy of ' + f'the FES2014 NetCDF file (specifically `{fname}`) ' + 'from https://www.aviso.altimetry.fr/en/data/products/auxiliary-products/global-tide-fes.html', + 'Once you obtain netCDF files, you can follow one of the ' + 'following options: ', + f'1) copy or symlink the file to "{fpath}"', + f'2) set the environment variable `{fname}` to point' + ' to the file', + ])) + + +class FES2014(TidalDataProvider): + + def __init__(self, resource=None): + self.resource = resource + + def get_elevation(self, constituent, vertices): + logger.info('Querying FES2014 for elevation constituent ' + f'{constituent}.') + amp = self._get_interpolation( + 'elevation', 'amplitude', constituent, vertices) + phase = self._get_interpolation( + 'elevation', 'phase', constituent, vertices) + return amp, phase + + def get_velocity(self, constituent, vertices): + logger.info('Querying FES2014 for velocity constituent ' + f'{constituent}.') + uamp = self._get_interpolation( + 'eastward_vel', 'Ua', constituent, vertices) + uphase = self._get_interpolation( + 'eastward_vel', 'Ug', constituent, vertices) + vamp = self._get_interpolation( + 'northward_vel', 'Va', constituent, vertices) + vphase = self._get_interpolation( + 'northward_vel', 'Vg', constituent, vertices) + return uamp, uphase, vamp, vphase + + @property + def constituents(self): + return ['M2', 'S2', 'N2', 'K2', 'K1', 'O1', 'P1', + 'Q1', 'Mm', 'Mf', 'M4', 'MN4', 'MS4', '2N2', 'S1'] + if not hasattr(self, '_constituents'): + files = glob.glob(appdirs.user_data_dir('fes2014') + '/ocean_tide_extrapolated/*.nc') + self._constituents = [f.split('/')[-1].split('.')[0] for f in files] + return self._constituents + + def _get_resource(self, variable, constituent) -> Dataset: + resource = self._resource[variable][constituent] + if resource is not None: + return Dataset(resource) + datapath = pathlib.Path(appdirs.user_data_dir('fes2014')) + if variable == 'elevation': + fname = datapath / f'{FES2014_TIDES_EXTRA}/{constituent.lower()}.nc' + if variable == 'eastward_vel': + fname = datapath / f'{FES2014_EASTWARD_VEL}/{constituent.lower()}.nc' + if variable == 'northward_vel': + fname = datapath / f'{FES2014_NORTHWARD_VEL}/{constituent.lower()}.nc' + return Dataset(fname) + + def _get_interpolation(self, phys_var, ncvar, constituent, vertices): + ds = self._get_resource(phys_var, constituent) + lon = ds['lon'][:] + lat = ds['lat'][:] + dxs = np.unique(np.diff(lon)) + dys = np.unique(np.diff(lat)) + if len(dxs) !=1 or len(dys) != 1: + raise ValueError(f'{phys_var}: lon, lat of {constituent}.nc not uniform! ') + dx = dxs[0] + dy = dys[0] + + #get interp index + xi = np.asarray( + [x + 360. if x < 0. else x for x in vertices[:, 0]]).flatten() + yi = vertices[:, 1].flatten() + idx = np.floor((xi - lon[0]) / dx).astype('int') + mask = np.nonzero((lon[idx] - xi) > 0)[0] + idx[mask] = idx[mask] - 1 + idy = np.floor((yi - lat[0]) / dy).astype('int') + mask = np.nonzero((lat[idy] - yi) > 0)[0] + idy[mask] = idy[mask] - 1 + xrat = (xi - lon[idx]) / (lon[idx+1] - lon[idx]) + yrat = (yi - lat[idy]) / (lat[idy+1] - lat[idy]) + if np.sum((xrat > 1) | (xrat < 0) | (yrat > 1) | (yrat < 0)) != 0: + raise ValueError(f'xrat or yrat > 1 or < 0') + + zi = ds[ncvar][:,:] + # vm = 100 junk for amplitude, vm = 370 junk for phase + if ncvar == 'amplitude' or ncvar == 'Ua' or ncvar == 'Va': + vm = 100 + zi = zi / 100 + if ncvar == 'phase' or ncvar == 'Ug' or ncvar == 'Vg': + vm = 370 + #Put phases into [0, 360) + mask = zi<0 + zi[mask] = zi[mask] + 360 + v0 = np.c_[zi[idy, idx], zi[idy, idx+1], zi[idy+1, idx], zi[idy+1, idx+1]].T + vmax = v0.max(axis=0) + vmin = v0.min(axis=0) + #handle phase jump + if ncvar == 'phase' or ncvar == 'Ug' or ncvar == 'Vg': + for kk in np.nonzero(abs(vmax - vmin) > 180)[0]: + mask = abs(v0[:, kk] - vmax[kk]) > 180 + v0[mask, kk] = v0[mask, kk] + 360 + v1 = v0[0] * (1 - xrat) + v0[1] * xrat + v2 = v0[2] * (1 - xrat) + v0[3] * xrat + values = v1 * (1 - yrat) + v2 * yrat + mask = np.nonzero((vmax > vm) * (vmin <= vm) * (vmin >= 0))[0] + values[mask] = vmin[mask] + if sum((vmax > vm) * ((vmin > vm) | (vmin < 0))) != 0: + raise ValueError('All junk values for {phys_var} {constituent}') + + ds.close() + return values + + @property + def resource(self): + return self._resource + + @resource.setter + def resource(self, resource): + if resource is None: + resource = {'elevation': {}, 'eastward_vel': {}, 'northward_vel': {}} + for constituent in self.constituents: + for key in resource.keys(): + resource[key].update({constituent: None}) + else: + raise NotImplementedError('Check that static files exist.') + for file in pathlib.Path(resource).glob('*.nc'): + print(file) + self._resource = resource diff --git a/pyschism/forcing/bctides/tides.py b/pyschism/forcing/bctides/tides.py index 65503133..0fa07ac6 100644 --- a/pyschism/forcing/bctides/tides.py +++ b/pyschism/forcing/bctides/tides.py @@ -9,6 +9,7 @@ # from pyschism.forcing.tides import bctypes from pyschism.forcing.bctides.tpxo import TPXO +from pyschism.forcing.bctides.fes2014 import FES2014 from pyschism.forcing.bctides.hamtide import HAMTIDE # from pyschism.forcing.baroclinic.gofs import GOFS # from pyschism.forcing.baroclinic.rtofs import RTOFS @@ -22,6 +23,7 @@ class TidalDatabase(Enum): TPXO = TPXO + FES2014 = FES2014 HAMTIDE = HAMTIDE def _missing_(self, name): @@ -573,7 +575,7 @@ def active_constituents(self): def all_constituents(self): if isinstance(self.tidal_database, HAMTIDE): return self.major_constituents - elif isinstance(self.tidal_database, TPXO): + elif isinstance(self.tidal_database, TPXO) or isinstance(self.tidal_database, FES2014): return (*self.major_constituents, *self.minor_constituents) else: raise NotImplementedError( From e22eb6c91478ceba6dda6ecd60e11b7a81da6c4a Mon Sep 17 00:00:00 2001 From: cuill Date: Tue, 14 Mar 2023 14:36:29 -0400 Subject: [PATCH 22/31] enable nudging over multiple ocean boundaries --- pyschism/forcing/hycom/hycom2schism.py | 256 +++++++++++++------------ 1 file changed, 132 insertions(+), 124 deletions(-) diff --git a/pyschism/forcing/hycom/hycom2schism.py b/pyschism/forcing/hycom/hycom2schism.py index ebf07ff0..5c30e6cb 100644 --- a/pyschism/forcing/hycom/hycom2schism.py +++ b/pyschism/forcing/hycom/hycom2schism.py @@ -472,74 +472,66 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, ocean_b class Nudge: - def __init__(self, ocean_bnd_ids=[0]): - - self.include = None - self.ocean_bnd_ids = ocean_bnd_ids + def __init__(self, hgrid=None, ocean_bnd_ids=None): + if hgrid is None: + raise ValueError('No hgrid information!') + else: + self.hgrid = hgrid - def gen_nudge(self, outdir: Union[str, os.PathLike], hgrid, rlmax = 1.5, rnu_day=0.25): + + if ocean_bnd_ids is None: + raise ValueError('Please specify indexes for ocean boundaries!') + else: + self.ocean_bnd_ids = ocean_bnd_ids - @jit(nopython=True, parallel=True) - def compute_nudge(lon, lat, nnode, opbd, out): - - rnu_max = 1.0 / rnu_day / 86400.0 - rnu = 0 - for idn in prange(nnode): - if idn in opbd: - rnu = rnu_max - distmin = 0. - else: - distmin = np.finfo(np.float64).max - for j in opbd: - tmp = np.square(lon[idn]-lon[j]) + np.square(lat[idn]-lat[j]) - rl2 = np.sqrt(tmp) - if rl2 < distmin: - distmin=rl2 - rnu = 0. - if distmin <= rlmax: - rnu = (1-distmin/rlmax)*rnu_max - #idxs_nudge[idn]=1 #idn - out[idn] = rnu + def gen_nudge(self, outdir: Union[str, os.PathLike], rlmax = 1.5, rnu_day=0.25): + outdir = pathlib.Path(outdir) - #get nudge zone - lon=hgrid.coords[:,0] - lat=hgrid.coords[:,1] + rnu_max = 1.0 / rnu_day / 86400.0 - #Get open boundary - gdf=hgrid.boundaries.open.copy() - opbd=[] - #for boundary in gdf.itertuples(): - for ibnd in self.ocean_bnd_ids: - opbd.extend(list(gdf.iloc[ibnd].indexes)) - opbd = np.array(opbd) + #get nudge zone + lon = self.hgrid.coords[:,0] + lat = self.hgrid.coords[:,1] + gdf = self.hgrid.boundaries.open.copy() + elnode = self.hgrid.elements.array + NE, NP = elnode.shape[0],len(lon) + nudge_coeff = np.zeros(NP, dtype=float) - elnode=hgrid.elements.array - NE, NP = len(elnode), len(lon) + global_idxs = {} - out = np.zeros([NP]) - idxs_nudge=np.zeros(NP, dtype=int) t0 = time() - #compute_nudge(lon, lat, NP, opbd2, idxs_nudge, out) - compute_nudge(lon, lat, NP, opbd, out) + for i in self.ocean_bnd_ids: + print(f'boundary {i}') + bnd_idxs = gdf.iloc[i].indexes + + dis = abs((lon + 1j*lat)[:, None] - (lon[bnd_idxs] + 1j*lat[bnd_idxs])[None, :]).min(axis=1) + out = (1-dis/rlmax)*rnu_max + out[out<0] = 0 + out[out>rnu_max] = rnu_max + fp = out>0 + + nudge_coeff[fp] = -out[fp] - idxs=np.where(out > 0)[0] - idxs_nudge[idxs]=1 + idxs_nudge=np.zeros(NP, dtype=int) + idxs=np.where(out > 0)[0] + idxs_nudge[idxs]=1 - #expand nudging marker to neighbor nodes - i34 = hgrid.elements.i34 - fp = i34==3 - idxs=np.where(np.max(out[elnode[fp, 0:3]], axis=1) > 0)[0] - idxs_nudge[elnode[fp,0:3][idxs,:]]=1 - idxs=np.where(np.max(out[elnode[~fp, :]], axis=1) > 0)[0] - idxs_nudge[elnode[~fp,:][idxs,:]]=1 - - #idxs_nudge=np.delete(idxs_nudge, np.where(idxs_nudge == -99)) - idxs=np.where(idxs_nudge == 1)[0] - self.include=idxs + #expand nudging marker to neighbor nodes + i34 = self.hgrid.elements.i34 + fp = i34==3 + idxs=np.where(np.max(out[elnode[fp, 0:3]], axis=1) > 0)[0] + idxs_nudge[elnode[fp,0:3][idxs,:]]=1 + idxs=np.where(np.max(out[elnode[~fp, :]], axis=1) > 0)[0] + idxs_nudge[elnode[~fp,:][idxs,:]]=1 + + idxs=np.where(idxs_nudge == 1)[0] + global_idxs[i] = idxs + + #logger.info(f'len of nudge idxs is {len(idxs)}') logger.info(f'It took {time() -t0} sencods to calcuate nudge coefficient') @@ -547,13 +539,13 @@ def compute_nudge(lon, lat, nnode, opbd, out): nudge.extend("\n") nudge.append(f"{NE} {NP}") nudge.extend("\n") - hgrid = hgrid.to_dict() + hgrid = self.hgrid.to_dict() nodes = hgrid['nodes'] elements = hgrid['elements'] for idn, (coords, values) in nodes.items(): line = [f"{idn}"] line.extend([f"{x:<.7e}" for x in coords]) - line.extend([f"{out[int(idn)-1]:<.7e}"]) + line.extend([f"{nudge_coeff[int(idn)-1]:<.7e}"]) line.extend("\n") nudge.append(" ".join(line)) @@ -569,9 +561,9 @@ def compute_nudge(lon, lat, nnode, opbd, out): shutil.copy2(outdir / 'TEM_nudge.gr3', outdir / 'SAL_nudge.gr3') - return self.include + return global_idxs - def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, rnday, restart = False, rlmax = None, rnu_day=None): + def fetch_data(self, outdir: Union[str, os.PathLike], vgrid, start_date, rnday, restart = False, rlmax = None, rnu_day=None): outdir = pathlib.Path(outdir) @@ -590,41 +582,31 @@ def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, rnu_day = 0.25 if rnu_day is None else rnu_day logger.info(f'Max relax distance is {rlmax} degree, max relax strengh is {rnu_day} days.') #Get the index for nudge - include = self.gen_nudge(outdir,hgrid, rlmax = rlmax, rnu_day=rnu_day) + global_idxs = self.gen_nudge(outdir, rlmax = rlmax, rnu_day=rnu_day) - #get coords of SCHISM - loni=hgrid.nodes.coords[:,0] - lati=hgrid.nodes.coords[:,1] #get bathymetry - depth = hgrid.values + depth = self.hgrid.values #compute zcor zcor = depth[:,None]*sigma nvrt=zcor.shape[1] - #Get open nudge array - nlon = hgrid.coords[include, 0] - nlat = hgrid.coords[include, 1] - nlonc = nlon.mean() - nlatc = nlat.mean() - xi,yi = transform_ll_to_cpp(nlon, nlat, nlonc, nlatc) - bxy = np.c_[yi, xi] - - zcor2=zcor[include,:] - idxs=np.where(zcor2 > 5000) - zcor2[idxs]=5000.0-1.0e-6 - - #construct schism grid - x2i=np.tile(xi,[nvrt,1]).T - y2i=np.tile(yi,[nvrt,1]).T - bxyz=np.c_[zcor2.reshape(np.size(zcor2)),y2i.reshape(np.size(y2i)),x2i.reshape(np.size(x2i))] - logger.info('Computing SCHISM zcor is done!') #allocate output variables - nNode=len(include) - one=1 - ntimes=self.rnday+1 + include = global_idxs[0] + nbnd = len(self.ocean_bnd_ids) + if nbnd > 1: + include = np.concatenate((global_idxs[0], global_idxs[1],)) + if nbnd > 2: + for i in self.ocean_bnd_ids[2:]: + include = np.concatenate((include, global_idxs[i])) + else: + include = global_idxs[0] + + nNode = include.shape[0] + one = 1 + ntimes = self.rnday+1 #timeseries_s=np.zeros([ntimes,nNode,nvrt,one]) #timeseries_t=np.zeros([ntimes,nNode,nvrt,one]) @@ -684,55 +666,81 @@ def fetch_data(self, outdir: Union[str, os.PathLike], hgrid, vgrid, start_date, database=get_database(date) logger.info(f'Fetching data for {date} from database {database}') - xmin, xmax = np.min(nlon), np.max(nlon) - ymin, ymax = np.min(nlat), np.max(nlat) - bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) - #logger.info(f'xmin is {xmin}, xmax is {xmax}') - time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2, _ = get_idxs(date, database, bbox, lonc=nlonc, latc=nlatc) + ind1 = 0 + ind2 = 0 + for ibnd in self.ocean_bnd_ids: + include = global_idxs[ibnd] - if date >= datetime.utcnow(): - date2 = datetime.utcnow() - timedelta(days=1) - url = f'https://tds.hycom.org/thredds/dodsC/{database}/FMRC/runs/GLBy0.08_930_FMRC_RUN_' + \ - f'{date2.strftime("%Y-%m-%dT12:00:00Z")}?depth[0:1:-1],lat[{lat_idx1}:1:{lat_idx2}],' + \ - f'lon[{lon_idx1}:1:{lon_idx2}],time[{time_idx}],' + \ - f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ - f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' + ind1 = ind2 + ind2 = ind1 + include.shape[0] + #Get open nudge array + nlon = self.hgrid.coords[include, 0] + nlat = self.hgrid.coords[include, 1] + nlonc = nlon.mean() + nlatc = nlat.mean() + xi,yi = transform_ll_to_cpp(nlon, nlat, nlonc, nlatc) + bxy = np.c_[yi, xi] - else: - url=f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ - f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \ - f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ - f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' - #logger.info(url) + zcor2=zcor[include,:] + idxs=np.where(zcor2 > 5000) + zcor2[idxs]=5000.0-1.0e-6 - ds=Dataset(url) - salt=np.squeeze(ds['salinity'][:,:,:]) - temp=np.squeeze(ds['water_temp'][:,:,:]) - #logger.info(f'The shape of temp is {temp.shape}') + #construct schism grid + x2i=np.tile(xi,[nvrt,1]).T + y2i=np.tile(yi,[nvrt,1]).T + bxyz=np.c_[zcor2.reshape(np.size(zcor2)),y2i.reshape(np.size(y2i)),x2i.reshape(np.size(x2i))] + logger.info('Computing SCHISM zcor is done!') - #Convert temp to potential temp - dep=ds['depth'][:] - ptemp = ConvertTemp(salt, temp, dep) + xmin, xmax = np.min(nlon), np.max(nlon) + ymin, ymax = np.min(nlat), np.max(nlat) + bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) + #logger.info(f'xmin is {xmin}, xmax is {xmax}') - logger.info('****Interpolation starts****') + time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2, _ = get_idxs(date, database, bbox, lonc=nlonc, latc=nlatc) - #ndt[it]=it - #salt - dst_salt['time'][it] = it - salt_int = interp_to_points_3d(dep, y2, x2, bxyz, salt) - salt_int = salt_int.reshape(zcor2.shape) - #timeseries_s[it,:,:,0]=salt_int - dst_salt['tracer_concentration'][it,:,:,0] = salt_int + if date >= datetime.utcnow(): + date2 = datetime.utcnow() - timedelta(days=1) + url = f'https://tds.hycom.org/thredds/dodsC/{database}/FMRC/runs/GLBy0.08_930_FMRC_RUN_' + \ + f'{date2.strftime("%Y-%m-%dT12:00:00Z")}?depth[0:1:-1],lat[{lat_idx1}:1:{lat_idx2}],' + \ + f'lon[{lon_idx1}:1:{lon_idx2}],time[{time_idx}],' + \ + f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ + f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' - #temp - dst_temp['time'][it] = it - temp_int = interp_to_points_3d(dep, y2, x2, bxyz, ptemp) - temp_int = temp_int.reshape(zcor2.shape) - #timeseries_t[it,:,:,0]=temp_int - dst_temp['tracer_concentration'][it,:,:,0] = temp_int + else: + url=f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ + f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \ + f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ + f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' + #logger.info(url) + + ds=Dataset(url) + salt=np.squeeze(ds['salinity'][:,:,:]) + temp=np.squeeze(ds['water_temp'][:,:,:]) + #logger.info(f'The shape of temp is {temp.shape}') - ds.close() + #Convert temp to potential temp + dep=ds['depth'][:] + ptemp = ConvertTemp(salt, temp, dep) + + logger.info('****Interpolation starts****') + + #ndt[it]=it + #salt + dst_salt['time'][it] = it + salt_int = interp_to_points_3d(dep, y2, x2, bxyz, salt) + salt_int = salt_int.reshape(zcor2.shape) + #timeseries_s[it,:,:,0]=salt_int + dst_salt['tracer_concentration'][it,ind1:ind2,:,0] = salt_int + + #temp + dst_temp['time'][it] = it + temp_int = interp_to_points_3d(dep, y2, x2, bxyz, ptemp) + temp_int = temp_int.reshape(zcor2.shape) + #timeseries_t[it,:,:,0]=temp_int + dst_temp['tracer_concentration'][it,ind1:ind2,:,0] = temp_int + + ds.close() #dst_temp.close() #dst_salt.close() From ab39ef6a1e2a77754c314c7ae20a9e8ce0e6ceb1 Mon Sep 17 00:00:00 2001 From: cuill Date: Tue, 14 Mar 2023 17:44:12 -0400 Subject: [PATCH 23/31] Remove longitude conversion in ERA5 --- pyschism/forcing/nws/nws2/era5.py | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/pyschism/forcing/nws/nws2/era5.py b/pyschism/forcing/nws/nws2/era5.py index 5da23da8..f313d79d 100644 --- a/pyschism/forcing/nws/nws2/era5.py +++ b/pyschism/forcing/nws/nws2/era5.py @@ -54,7 +54,8 @@ def __init__(self, start_date=None, rnday: Union[float, timedelta] = 4, bbox=Non }) filename = self.tmpdir / f"era5_{self.start_date.strftime('%Y%m%d')}.nc" - r.download(filename) + if filename.is_file() == False: + r.download(filename) @property def tmpdir(self): @@ -87,13 +88,7 @@ def lat(self): return self._lat def xy_grid(self): - lon = [] - for x in self.lon: - if x > 180: - lon.append(x-360) - else: - lon.append(x) - return np.meshgrid(np.array(lon), self.lat) + return np.meshgrid(self.lon, self.lat) def put_sflux_fields(iday, date, timevector, ds, nx_grid, ny_grid, air, rad, prc, output_interval, OUTDIR): rt=pd.to_datetime(str(date)) From 8c91c7fe90fe3b943c994d3260329b87a192fc2f Mon Sep 17 00:00:00 2001 From: cuill Date: Wed, 15 Mar 2023 11:30:37 -0400 Subject: [PATCH 24/31] update example scripts --- examples/Gr3/Drag/GoME+0.001.reg | 23 ++++ examples/Gr3/Drag/Lake_Charles_0.reg | 20 ++++ .../01_gen_drag.py => Gr3/Drag/gen_drag.py} | 8 +- examples/Gr3/Shapiro/coastal_0.2.cpp.reg | 100 ++++++++++++++++++ examples/Gr3/Shapiro/coastal_0.5_1.cpp.reg | 91 ++++++++++++++++ examples/Gr3/Shapiro/coastal_0.5_2.cpp.reg | 49 +++++++++ .../Shapiro/gen_shapiro.py} | 13 ++- examples/Gr3/gen_gr3_constant.py | 14 +++ examples/HYCOM/download_hycom.py | 37 +++++++ examples/HYCOM/gen_bnd.py | 33 ++++++ examples/HYCOM/gen_nudge.py | 37 +++++++ examples/NWM/gen_sourcesink.py | 51 +++++++++ examples/Sflux/gen_sflux_era5.py | 37 +++++++ .../gen_sflux_gfs2.py | 0 .../gen_sflux_hrrr3.py | 0 .../{examples_sflux => Sflux}/link_sflux.py | 0 examples/example_4/test_hycom.py | 21 ---- examples/examples_gridgr3/02_gen_shapiro.py | 15 --- examples/examples_hindcast/01_nwm.py | 32 ------ examples/examples_hindcast/02_era5.py | 18 ---- examples/examples_sflux/gen_sflux_era5.py | 31 ------ 21 files changed, 503 insertions(+), 127 deletions(-) create mode 100644 examples/Gr3/Drag/GoME+0.001.reg create mode 100644 examples/Gr3/Drag/Lake_Charles_0.reg rename examples/{examples_gridgr3/01_gen_drag.py => Gr3/Drag/gen_drag.py} (64%) create mode 100644 examples/Gr3/Shapiro/coastal_0.2.cpp.reg create mode 100644 examples/Gr3/Shapiro/coastal_0.5_1.cpp.reg create mode 100644 examples/Gr3/Shapiro/coastal_0.5_2.cpp.reg rename examples/{examples_gridgr3/03_gen_shapiro2.py => Gr3/Shapiro/gen_shapiro.py} (55%) create mode 100644 examples/Gr3/gen_gr3_constant.py create mode 100644 examples/HYCOM/download_hycom.py create mode 100644 examples/HYCOM/gen_bnd.py create mode 100644 examples/HYCOM/gen_nudge.py create mode 100644 examples/NWM/gen_sourcesink.py create mode 100644 examples/Sflux/gen_sflux_era5.py rename examples/{examples_sflux => Sflux}/gen_sflux_gfs2.py (100%) rename examples/{examples_sflux => Sflux}/gen_sflux_hrrr3.py (100%) rename examples/{examples_sflux => Sflux}/link_sflux.py (100%) delete mode 100644 examples/example_4/test_hycom.py delete mode 100644 examples/examples_gridgr3/02_gen_shapiro.py delete mode 100644 examples/examples_hindcast/01_nwm.py delete mode 100644 examples/examples_hindcast/02_era5.py delete mode 100644 examples/examples_sflux/gen_sflux_era5.py diff --git a/examples/Gr3/Drag/GoME+0.001.reg b/examples/Gr3/Drag/GoME+0.001.reg new file mode 100644 index 00000000..01133ab2 --- /dev/null +++ b/examples/Gr3/Drag/GoME+0.001.reg @@ -0,0 +1,23 @@ +Region written by ACE/gredit +1 +20 1 +-70.711855 41.847215 +-70.588878 41.984546 +-70.567176 42.172472 +-70.473134 42.729025 +-70.205477 43.220526 +-69.474846 43.632519 +-68.563365 43.950549 +-67.196145 44.427594 +-66.819978 44.832360 +-66.805510 44.933551 +-67.268484 45.092566 +-68.353580 45.186529 +-69.735269 45.034742 +-71.804185 43.632519 +-72.079076 43.090423 +-71.898226 42.541098 +-71.138659 42.201384 +-70.899938 42.078509 +-70.827599 42.042369 +-70.762493 41.883354 diff --git a/examples/Gr3/Drag/Lake_Charles_0.reg b/examples/Gr3/Drag/Lake_Charles_0.reg new file mode 100644 index 00000000..607ea60e --- /dev/null +++ b/examples/Gr3/Drag/Lake_Charles_0.reg @@ -0,0 +1,20 @@ +Region written by ACE/gredit +1 +17 1 +-93.248568 30.287045 +-93.189017 30.265876 +-93.200469 30.222967 +-93.231390 30.185778 +-93.270327 30.117695 +-93.296094 30.045607 +-93.316708 30.037597 +-93.324152 29.986106 +-93.323579 29.943769 +-93.317853 29.807030 +-93.323579 29.773275 +-93.363089 29.773275 +-93.363089 29.885412 +-93.348201 29.988395 +-93.341330 30.082795 +-93.350492 30.142869 +-93.321862 30.214957 diff --git a/examples/examples_gridgr3/01_gen_drag.py b/examples/Gr3/Drag/gen_drag.py similarity index 64% rename from examples/examples_gridgr3/01_gen_drag.py rename to examples/Gr3/Drag/gen_drag.py index fdb14016..de9a63e9 100644 --- a/examples/examples_gridgr3/01_gen_drag.py +++ b/examples/Gr3/Drag/gen_drag.py @@ -1,16 +1,18 @@ from pyschism.mesh.hgrid import Hgrid from pyschism.mesh.fgrid import DragCoefficient -hgrid=Hgrid.open('/sciclone/schism10/whuang07/NWM/Case1/RUN07b_ZG/hgrid.gr3', crs='epsg:4326') +hgrid=Hgrid.open('./hgrid.gr3', crs='epsg:4326') depth1=-1.0 depth2=-3.0 bfric_river=0.0025 bfric_land=0.025 fgrid=DragCoefficient.linear_with_depth(hgrid, depth1, depth2, bfric_river, bfric_land) + +#modify values in regions regions=['GoME+0.001.reg', 'Lake_Charles_0.reg'] values=[0.001, 0.0] -flags=[1, 0] +flags=[1, 0] # 0: reset to given value, 1: add given value to the current value for reg, value, flag in zip(regions, values, flags): - fgrid.modify_by_region(hgrid, f'Drag_Python/{reg}', 0.001,depth1, 1) + fgrid.modify_by_region(hgrid, f'./{reg}', value, depth1, flag) fgrid.write('drag.gr3', overwrite=True) diff --git a/examples/Gr3/Shapiro/coastal_0.2.cpp.reg b/examples/Gr3/Shapiro/coastal_0.2.cpp.reg new file mode 100644 index 00000000..a5b996a1 --- /dev/null +++ b/examples/Gr3/Shapiro/coastal_0.2.cpp.reg @@ -0,0 +1,100 @@ + +1 +97 1 +-2087974.19892804 2735461.69114659 +-2046267.92519898 2858522.40115863 +-2035841.34196965 2913215.91823516 +-2050438.86055683 3011208.43213511 +-2027500.33426503 3093248.84537155 +-1967025.86885042 3141105.86693440 +-1848162.52744231 3200357.41374636 +-1754323.30236814 3264166.71804448 +-1664654.96156731 3280119.07046258 +-1583328.03311239 3245935.65378028 +-1479062.77871217 3216310.14785139 +-1443612.61986214 3191242.29593169 +-1395650.62825981 3182126.70366696 +-1328920.94616267 3223146.84159508 +-1295556.11362641 3248214.71207795 +-1245508.97358513 3227704.65888461 +-1220485.30601345 3236820.15783981 +-1189205.77591216 3286955.88916706 +-1157926.34976052 3327976.01260565 +-1095367.30002346 3327976.01664844 +-1034893.64740146 3341649.43003119 +-1009869.99155960 3359880.53370987 +-957737.49618650 3364438.33868968 +-916031.54254958 3357601.68995285 +-870154.91374796 3316581.57039991 +-863899.05146884 3289234.86073454 +-813851.84959525 3257330.23677857 +-767975.22361845 3259609.19461440 +-724183.98731015 3225425.72426741 +-659539.67575409 3143385.59543170 +-605321.90278038 3029440.84194120 +-544848.26862581 2865360.36119174 +-536507.01806169 2774204.62507747 +-567786.52933550 2737742.30829268 +-513568.75721901 2728626.70139270 +-413474.36214879 2749136.76160622 +-327977.07926365 2790156.88185876 +-313380.06901297 2835734.69426225 +-311294.68014425 2906380.48159590 +-309209.39297320 2999815.17512957 +-323806.50493764 3100086.51825603 +-353000.72887166 3207194.62229566 +-394706.67592320 3284677.05921303 +-419730.32554551 3398621.70152670 +-388450.81439184 3503450.95949357 +-286271.13232808 3606001.14873164 +-213285.67427848 3665252.48296850 +-123617.81714054 3713109.25278223 +-73570.61978679 3724503.70588242 +-42291.10873009 3715388.09887429 +-38120.53437502 3767802.67219048 +-25608.70961326 3792870.53580273 +20267.81168508 3817938.39941500 +49462.03556304 3820217.24550671 +76570.97226339 3820217.24550672 +78656.25944087 3854400.71612716 +143300.46703233 3904536.44335163 +176665.26526503 3929604.19564317 +166238.82937764 3995692.17947167 +128703.45678998 4095963.52259992 +168324.11655511 4189398.21613237 +205860.12695341 4280554.58522318 +235053.71302052 4330689.67947690 +272589.08560818 4364873.15009736 +312209.74537334 4430961.13392587 +322636.28295726 4494770.16034198 +360171.65554500 4510722.41694602 +504057.28436660 4549463.57974983 +593725.14149262 4560858.03284954 +718843.08408122 4565415.83635253 +760549.13276806 4576810.40077240 +752207.88234603 4640619.42718910 +731354.90884917 4679360.58999345 +687563.57299579 4738611.92422672 +702160.68494185 4800141.99323051 +773060.85582705 4845719.91695078 +873155.25071970 4889018.88325953 +987846.65615984 4934596.80698683 +1067088.07797614 4996126.98733013 +1213058.99619721 5041704.91115431 +1206803.13460117 5025752.65453882 +1081685.19017228 4945991.26010490 +1073343.93974961 4900413.44769795 +1104623.45139774 4829767.66036284 +1146329.39920358 4804699.79675655 +1225570.82246726 4836604.31001547 +1315238.68271874 4909528.94360159 +1375712.42082002 4934596.80740049 +1525853.92409735 4977895.88606557 +1659313.13598129 5025752.65830226 +1782345.81095815 5110071.75004494 +1402821.25358723 5155649.66541496 +356001.08119003 4957385.82450028 +-1476977.37949447 4209908.26326187 +-2146364.75561788 3555864.56100148 +-2219352.08710189 3104642.21781792 +-2198496.24906293 2544033.50007754 diff --git a/examples/Gr3/Shapiro/coastal_0.5_1.cpp.reg b/examples/Gr3/Shapiro/coastal_0.5_1.cpp.reg new file mode 100644 index 00000000..9afe01cb --- /dev/null +++ b/examples/Gr3/Shapiro/coastal_0.5_1.cpp.reg @@ -0,0 +1,91 @@ + +1 +88 1 +-2110399.18520885 2759473.72383486 +-2081889.68589075 2785073.56786075 +-2063382.21061246 2828875.00377477 +-2045376.06793056 2860259.28829079 +-2041631.15575467 2889859.31897550 +-2050273.45831639 2944864.84574453 +-2066099.21736074 2999260.29785134 +-2059761.59013415 3050835.65094407 +-2045204.47150868 3096399.91130284 +-2008077.08435444 3136493.72825179 +-1978150.61985263 3155831.32141113 +-1921370.89089629 3186179.04642157 +-1863396.84021269 3212362.60420253 +-1814267.94187892 3257106.74178290 +-1768008.39598740 3283788.44727055 +-1716889.84195376 3301314.99811018 +-1678640.52151209 3309575.78697443 +-1637781.41386630 3311793.30771635 +-1571558.69822787 3290215.93267294 +-1540943.09169962 3287693.05928811 +-1528630.61329063 3293718.58695410 +-1518398.13790829 3303665.52616034 +-1490368.31685577 3299219.50771532 +-1474697.10182308 3283786.10215843 +-1441844.31707609 3251616.99261998 +-1401353.14921681 3241526.64277997 +-1373217.96164349 3250890.79753969 +-1350856.85789183 3254498.93996371 +-1329013.18087201 3245960.57331146 +-1314354.47629823 3268355.32456693 +-1299741.32237654 3275148.86019249 +-1256432.42564340 3244295.82893895 +-1234975.04701564 3238344.11763361 +-1224678.33033147 3244528.12009131 +-1240677.70001922 3265343.46807355 +-1268144.46963296 3282599.09052713 +-1273006.93803143 3291224.62897666 +-1252962.16745369 3313543.91366011 +-1247674.62477571 3324124.55235404 +-1253138.00288865 3342626.38854693 +-1267743.33668832 3345798.48187961 +-1280910.61764761 3338486.47622982 +-1293391.19487411 3353238.48602371 +-1260905.06140930 3355870.47426355 +-1244134.21608518 3369999.57121295 +-1213750.08526171 3377935.87525777 +-1170399.47392774 3373085.21418655 +-1121329.95529754 3376528.37073100 +-1114086.61240031 3402970.49031818 +-1101245.90130761 3378039.99759169 +-1086464.84718525 3367543.66329230 +-1030836.07351491 3373442.46364716 +-978429.58011614 3382330.04899383 +-920653.90342542 3370067.21322368 +-864055.36757060 3345945.81026967 +-836220.23518445 3306187.48786777 +-808123.82362033 3302131.02595357 +-744705.32699062 3323022.98375028 +-705832.26634758 3334614.91226253 +-663671.58049660 3315495.13165063 +-630645.89563965 3266171.57453546 +-584275.61149532 3234802.93564763 +-572941.16499378 3188870.08573578 +-586700.64540248 3128440.66283929 +-584973.56490607 3091235.73558275 +-560978.75302966 3026675.09217801 +-527341.88613622 2953960.51839202 +-514483.51927664 2933130.78254744 +-496156.57358228 2931248.33890320 +-472456.63392595 2868545.87998493 +-446799.36500404 2863355.00193600 +-422835.81610921 2817054.59683932 +-408438.40169585 2786797.55494828 +-361811.92086997 2793852.43468602 +-328531.34296750 2811030.89660020 +-318210.00837721 2843101.23329873 +-306038.64532001 2976010.86865396 +-316986.67249013 3027394.31296513 +-352790.65305549 3133808.23023285 +-352440.33273801 3173393.35272890 +-400717.80498957 3255997.68513505 +-428030.67124798 3327950.82494407 +-446629.96550044 3411786.06908229 +-436912.32113368 3464794.02532570 +-403951.72282368 3541328.20505098 +-348470.58474651 3597329.59237861 +-390979.85246274 4257886.12245173 +-3112804.85869375 3441796.78992251 diff --git a/examples/Gr3/Shapiro/coastal_0.5_2.cpp.reg b/examples/Gr3/Shapiro/coastal_0.5_2.cpp.reg new file mode 100644 index 00000000..94464e2e --- /dev/null +++ b/examples/Gr3/Shapiro/coastal_0.5_2.cpp.reg @@ -0,0 +1,49 @@ + +1 +46 1 +-347729.54378964 3592328.90283183 +-287049.50749085 3633932.26631487 +-238940.60829766 3671166.00903764 +-220603.01641726 3680502.99873341 +-213230.99489887 3690855.45476839 +-207335.79283072 3717710.72795215 +-184848.75778042 3743707.79989715 +-161089.81352332 3760345.03390648 +-123576.63259316 3770143.90830229 +-96616.52158290 3766577.74763832 +-82464.16923896 3768389.45051966 +-85053.83281319 3775526.16922139 +-70136.31755149 3805436.46811014 +-42335.13910624 3832741.72637059 +-5402.84839885 3851514.81537110 +37485.50162728 3858735.17772310 +47975.08821008 3851685.63446105 +58366.43239770 3849668.35407018 +65331.74771613 3862728.97462225 +102467.88175704 3893406.54560510 +143119.78886254 3916779.50296324 +163160.45588122 3921204.80916108 +169367.14931525 3968781.87274773 +154625.28775288 4002578.92647125 +125597.78017109 4071551.94429403 +114612.51529213 4104596.88269993 +123386.81656259 4128439.25114752 +146962.08637912 4163893.71614900 +176439.25891736 4223803.51916991 +202178.90905483 4264183.69741783 +249443.33387059 4359889.60984369 +293867.24060720 4401399.22129333 +314166.55820454 4456659.46172878 +319013.84200704 4496205.91718178 +325248.06109014 4509818.55177440 +362470.87894428 4516554.51944404 +451208.41200986 4537963.22321620 +531323.70942565 4568968.81645830 +676680.41556462 4599593.57973840 +725380.40754349 4634044.56243357 +735184.97725482 4660183.65325188 +715142.18969791 4687771.05611060 +668300.75587200 4755359.15678493 +158276.73532632 4703201.81839697 +-386028.18928844 4529800.38140454 +-576207.03719000 3910267.46264677 diff --git a/examples/examples_gridgr3/03_gen_shapiro2.py b/examples/Gr3/Shapiro/gen_shapiro.py similarity index 55% rename from examples/examples_gridgr3/03_gen_shapiro2.py rename to examples/Gr3/Shapiro/gen_shapiro.py index ac016311..22a4185b 100644 --- a/examples/examples_gridgr3/03_gen_shapiro2.py +++ b/examples/Gr3/Shapiro/gen_shapiro.py @@ -2,22 +2,21 @@ from pyschism.mesh.gridgr3 import Shapiro if __name__ == '__main__': - hgrid = Hgrid.open('/sciclone/schism10/whuang07/NWM/Case1/RUN07b_ZG/hgrid.gr3', crs='EPSG:4326') + hgrid = Hgrid.open('./hgrid.gr3', crs='EPSG:4326') shapiro_max = 0.5 threshold_slope = 0.5 depths = [-99999, 20, 50] # tweaks in shallow waters shapiro_vals1 = [0.2, 0.2, 0.05] # tweaks in shallow waters - regdir = '/sciclone/data10/lcui01/schism/src/Utility/Pre-Processing/NWM/Shapiro_python' regions = [ - f'{regdir}/coastal_0.2.cpp.reg', - f'{regdir}/coastal_0.5_1.cpp.reg', - f'{regdir}/coastal_0.5_2.cpp.reg' + f'./coastal_0.2.cpp.reg', + f'./coastal_0.5_1.cpp.reg', + f'./coastal_0.5_2.cpp.reg' ] # tweaks in regions, the order matters shapiro_vals2 = [0.2, 0.5, 0.5] # tweaks in regions, the order matters i_set_add_s = [0, 0, 0] shapiro = Shapiro.slope_filter( hgrid, shapiro_vals1, depths, shapiro_max, threshold_slope, - regions, shapiro_vals2, i_set_add_s, lonc=77.07, latc=24.0) - shapiro.write('shapiro_pyschism.gr3', overwrite=True) + regions, shapiro_vals2, i_set_add_s, lonc=-77.07, latc=24.0) + shapiro.write('shapiro.gr3', overwrite=True) diff --git a/examples/Gr3/gen_gr3_constant.py b/examples/Gr3/gen_gr3_constant.py new file mode 100644 index 00000000..2ce3ad57 --- /dev/null +++ b/examples/Gr3/gen_gr3_constant.py @@ -0,0 +1,14 @@ +from pyschism.mesh.hgrid import Gr3 + +if __name__ == '__main__': + hgrid = Gr3.open('./hgrid.gr3', crs='epsg:4326') + grd = hgrid.copy() + + gr3_names = ['albedo', 'diffmax', 'diffmin', 'watertype', 'windrot_geo2proj'] + values = [0.1, 1.0, 1e-6, 1.0, 0.0] + + for name, value in zip(gr3_names, values): + grd.description = name + grd.nodes.values[:] = value + grd.write(f'{name}.gr3', overwrite=True) + diff --git a/examples/HYCOM/download_hycom.py b/examples/HYCOM/download_hycom.py new file mode 100644 index 00000000..0c0e5732 --- /dev/null +++ b/examples/HYCOM/download_hycom.py @@ -0,0 +1,37 @@ +from datetime import datetime +import logging +from time import time + +from matplotlib.transforms import Bbox + +from pyschism.forcing.hycom.hycom2schism import DownloadHycom +from pyschism.mesh.hgrid import Hgrid + +''' +Download hycom data for Fortran scripts. +Default is to download data for generating initial condition (use hgrid + as parameter to cover the whole domain). +Optionally, download data for generating th.nc (bnd=True) and nu.nc (nudge=True) + (use bbox as parameter to cut small domain). +''' +logging.basicConfig( + format="[%(asctime)s] %(name)s %(levelname)s: %(message)s", + force=True, +) +logger = logging.getLogger('pyschism') +logger.setLevel(logging.INFO) + +if __name__ == '__main__': + hgrid = Hgrid.open('./hgrid.gr3', crs='epsg:4326') + date=datetime(2018, 1, 1) + + #xmin = -65 + #xmax = -50 + #ymin = 0 + #ymax = 53 + #bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) + hycom = DownloadHycom(hgrid=hgrid) + t0 = time() + hycom.fetch_data(date, rnday=1, bnd=False, nudge=False, outdir='./') + print(f'It took {(time()-t0)/60} mins to download') + diff --git a/examples/HYCOM/gen_bnd.py b/examples/HYCOM/gen_bnd.py new file mode 100644 index 00000000..681eb6c1 --- /dev/null +++ b/examples/HYCOM/gen_bnd.py @@ -0,0 +1,33 @@ +from datetime import datetime +import logging + +from pyschism.mesh.hgrid import Hgrid +from pyschism.forcing.hycom.hycom2schism import OpenBoundaryInventory + +''' +outputs: + elev2D.th.nc (elev=True) + SAL_3D.th.nc (TS=True) + TEM_3D.th.nc (TS=True) + uv3D.th.nc (UV=True) +''' + +logging.basicConfig( + format="[%(asctime)s] %(name)s %(levelname)s: %(message)s", + force=True, +) +logger = logging.getLogger('pyschism') +logger.setLevel(logging.INFO) + +if __name__ == "__main__": + hgrid=Hgrid.open('./hgrid.gr3', crs='epsg:4326') + vgrid='./vgrid.in' + outdir='./' + start_date=datetime(2018, 8, 1) + rnday=20 + + #boundary + bnd=OpenBoundaryInventory(hgrid, vgrid) + + #ocean_bnd_ids - segment indices, starting from zero. + bnd.fetch_data(outdir, start_date, rnday, elev2D=True, TS=True, UV=True, ocean_bnd_ids=[0,1,2]) diff --git a/examples/HYCOM/gen_nudge.py b/examples/HYCOM/gen_nudge.py new file mode 100644 index 00000000..a430ea42 --- /dev/null +++ b/examples/HYCOM/gen_nudge.py @@ -0,0 +1,37 @@ +from datetime import datetime, timedelta +import logging + +from pyschism.mesh import Hgrid, Vgrid +from pyschism.forcing.hycom.hycom2schism import Nudge + +''' +outputs: + SAL_nudge.gr3 + SAL_nu.nc + TEM_nudge.gr3 + TEM_nu.nc +''' + +logging.basicConfig( + format="[%(asctime)s] %(name)s %(levelname)s: %(message)s", + force=True, +) +logger = logging.getLogger('pyschism') +logger.setLevel(logging.INFO) + +if __name__ == '__main__': + + + hgrid=Hgrid.open('./hgrid.gr3', crs='epsg:4326') + vgrid='./vgrid.in' + outdir='./' + start_date = datetime(2018, 1, 1) + rnday=30 + + #ocean_bnd_ids - segment indices, starting from zero. + nudge=Nudge(hgrid=hgrid, ocean_bnd_ids=[0, 1, 9]) + + #rlmax - max relax distance in m or degree + #rnu_day - max relax strength in days + #restart = True will append to the existing nc file, works when first try doesn't break. + nudge.fetch_data(outdir, vgrid, start_date, rnday, restart=False, rlmax=1, rnu_day=1) diff --git a/examples/NWM/gen_sourcesink.py b/examples/NWM/gen_sourcesink.py new file mode 100644 index 00000000..dd83a4dc --- /dev/null +++ b/examples/NWM/gen_sourcesink.py @@ -0,0 +1,51 @@ +from datetime import datetime +from time import time +import pathlib +import logging + +from pyschism.forcing.source_sink.nwm import NationalWaterModel, NWMElementPairings +from pyschism.mesh import Hgrid + +logging.basicConfig( + format="[%(asctime)s] %(name)s %(levelname)s: %(message)s", + force=True, +) +logging.captureWarnings(True) + +log_level = logging.DEBUG +logging.getLogger('pyschism').setLevel(log_level) + +if __name__ == '__main__': + + startdate = datetime(2017, 12, 1) + rnday = 396 + hgrid = Hgrid.open("./hgrid.gr3", crs="epsg:4326") + + t0 = time() + + #source/sink json files, if not exist, it will call NWMElementPairings to generate. + sources_pairings = pathlib.Path('./sources.json') + sinks_pairings = pathlib.Path('./sinks.json') + output_directory = pathlib.Path('./') + + #input directory which saves downloaded nc files + #cache = pathlib.Path(f'./{startdate.strftime("%Y%m%d")}') + #cache.mkdir(exist_ok=True, parents=True) + cache = pathlib.Path('./NWM_v2.1') + + # check if source/sink json file exists + if all([sources_pairings.is_file(), sinks_pairings.is_file()]) is False: + pairings = NWMElementPairings(hgrid) + sources_pairings.parent.mkdir(exist_ok=True, parents=True) + pairings.save_json(sources=sources_pairings, sinks=sinks_pairings) + else: + pairings = NWMElementPairings.load_json( + hgrid, + sources_pairings, + sinks_pairings) + + #check nc files, if not exist will download + nwm=NationalWaterModel(pairings=pairings, cache=cache) + + nwm.write(output_directory, hgrid, startdate, rnday, overwrite=True) + print(f'It took {time()-t0} seconds to generate source/sink') diff --git a/examples/Sflux/gen_sflux_era5.py b/examples/Sflux/gen_sflux_era5.py new file mode 100644 index 00000000..6b1c9b9f --- /dev/null +++ b/examples/Sflux/gen_sflux_era5.py @@ -0,0 +1,37 @@ +from datetime import datetime +from time import time +import pathlib + +from pyschism.forcing.nws.nws2.era5 import ERA5 + +from pyschism.mesh.hgrid import Hgrid + +if __name__ == '__main__': + + startdate=datetime(2018, 1, 1) + rnday=366 + + t0=time() + hgrid=Hgrid.open('./hgrid.gr3',crs='EPSG:4326') + bbox = hgrid.get_bbox('EPSG:4326', output_type='bbox') + + er=ERA5() + outdir = pathlib.Path('./') + interval = 1 + er.write( + outdir=outdir, + start_date=startdate, + rnday=rnday, + air=True, #sflux_air_1 + rad=True, #sflux_rad_1 + prc=True, #sflux_prc_1 + bbox=bbox, + output_interval=interval, #raw data is hourly + overwrite=True, + tmpdir=outdir) #default tmpdir is system's tmp (/tmp), it is possbile /tmp has no enough space for large dataset. + + #write sflux_inputs.txt + with open("./sflux_inputs.txt", "w") as f: + f.write("&sflux_inputs\n/\n") + + print(f'It took {(time()-t0)/60} minutes to generate {rnday} days') diff --git a/examples/examples_sflux/gen_sflux_gfs2.py b/examples/Sflux/gen_sflux_gfs2.py similarity index 100% rename from examples/examples_sflux/gen_sflux_gfs2.py rename to examples/Sflux/gen_sflux_gfs2.py diff --git a/examples/examples_sflux/gen_sflux_hrrr3.py b/examples/Sflux/gen_sflux_hrrr3.py similarity index 100% rename from examples/examples_sflux/gen_sflux_hrrr3.py rename to examples/Sflux/gen_sflux_hrrr3.py diff --git a/examples/examples_sflux/link_sflux.py b/examples/Sflux/link_sflux.py similarity index 100% rename from examples/examples_sflux/link_sflux.py rename to examples/Sflux/link_sflux.py diff --git a/examples/example_4/test_hycom.py b/examples/example_4/test_hycom.py deleted file mode 100644 index b9aa1949..00000000 --- a/examples/example_4/test_hycom.py +++ /dev/null @@ -1,21 +0,0 @@ -import numpy as np -from datetime import datetime, timedelta -import logging -from time import time - -from netCDF4 import Dataset -from matplotlib.transforms import Bbox - -from pyschism.mesh import Hgrid -from pyschism.forcing.hycom.hycom import HotStartInventory - -logger = logging.getLogger(__name__) - -if __name__ == '__main__': - logger.info('Starting!') - start_date = datetime.now() - hgrid = Hgrid.open('../ICOGS3D_dep/hgrid.gr3',crs='EPSG:4326') - hycom = HotStartInventory() - t0 = time() - hycom.fetch_data(hgrid, start_date) - print(f'Took {time()-t0} seconds') diff --git a/examples/examples_gridgr3/02_gen_shapiro.py b/examples/examples_gridgr3/02_gen_shapiro.py deleted file mode 100644 index f8d41516..00000000 --- a/examples/examples_gridgr3/02_gen_shapiro.py +++ /dev/null @@ -1,15 +0,0 @@ -from pyschism.mesh.hgrid import Hgrid -from pyschism.mesh.base import Gr3 -from pyschism.mesh.gridgr3 import Shapiro - -hgrid=Hgrid.open('/sciclone/schism10/whuang07/NWM/Case1/RUN07b_ZG/hgrid.gr3', crs='epsg:4326') -shapiro=Shapiro.from_binary(hgrid,lonc = -77.07, latc = 24.0) -regions=['coastal_0.2.cpp.reg', 'coastal_0.5_1.cpp.reg', 'coastal_0.5_2.cpp.reg'] -values=[0.2, 0.5, 0.5] -flags=[0, 0, 0] -depth1=-99999 -#hgrid.write('hgrid_after_cpp.gr3', overwrite=True) - -for reg, value, flag in zip(regions, values, flags): - shapiro.modify_by_region(hgrid, f'Shapiro_python/{reg}', value, depth1, flag) -shapiro.write('shapiro_region.gr3', overwrite=True) diff --git a/examples/examples_hindcast/01_nwm.py b/examples/examples_hindcast/01_nwm.py deleted file mode 100644 index f5b25aba..00000000 --- a/examples/examples_hindcast/01_nwm.py +++ /dev/null @@ -1,32 +0,0 @@ -from datetime import datetime -from time import time - -import logging - -from pyschism.forcing import source_sink -from pyschism.mesh import Hgrid - - -logging.basicConfig( - format="[%(asctime)s] %(name)s %(levelname)s: %(message)s", - force=True, -) -logging.captureWarnings(True) - -log_level = logging.DEBUG -source_sink.nwm.logger.setLevel(log_level) -source_sink.base.logger.setLevel(log_level) - -startdate=datetime(2018, 8, 24) -rnday=40 -hgrid = Hgrid.open("./hgrid.gr3", crs="epsg:4326") - -t0=time() -nwm=source_sink.nwm.NationalWaterModel() -nwm.write('./', hgrid, startdate, rnday, overwrite=True) -print(f'It took {time()-t0} seconds to generate source/sink') -nwm.pairings.save_json(sources=True, sinks=True) -nwm.pairings.sources_gdf.to_file("sources_run07b_new.shp") -nwm.pairings.sinks_gdf.to_file("sinks_run07b_new.shp") -nwm.pairings.intersection.to_file("intersection_run07b_new.shp") -nwm.pairings.reaches.to_file("reaches_run07b_new.shp") diff --git a/examples/examples_hindcast/02_era5.py b/examples/examples_hindcast/02_era5.py deleted file mode 100644 index 8fecf832..00000000 --- a/examples/examples_hindcast/02_era5.py +++ /dev/null @@ -1,18 +0,0 @@ -from datetime import datetime -from time import time -import pathlib - -from pyschism.forcing.nws.nws2.era5 import ERA5 - -from pyschism.mesh.hgrid import Hgrid - -startdate=datetime(2018,7,1) - -t0=time() -hgrid=Hgrid.open('hgrid.gr3',crs='EPSG:4326') -bbox = hgrid.get_bbox('EPSG:4326', output_type='bbox') - -er=ERA5() -outdir = pathlib.Path('./ERA5') -er.write(outdir=outdir, start_date=startdate, rnday=10, air=True, rad=True, prc=True, bbox=bbox, overwrite=True) -print(f'It took {(time()-t0)/60} minutes to generate 10 days') diff --git a/examples/examples_sflux/gen_sflux_era5.py b/examples/examples_sflux/gen_sflux_era5.py deleted file mode 100644 index 4f127dc9..00000000 --- a/examples/examples_sflux/gen_sflux_era5.py +++ /dev/null @@ -1,31 +0,0 @@ -from datetime import datetime -from time import time -import pathlib -import logging - -from pyschism.forcing.nws.nws2.era5 import ERA5 - -from pyschism.mesh.hgrid import Hgrid - -logging.basicConfig( - format="[%(asctime)s] %(name)s %(levelname)s: %(message)s", - force=True, -) -logging.captureWarnings(True) - -log_level = logging.DEBUG -logging.getLogger('pyschism').setLevel(log_level) - -if __name__ == '__main__': - - startdate=datetime(2017, 4, 26) - rnday=18 - - t0=time() - hgrid=Hgrid.open('./hgrid.gr3',crs='EPSG:4326') - bbox = hgrid.get_bbox('EPSG:4326', output_type='bbox') - - er=ERA5() - outdir = pathlib.Path('./') - er.write(outdir=outdir, start_date=startdate, rnday=rnday, air=True, rad=True, prc=True, bbox=bbox, overwrite=True) - print(f'It took {(time()-t0)/60} minutes to generate {rnday} days') From 7dac44ab4cd4c523d174a45d809cb0134fdd0b57 Mon Sep 17 00:00:00 2001 From: cuill Date: Wed, 15 Mar 2023 14:24:39 -0400 Subject: [PATCH 25/31] Added an example script for bctides --- examples/Bctides/gen_bctides.py | 70 +++++++++++++++++++++++++++++++++ 1 file changed, 70 insertions(+) create mode 100644 examples/Bctides/gen_bctides.py diff --git a/examples/Bctides/gen_bctides.py b/examples/Bctides/gen_bctides.py new file mode 100644 index 00000000..39097a3c --- /dev/null +++ b/examples/Bctides/gen_bctides.py @@ -0,0 +1,70 @@ +from time import time +import os +from datetime import datetime, timedelta +import logging +import pathlib +import f90nml + +from pyschism.mesh import Hgrid +from pyschism.forcing.bctides import Bctides, iettype, ifltype, isatype, itetype +from pyschism.forcing.source_sink.nwm import NationalWaterModel, NWMElementPairings +from pyschism.forcing.nws import NWS2, GFS, HRRR + +if __name__ == "__main__": + + ''' + Assume files are already located in: + database='fes2014' + ~/.local/share/fes2014/eastward_velocity/ + ~/.local/share/fes2014/northward_velocity/ + ~/.local/share/fes2014/ocean_tide_extrapolated/ + database = 'tpxo' + ~/.local/share/tpxo/ + ''' + #setup logging + logging.basicConfig( + format = "[%(asctime)s] %(name)s %(levelname)s: %(message)s", + force=True, + ) + logging.getLogger("pyschism").setLevel(logging.DEBUG) + + start_date = datetime(2017, 12, 1) + rnday = 396 + end_date = start_date + timedelta(days=61) + outdir = './' + hgrid = Hgrid.open("./hgrid.gr3", crs="epsg:4326") + + #elevation + iet3 = iettype.Iettype3(constituents='major', database='fes2014') + iet4 = iettype.Iettype4() + iet5 = iettype.Iettype5(iettype3=iet3, iettype4=iet4) + + #velocity + ifl1 = ifltype.Ifltype1() + ifl3 = ifltype.Ifltype3(constituents='major', database='fes2014') + ifl4 = ifltype.Ifltype4() + ifl5 = ifltype.Ifltype5(ifltype3=ifl3, ifltype4=ifl4) + + #salinity + isa4 = isatype.Isatype4() + isa2 = isatype.Isatype2(5, 1) + + #temperature + ite4 = itetype.Itetype4() + ite2 = itetype.Itetype2(0, 1) + + bctides=Bctides( + hgrid, + iettype={'1': iet5, '2': iet5}, + ifltype={'1': ifl5, '2': ifl5, '3': ifl1}, + isatype={'1': isa4, '2':isa4, '3': isa2}, + itetype={'1': ite4, '2':ite4, '3': ite2} + ) + + bctides.write( + outdir, + start_date=start_date, + end_date=end_date, + bctides=True, + overwrite=True + ) From 7ef1d94f5a2e6fdf0aa39ee37f083cddaf7813a8 Mon Sep 17 00:00:00 2001 From: cuill Date: Thu, 16 Mar 2023 09:55:15 -0400 Subject: [PATCH 26/31] fixed negative nudge coeff --- pyschism/forcing/hycom/hycom2schism.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pyschism/forcing/hycom/hycom2schism.py b/pyschism/forcing/hycom/hycom2schism.py index 5c30e6cb..05b100ec 100644 --- a/pyschism/forcing/hycom/hycom2schism.py +++ b/pyschism/forcing/hycom/hycom2schism.py @@ -513,8 +513,7 @@ def gen_nudge(self, outdir: Union[str, os.PathLike], rlmax = 1.5, rnu_day=0.25): out[out<0] = 0 out[out>rnu_max] = rnu_max fp = out>0 - - nudge_coeff[fp] = -out[fp] + nudge_coeff[fp] = out[fp] idxs_nudge=np.zeros(NP, dtype=int) idxs=np.where(out > 0)[0] From d07c11bf4b0349c52b96dc1a6db5f4b832a84e7f Mon Sep 17 00:00:00 2001 From: cuill Date: Thu, 16 Mar 2023 17:10:12 -0400 Subject: [PATCH 27/31] fixed time_index issue for restarting (bnd and nudge) --- pyschism/forcing/hycom/hycom2schism.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/pyschism/forcing/hycom/hycom2schism.py b/pyschism/forcing/hycom/hycom2schism.py index 05b100ec..941fa9fa 100644 --- a/pyschism/forcing/hycom/hycom2schism.py +++ b/pyschism/forcing/hycom/hycom2schism.py @@ -11,7 +11,7 @@ import numpy as np import scipy as sp -from numba import jit, prange +#from numba import jit, prange import netCDF4 as nc from netCDF4 import Dataset from matplotlib.transforms import Bbox @@ -345,8 +345,9 @@ def fetch_data(self, outdir: Union[str, os.PathLike], start_date, rnday, ocean_b timevector = self.timevector it0 = 0 elif restart: - timevector = self.timevector[time_idx_restart+1:] - it0 = time_idx_restart+1 + #restart from one day earlier + timevector = self.timevector[time_idx_restart-1:] + it0 = time_idx_restart-1 for it1, date in enumerate(timevector): @@ -651,8 +652,9 @@ def fetch_data(self, outdir: Union[str, os.PathLike], vgrid, start_date, rnday, logger.info('**** Accessing GOFS data*****') if restart: - timevector = self.timevector[time_idx_restart+1:] - it0 = time_idx_restart+1 + #restart from one day earlier to make sure all files consistant + timevector = self.timevector[time_idx_restart-1:] + it0 = time_idx_restart-1 else: timevector = self.timevector it0 = 0 From fa23ce69b2a88f65f8312194966089a5c0bbed41 Mon Sep 17 00:00:00 2001 From: cuill Date: Wed, 29 Mar 2023 10:04:55 -0400 Subject: [PATCH 28/31] added sub-sampling for downloading hycom --- examples/HYCOM/download_hycom.py | 2 +- pyschism/forcing/hycom/hycom2schism.py | 16 ++++++++-------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/examples/HYCOM/download_hycom.py b/examples/HYCOM/download_hycom.py index 0c0e5732..b97d242f 100644 --- a/examples/HYCOM/download_hycom.py +++ b/examples/HYCOM/download_hycom.py @@ -32,6 +32,6 @@ #bbox = Bbox.from_extents(xmin, ymin, xmax, ymax) hycom = DownloadHycom(hgrid=hgrid) t0 = time() - hycom.fetch_data(date, rnday=1, bnd=False, nudge=False, outdir='./') + hycom.fetch_data(date, rnday=1, bnd=False, nudge=False, sub_sample=3, outdir='./') print(f'It took {(time()-t0)/60} mins to download') diff --git a/pyschism/forcing/hycom/hycom2schism.py b/pyschism/forcing/hycom/hycom2schism.py index 941fa9fa..0fa8e911 100644 --- a/pyschism/forcing/hycom/hycom2schism.py +++ b/pyschism/forcing/hycom/hycom2schism.py @@ -762,7 +762,7 @@ def __init__(self, hgrid=None, bbox=None): elif bbox is not None: self.bbox = bbox - def fetch_data(self, start_date, rnday=1, fmt='schism', bnd=False, nudge=False, outdir=None): + def fetch_data(self, start_date, rnday=1, fmt='schism', bnd=False, nudge=False, sub_sample=1, outdir=None): ''' start_date: datetime.datetime rnday: integer @@ -785,13 +785,13 @@ def fetch_data(self, start_date, rnday=1, fmt='schism', bnd=False, nudge=False, time_idx, lon_idx1, lon_idx2, lat_idx1, lat_idx2, x2, y2, isLonSame = get_idxs(date, database, self.bbox) if fmt == 'schism': - url_ssh = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:1:{lat_idx2}],' + \ - f'lon[{lon_idx1}:1:{lon_idx2}],depth[0:1:-1],time[{time_idx}:8:{time_idx}],' + \ - f'surf_el[{time_idx}][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ - f'water_u[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ - f'water_v[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ - f'water_temp[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}],' + \ - f'salinity[{time_idx}][0:1:39][{lat_idx1}:1:{lat_idx2}][{lon_idx1}:1:{lon_idx2}]' + url_ssh = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:{sub_sample}:{lat_idx2}],' + \ + f'lon[{lon_idx1}:{sub_sample}:{lon_idx2}],depth[0:1:-1],time[{time_idx}:8:{time_idx}],' + \ + f'surf_el[{time_idx}][{lat_idx1}:{sub_sample}:{lat_idx2}][{lon_idx1}:{sub_sample}:{lon_idx2}],' + \ + f'water_u[{time_idx}][0:1:39][{lat_idx1}:{sub_sample}:{lat_idx2}][{lon_idx1}:{sub_sample}:{lon_idx2}],' + \ + f'water_v[{time_idx}][0:1:39][{lat_idx1}:{sub_sample}:{lat_idx2}][{lon_idx1}:{sub_sample}:{lon_idx2}],' + \ + f'water_temp[{time_idx}][0:1:39][{lat_idx1}:{sub_sample}:{lat_idx2}][{lon_idx1}:{sub_sample}:{lon_idx2}],' + \ + f'salinity[{time_idx}][0:1:39][{lat_idx1}:{sub_sample}:{lat_idx2}][{lon_idx1}:{sub_sample}:{lon_idx2}]' foutname = f'hycom_{date.strftime("%Y%m%d")}.nc' #foutname = f'TS_{i+1}.nc' From 503c130af2fb63e75ba17ca1808f029ef547bc47 Mon Sep 17 00:00:00 2001 From: cuill Date: Fri, 31 Mar 2023 12:41:31 -0400 Subject: [PATCH 29/31] fixed convert_longitude --- pyschism/forcing/hycom/hycom2schism.py | 37 +++++++++++++++++----- pyschism/forcing/nws/nws2/era5.py | 43 ++++++++++++++------------ 2 files changed, 52 insertions(+), 28 deletions(-) diff --git a/pyschism/forcing/hycom/hycom2schism.py b/pyschism/forcing/hycom/hycom2schism.py index 0fa8e911..fbb2eb25 100644 --- a/pyschism/forcing/hycom/hycom2schism.py +++ b/pyschism/forcing/hycom/hycom2schism.py @@ -23,18 +23,35 @@ logger = logging.getLogger(__name__) -def convert_longitude(ds): - lon_name = 'lon' - ds['_lon_adjusted'] = xr.where( - ds[lon_name] > 180, - ds[lon_name] - 360, - ds[lon_name]) +def convert_longitude(ds, bbox): +#https://stackoverflow.com/questions/53345442/about-changing-longitude-array-from-0-360-to-180-to-180-with-python-xarray +#Light_B's solution didn't generate the correct result +#Michael's solution works, but it takes significantly longer to write nc file (~30 mins compared with 5 mins) +#TODO: figure out why it takes much longer with the second method + #lon_attr = ds.coords['lon'].attrs + if bbox.xmin < 0: + logger.info(f'Convert HYCOM longitude from [0, 360) to [-180, 180):') + #ds.coords['lon'] = (ds.coords['lon'] + 180) % 360 - 180 + ds['_lon_adjusted'] = xr.where(ds['lon'] > 180, ds['lon'] - 360, ds['lon']) + elif bbox.xmin > 0: + logger.info(f'Convert HYCOM longitude from [-180, 180) to [0, 360): ') + #ds.coords['lon'] = (ds.coords['lon'] + 360) % 360 - 180 + ds['_lon_adjusted'] = xr.where(ds['lon'] < 0, ds['lon'] + 360, ds['lon']) + + t0 = time() ds = ( ds.swap_dims({lon_name: '_lon_adjusted'}) .sel(**{'_lon_adjusted': sorted(ds._lon_adjusted)}) .drop(lon_name) ) ds = ds.rename({'_lon_adjusted': lon_name}) + #ds = ds.sortby(ds.lon) + #ds.coords['lon'].attrs = lon_attr + logger.info(f'swap dims took {time()-t0} seconds!') + + #make sure it is clipped to the bbox + ds = ds.sel(lon=slice(bbox.xmin - 0.5, bbox.xmax + 0.5)) + return ds def get_database(date, Bbox=None): @@ -786,7 +803,7 @@ def fetch_data(self, start_date, rnday=1, fmt='schism', bnd=False, nudge=False, if fmt == 'schism': url_ssh = f'https://tds.hycom.org/thredds/dodsC/{database}?lat[{lat_idx1}:{sub_sample}:{lat_idx2}],' + \ - f'lon[{lon_idx1}:{sub_sample}:{lon_idx2}],depth[0:1:-1],time[{time_idx}:8:{time_idx}],' + \ + f'lon[{lon_idx1}:{sub_sample}:{lon_idx2}],depth[0:1:-1],time[{time_idx}],' + \ f'surf_el[{time_idx}][{lat_idx1}:{sub_sample}:{lat_idx2}][{lon_idx1}:{sub_sample}:{lon_idx2}],' + \ f'water_u[{time_idx}][0:1:39][{lat_idx1}:{sub_sample}:{lat_idx2}][{lon_idx1}:{sub_sample}:{lon_idx2}],' + \ f'water_v[{time_idx}][0:1:39][{lat_idx1}:{sub_sample}:{lat_idx2}][{lon_idx1}:{sub_sample}:{lon_idx2}],' + \ @@ -815,14 +832,18 @@ def fetch_data(self, start_date, rnday=1, fmt='schism', bnd=False, nudge=False, if not isLonSame: logger.info('Lon is not the same!') - ds = convert_longitude(ds) + ds = convert_longitude(ds, self.bbox) ds = ds.rename_dims({'lon':'xlon'}) ds = ds.rename_dims({'lat':'ylat'}) ds = ds.rename_vars({'lat':'ylat'}) ds = ds.rename_vars({'lon':'xlon'}) + + t0 = time() + logger.info(f'Start writing nc file!') ds.to_netcdf(foutname, 'w', unlimited_dims='time', encoding={'temperature':{'dtype': 'h', '_FillValue': -30000.,'scale_factor': 0.001, 'add_offset': 20., 'missing_value': -30000.}}) ds.close() + logger.info(f'It took {time()-t0} seconds to write nc file!') #link output file to Fortran required files dir_path = os.path.abspath(outdir) diff --git a/pyschism/forcing/nws/nws2/era5.py b/pyschism/forcing/nws/nws2/era5.py index f313d79d..930d085c 100644 --- a/pyschism/forcing/nws/nws2/era5.py +++ b/pyschism/forcing/nws/nws2/era5.py @@ -33,28 +33,31 @@ def __init__(self, start_date=None, rnday: Union[float, timedelta] = 4, bbox=Non if tmpdir is not None: self.tmpdir = tmpdir - r = self.client.retrieve( - 'reanalysis-era5-single-levels', - { - 'variable':[ - '10m_u_component_of_wind','10m_v_component_of_wind','mean_sea_level_pressure', - '2m_dewpoint_temperature','2m_temperature','mean_total_precipitation_rate', - 'mean_surface_downward_long_wave_radiation_flux','mean_surface_downward_short_wave_radiation_flux' - ], - 'product_type':'reanalysis', - 'date':f"{self.start_date.strftime('%Y-%m-%d')}/{self.end_date.strftime('%Y-%m-%d')}", - 'time':[ - '00:00','01:00','02:00','03:00','04:00','05:00', - '06:00','07:00','08:00','09:00','10:00','11:00', - '12:00','13:00','14:00','15:00','16:00','17:00', - '18:00','19:00','20:00','21:00','22:00','23:00' - ], - 'area': [self._bbox.ymax+0.5, self._bbox.xmin-0.5, self._bbox.ymin-0.5, self._bbox.xmax+0.5], # North, West, South, East. Default: global - 'format':'netcdf' - }) - filename = self.tmpdir / f"era5_{self.start_date.strftime('%Y%m%d')}.nc" + if filename.is_file() == False: + + r = self.client.retrieve( + 'reanalysis-era5-single-levels', + { + 'variable':[ + '10m_u_component_of_wind','10m_v_component_of_wind','mean_sea_level_pressure', + '2m_dewpoint_temperature','2m_temperature','mean_total_precipitation_rate', + 'mean_surface_downward_long_wave_radiation_flux','mean_surface_downward_short_wave_radiation_flux' + ], + 'product_type':'reanalysis', + 'date':f"{self.start_date.strftime('%Y-%m-%d')}/{self.end_date.strftime('%Y-%m-%d')}", + 'time':[ + '00:00','01:00','02:00','03:00','04:00','05:00', + '06:00','07:00','08:00','09:00','10:00','11:00', + '12:00','13:00','14:00','15:00','16:00','17:00', + '18:00','19:00','20:00','21:00','22:00','23:00' + ], + 'area': [self._bbox.ymax+0.5, self._bbox.xmin-0.5, self._bbox.ymin-0.5, self._bbox.xmax+0.5], # North, West, South, East. Default: global + 'format':'netcdf' + } + ) + r.download(filename) @property From 5b5ded0d98ffb0776a1a99538cb33f99d1ffa1a8 Mon Sep 17 00:00:00 2001 From: cuill Date: Tue, 11 Apr 2023 17:35:21 -0400 Subject: [PATCH 30/31] update elev_ic --- examples/Gr3/Elev_IC/gen_elev.py | 9 ++++++ pyschism/mesh/gridgr3.py | 53 +++----------------------------- 2 files changed, 14 insertions(+), 48 deletions(-) create mode 100644 examples/Gr3/Elev_IC/gen_elev.py diff --git a/examples/Gr3/Elev_IC/gen_elev.py b/examples/Gr3/Elev_IC/gen_elev.py new file mode 100644 index 00000000..40f96915 --- /dev/null +++ b/examples/Gr3/Elev_IC/gen_elev.py @@ -0,0 +1,9 @@ +from pyschism.mesh.base import Gr3 +from pyschism.mesh.gridgr3 import ElevIc + +if __name__ == '__main__': + + hgrid = Gr3.open('./hgrid.gr3', crs='epsg:4326') + + elevic = ElevIc.default(hgrid=hgrid, offset=0.1) + elevic.write('elev.ic', overwrite=True) diff --git a/pyschism/mesh/gridgr3.py b/pyschism/mesh/gridgr3.py index a6fe6ae0..a3606b21 100644 --- a/pyschism/mesh/gridgr3.py +++ b/pyschism/mesh/gridgr3.py @@ -209,57 +209,14 @@ class IcField(Gr3Field): pass -class ElevIc(IcField): +class ElevIc(Gr3): @classmethod - def default(cls, hgrid, offset=-0.1): - obj = cls.constant(hgrid, 0.0) - mask = np.logical_and(hgrid.values > 0.0, obj.values < hgrid.values) - idxs = np.where(mask) - obj.values[idxs] = hgrid.values[idxs] + offset + def default(cls, hgrid, offset=0.1): + #Use Gr3.open to read hgrid.gr3 thus keep original bathymetry + obj = hgrid.copy() + obj.values[:] = np.maximum(0, -hgrid.values - offset) return obj - @classmethod - def modify_by_region(cls, hgrid, region=None, lonc = -77.07, latc = 24.0, offset=0.1): - #convert hgrid from lon/lat to cpp - hgrid = hgrid.copy() - hgrid.nodes.transform_to_cpp(lonc, latc) - #xy = hgrid.nodes.coords - #x = xy[:, 0] - #y = xy[:, 1] - - obj = cls.constant(hgrid, 0) - - if region is not None: - lines=[line.strip().split() for line in open(region, 'r').readlines()] - data=np.squeeze(np.array([lines[3:]])).astype('float') - x=data[:,0] - y=data[:,1] - coords = list( zip(x, y)) - poly = Polygon(coords) - - #region is in cpp projection - gdf1 = gpd.GeoDataFrame( - {'geometry': [poly]}) - - points = [Point(*coord) for coord in hgrid.nodes.coords] - gdf2 = gpd.GeoDataFrame( - {'geometry': points, 'index': list(range(len(points)))}) - gdf_in = gpd.sjoin(gdf2, gdf1, op="within") - picks = [i.index for i in gdf_in.itertuples()] - - #generate include.gr3 - obj = cls.constant(hgrid, 0) - obj.values[picks] = 1 - obj.description = 'include.reg, in: 1, out: 0' - obj.write('include.gr3', overwrite=True) - - elevic = hgrid.copy() - idxs = np.where(obj.values == 0) - elevic.nodes.values[idxs] = np.maximum(0, -elevic.nodes.values[idxs]) - idxs = np.where(obj.values == 1) - elevic.nodes.values[idxs] = np.maximum(0, -elevic.nodes.values[idxs]-offset) - return elevic - class TempIc(IcField): @classmethod def from_hycom(cls, gr3: Gr3, hycom: Hycom, date): From 744502d807a366976505b6ae0ab4a6d0cc211740 Mon Sep 17 00:00:00 2001 From: cuill Date: Thu, 13 Apr 2023 10:11:10 -0400 Subject: [PATCH 31/31] only select conus layers --- pyschism/forcing/source_sink/nwm.py | 22 +++++++++++++--------- pyschism/utils/jsonencoder.py | 13 +++++++++++++ 2 files changed, 26 insertions(+), 9 deletions(-) create mode 100644 pyschism/utils/jsonencoder.py diff --git a/pyschism/forcing/source_sink/nwm.py b/pyschism/forcing/source_sink/nwm.py index 05ca8a68..c462854f 100644 --- a/pyschism/forcing/source_sink/nwm.py +++ b/pyschism/forcing/source_sink/nwm.py @@ -31,6 +31,7 @@ from pyschism import dates from pyschism.mesh.base import Gr3 +from pyschism.utils.jsonencoder import NpEncoder from pyschism.forcing.source_sink.base import SourceSink, Sources, Sinks @@ -39,7 +40,6 @@ logger = logging.getLogger(__name__) - class NWMElementPairings: def __init__(self, hgrid: Gr3, nwm_file=None, workers=-1): @@ -198,13 +198,13 @@ def save_json(self, sources=True, sinks=True): if sources: logger.info(f"Saving {sources}") with open(sources, "w") as fh: - json.dump(self.sources, fh) + json.dump(self.sources, fh, cls=NpEncoder) sinks = "sinks.json" if sinks is True else sinks if sinks: with open(sinks, "w") as fh: logger.info(f"Saving {sinks}") - json.dump(self.sinks, fh) + json.dump(self.sinks, fh, cls=NpEncoder) @staticmethod def load_json(hgrid, sources=None, sinks=None): @@ -261,11 +261,14 @@ def hgrid(self): def gdf(self): if not hasattr(self, "_gdf"): gdf_coll = [] - for reach_layer in [ - reach_layer - for reach_layer in fiona.listlayers(self.nwm_file) - if "reaches" in reach_layer - ]: + #for reach_layer in [ + # reach_layer + # for reach_layer in fiona.listlayers(self.nwm_file) + # if "reaches" in reach_layer + #]: + reach_layers = ['nwm_reaches_conus'] + for reach_layer in reach_layers: + logger.info(f'layer is {reach_layer}') layer_crs = gpd.read_file( self.nwm_file, rows=1, layer=reach_layer).crs bbox = self.hgrid.get_bbox(crs=layer_crs) @@ -276,7 +279,8 @@ def gdf(self): layer=reach_layer, ) ) - self._gdf = pd.concat(gdf_coll) + #self._gdf = pd.concat(gdf_coll) + self._gdf = gpd.GeoDataFrame(gdf_coll[0]) return self._gdf @property diff --git a/pyschism/utils/jsonencoder.py b/pyschism/utils/jsonencoder.py new file mode 100644 index 00000000..b2822b49 --- /dev/null +++ b/pyschism/utils/jsonencoder.py @@ -0,0 +1,13 @@ +import json +import numpy as np + +#https://java2blog.com/object-of-type-int64-is-not-json-serializable/ +class NpEncoder(json.JSONEncoder): + def default(self, obj): + if isinstance(obj, np.integer): + return int(obj) + if isinstance(obj, np.floating): + return float(obj) + if isinstance(obj, np.ndarray): + return obj.tolist() + return super(NpEncoder, self).default(obj)