From 9bb4558e4352a55eb7edd95a5a3d9abb4e58aee9 Mon Sep 17 00:00:00 2001
From: thopri <thopri@noc.ac.uk>
Date: Mon, 6 Apr 2020 11:56:36 +0100
Subject: [PATCH] Implemented FES hc extraction
---
inputs/namelist_remote.bdy | 3 +-
pynemo/tide/fes_extract_HC.py | 71 ++++++++-------
pynemo/tide/nemo_bdy_tide.py | 161 ---------------------------------
pynemo/tide/nemo_bdy_tide2.py | 1 -
pynemo/tide/nemo_bdy_tide3.py | 39 +++++---
pynemo/tide/tpxo_extract_HC.py | 4 +-
6 files changed, 65 insertions(+), 214 deletions(-)
delete mode 100644 pynemo/tide/nemo_bdy_tide.py
delete mode 100644 pynemo/tide/nemo_bdy_tide2.py
diff --git a/inputs/namelist_remote.bdy b/inputs/namelist_remote.bdy
index eef5337..7768005 100644
--- a/inputs/namelist_remote.bdy
+++ b/inputs/namelist_remote.bdy
@@ -70,8 +70,7 @@
ln_tide = .true. ! =T : produce bdy tidal conditions
sn_tide_model = 'fes' ! Name of tidal model (fes|tpxo)
clname(1) = 'M2' ! constituent name
- clname(2) = 'S2'
- clname(3) = 'K2'
+ clname(2) = 'S2'
ln_trans = .true. ! interpolate transport rather than
! velocities
!------------------------------------------------------------------------------
diff --git a/pynemo/tide/fes_extract_HC.py b/pynemo/tide/fes_extract_HC.py
index 4372cd5..d8cc3fb 100644
--- a/pynemo/tide/fes_extract_HC.py
+++ b/pynemo/tide/fes_extract_HC.py
@@ -12,7 +12,7 @@ from netCDF4 import Dataset
from scipy import interpolate
import numpy as np
-class FesExtract(object):
+class HcExtract(object):
""" This is FES model extract_hc.c implementation in python adapted from tpxo_extract_HC.py"""
def __init__(self, settings, lat, lon, grid_type):
"""initialises the Extract of tide information from the netcdf
@@ -40,11 +40,13 @@ class FesExtract(object):
#constituents = ['2N2', 'EPS2', 'J1', 'K1', 'K2', 'L2', 'LA2', 'M2', 'M3', 'M4', 'M6', 'M8', 'MF', 'MKS2',
#'MM', 'MN4', 'MS4', 'MSF', 'MSQM', 'MTM', 'MU2', 'N2', 'N4', 'NU2', 'O1', 'P1', 'Q1', 'R2',
#'S1', 'S2', 'S4', 'SA', 'SSA', 'T2']
- self.constituents = ['M2','S2']
+ constituents = ['M2','S2']
+
+ self.cons = constituents
# extract lon and lat z data
- lon_z = np.array(Dataset(settings['tide_fes']+self.constituents[0]+'_Z.nc').variables['lon'])
- lat_z = np.array(Dataset(settings['tide_fes']+self.constituents[0]+'_Z.nc').variables['lat'])
+ lon_z = np.array(Dataset(settings['tide_fes']+constituents[0]+'_Z.nc').variables['lon'])
+ lat_z = np.array(Dataset(settings['tide_fes']+constituents[0]+'_Z.nc').variables['lat'])
lon_resolution = lon_z[1] - lon_z[0]
data_in_km = 0 # added to maintain the reference to matlab tmd code
@@ -55,21 +57,21 @@ class FesExtract(object):
self.mask_dataset = {}
# extract example amplitude grid for Z, U and V and change NaNs to 0 (for land) and other values to 1 (for water)
- mask_z = np.array(np.rot90(Dataset(settings['tide_fes']+self.constituents[0]+'_Z.nc').variables['amplitude'][:]))
+ mask_z = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_Z.nc').variables['amplitude'][:]))
where_are_NaNs = np.isnan(mask_z)
where_no_NaNs = np.invert(np.isnan(mask_z))
mask_z[where_no_NaNs] = 1
mask_z[where_are_NaNs] = 0
self.mask_dataset[mz_name] = mask_z
- mask_u = np.array(np.rot90(Dataset(settings['tide_fes']+self.constituents[0]+'_U.nc').variables['Ua'][:]))
+ mask_u = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_U.nc').variables['Ua'][:]))
where_are_NaNs = np.isnan(mask_u)
where_no_NaNs = np.invert(np.isnan(mask_u))
mask_u[where_no_NaNs] = 1
mask_u[where_are_NaNs] = 0
self.mask_dataset[mu_name] = mask_u
- mask_v = np.array(np.rot90(Dataset(settings['tide_fes']+self.constituents[0]+'_V.nc').variables['Va'][:]))
+ mask_v = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_V.nc').variables['Va'][:]))
where_are_NaNs = np.isnan(mask_v)
where_no_NaNs = np.invert(np.isnan(mask_v))
mask_v[where_no_NaNs] = 1
@@ -77,33 +79,33 @@ class FesExtract(object):
self.mask_dataset[mv_name] = mask_v
#read and convert the height_dataset file to complex and store in dicts
- for ncon in range(len(self.constituents)):
- amp = np.array(np.rot90(Dataset(settings['tide_fes']+str(self.constituents[ncon])+'_Z.nc').variables['amplitude'][:]))
- phase = np.array(np.rot90(Dataset(settings['tide_fes']+self.constituents[ncon]+'_Z.nc').variables['phase'][:]))
- lat_z = np.array(Dataset(settings['tide_fes']+self.constituents[ncon]+'_Z.nc').variables['lat'][:])
- lon_z = np.array(Dataset(settings['tide_fes']+self.constituents[ncon]+'_Z.nc').variables['lon'][:])
+ for ncon in range(len(constituents)):
+ amp = np.array(np.rot90(Dataset(settings['tide_fes']+str(constituents[ncon])+'_Z.nc').variables['amplitude'][:]))
+ phase = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_Z.nc').variables['phase'][:]))
+ lat_z = np.array(Dataset(settings['tide_fes']+constituents[ncon]+'_Z.nc').variables['lat'][:])
+ lon_z = np.array(Dataset(settings['tide_fes']+constituents[ncon]+'_Z.nc').variables['lon'][:])
hRe = amp*np.sin(phase)
hIm = amp*np.cos(phase)
- self.height_dataset[self.constituents[ncon]] = {'lat_z':lat_z,'lon_z':lon_z,'hRe':hRe,'hIm':hIm}
+ self.height_dataset[constituents[ncon]] = {'lat_z':lat_z,'lon_z':lon_z,'hRe':hRe,'hIm':hIm}
#read and convert the velocity_dataset files to complex
- for ncon in range(len(self.constituents)):
- amp = np.array(np.rot90(Dataset(settings['tide_fes']+self.constituents[ncon]+'_U.nc').variables['Ua'][:]))
- phase = np.array(np.rot90(Dataset(settings['tide_fes']+self.constituents[ncon]+'_U.nc').variables['Ug'][:]))
- lat_u = np.array(Dataset(settings['tide_fes']+self.constituents[ncon]+'_U.nc').variables['lat'][:])
- lon_u = np.array(Dataset(settings['tide_fes']+self.constituents[ncon]+'_U.nc').variables['lon'][:])
+ for ncon in range(len(constituents)):
+ amp = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ua'][:]))
+ phase = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ug'][:]))
+ lat_u = np.array(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['lat'][:])
+ lon_u = np.array(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['lon'][:])
URe = amp*np.sin(phase)
UIm = amp*np.cos(phase)
- self.Uvelocity_dataset[self.constituents[ncon]] = {'lat_u':lat_u,'lon_u':lon_u,'URe':URe,'UIm':UIm}
+ self.Uvelocity_dataset[constituents[ncon]] = {'lat_u':lat_u,'lon_u':lon_u,'URe':URe,'UIm':UIm}
- for ncon in range(len(self.constituents)):
- amp = np.array(np.rot90(Dataset(settings['tide_fes']+self.constituents[ncon]+'_V.nc').variables['Va'][:]))
- phase = np.array(np.rot90(Dataset(settings['tide_fes']+self.constituents[ncon]+'_V.nc').variables['Vg'][:]))
- lat_v = np.array(Dataset(settings['tide_fes']+self.constituents[ncon]+'_V.nc').variables['lat'][:])
- lon_v = np.array(Dataset(settings['tide_fes']+self.constituents[ncon]+'_V.nc').variables['lon'][:])
+ for ncon in range(len(constituents)):
+ amp = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Va'][:]))
+ phase = np.array(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Vg'][:]))
+ lat_v = np.array(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['lat'][:])
+ lon_v = np.array(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['lon'][:])
VRe = amp*np.sin(phase)
VIm = amp*np.cos(phase)
- self.Vvelocity_dataset[self.constituents[ncon]] = {'lat_v':lat_v,'lon_v':lon_v,'VRe':VRe,'VIm':VIm}
+ self.Vvelocity_dataset[constituents[ncon]] = {'lat_v':lat_v,'lon_v':lon_v,'VRe':VRe,'VIm':VIm}
# open grid variables these are resampled TPXO parameters so may not work correctly.
@@ -169,11 +171,11 @@ class FesExtract(object):
hRe_name, hIm_name, lon_z_name, lat_z_name,
lon, lat, maskname=mz_name)
elif grid_type == 'u':
- self.amp, self.gph = self.interpolate_constituents(self.velocity_dataset,
+ self.amp, self.gph = self.interpolate_constituents(self.Uvelocity_dataset,
URe_name, UIm_name, lon_u_name, lat_u_name,
lon, lat, depth, maskname=mu_name)
elif grid_type == 'v':
- self.amp, self.gph = self.interpolate_constituents(self.velocity_dataset,
+ self.amp, self.gph = self.interpolate_constituents(self.Vvelocity_dataset,
VRe_name, VIm_name, lon_v_name, lat_v_name,
lon, lat, depth, maskname=mv_name)
else:
@@ -183,13 +185,16 @@ class FesExtract(object):
def interpolate_constituents(self, nc_dataset, real_var_name, img_var_name, lon_var_name,
lat_var_name, lon, lat, height_data=None, maskname=None):
""" Interpolates the tidal constituents along the given lat lon coordinates """
- amp = np.zeros((len(nc_dataset), nc_dataset[list(nc_dataset)[0]]['lon_z'].shape[0]))
- gph = np.zeros((len(nc_dataset), nc_dataset[list(nc_dataset)[0]]['lon_z'].shape[0]))
+ amp = np.zeros((len(nc_dataset), lon.shape[0]))
+ gph = np.zeros((len(nc_dataset), lon.shape[0]))
- data = np.array(np.ravel(nc_dataset['M2'][real_var_name]), dtype=complex)
- data.imag = np.array(np.ravel(nc_dataset['M2'][img_var_name]))
+ # TODO: need to sort multiple HC, at the momemnt it uses M2 for every harmonic
- data = data.reshape(nc_dataset['M2'][real_var_name].shape)
+ data = np.array((nc_dataset['M2'][real_var_name]), dtype=complex)
+ data.imag = np.array((nc_dataset['M2'][img_var_name]))
+ # add extra dim to be compatable with adapted code that expects a list of HC
+ data = np.expand_dims(data,axis=0)
+ #data = data.reshape(1,nc_dataset['M2'][real_var_name].shape)
# data[data==0] = np.NaN
# Lat Lon values
@@ -220,7 +225,7 @@ class FesExtract(object):
maskedpoints = interpolate.interpn((x_values, y_values), mask, lonlat)
data_temp = np.zeros((data.shape[0], lon.shape[0], 2, ))
- #check at same point in TPXO extract HC script.
+ #check at same point in TPXO extract HC sc
for cons_index in range(data.shape[0]):
#interpolate real values
data_temp[cons_index, :, 0] = interpolate_data(x_values, y_values,
diff --git a/pynemo/tide/nemo_bdy_tide.py b/pynemo/tide/nemo_bdy_tide.py
deleted file mode 100644
index 0fdf5ab..0000000
--- a/pynemo/tide/nemo_bdy_tide.py
+++ /dev/null
@@ -1,161 +0,0 @@
-'''
-
-
-
-'''
-import numpy as np
-import scipy.spatial as sp
-from netCDF4 import Dataset
-import copy # DEBUG ONLY- allows multiple runs without corruption
-import nemo_bdy_grid_angle
-#from nemo_bdy_extr_tm3 import rot_rep
-
-class Extract:
-
- def __init__(self, setup, DstCoord, Grid):
-
- self.g_type = Grid.grid_type
- DC = copy.deepcopy(DstCoord)
- dst_lon = DC.bdy_lonlat[self.g_type]['lon'][Grid.bdy_r == 0]
- dst_lat = DC.bdy_lonlat[self.g_type]['lat'][Grid.bdy_r == 0]
- self.dst_dep = DC.depths[self.g_type]['bdy_hbat'][Grid.bdy_r == 0]
- self.harm_Im = {} # tidal boundary data: Imaginary
- self.harm_Re = {} # tidal boundary data: Real
-
- # Modify lon for 0-360 TODO this needs to be auto-dectected
-
- dst_lon = np.array([x if x > 0 else x+360 for x in dst_lon])
-
- fileIDb = '/Users/jdha/Projects/pynemo_data/DATA/grid_tpxo7.2.nc' # TPX bathymetry file
- nb = Dataset(fileIDb) # Open the TPX bathybetry file using the NetCDF4-Python library
-
- # Open the TPX Datafiles using the NetCDF4-Python library
-# T_GridAngles = nemo_bdy_grid_angle.GridAngle(
-# self.settings['src_hgr'], imin, imax, jmin, jmax, 't')
-# RotStr_GridAngles = nemo_bdy_grid_angle.GridAngle(
-# self.settings['dst_hgr'], 1, maxI, 1, maxJ, self.rot_str)
-
-# self.gcos = T_GridAngles.cosval
-# self.gsin = T_GridAngles.sinval
-
- if self.g_type == 't':
- self.fileID = '/Users/jdha/Projects/pynemo_data/DATA/h_tpxo7.2.nc' # TPX sea surface height file
- self.var_Im = 'hIm'
- self.var_Re = 'hRe'
- nc = Dataset(self.fileID) # pass variable ids to nc
- lon = np.ravel(nc.variables['lon_z'][:,:]) # need to add in a east-west wrap-around
- lat = np.ravel(nc.variables['lat_z'][:,:])
- bat = np.ravel(nb.variables['hz'][:,:])
- msk = np.ravel(nb.variables['mz'][:,:])
- elif self.g_type == 'u':
- self.fileID = '/Users/jdha/Projects/pynemo_data/DATA/u_tpxo7.2.nc' # TPX velocity file
- self.var_Im = 'UIm'
- self.var_Re = 'URe'
- self.key_tr = setup['tide_trans']
- nc = Dataset(self.fileID) # pass variable ids to nc
- lon = np.ravel(nc.variables['lon_u'][:,:])
- lat = np.ravel(nc.variables['lat_u'][:,:])
- bat = np.ravel(nb.variables['hu'][:,:])
- msk = np.ravel(nb.variables['mu'][:,:])
- else:
- self.fileID = '/Users/jdha/Projects/pynemo_data/DATA/u_tpxo7.2.nc' # TPX velocity file
- self.var_Im = 'VIm'
- self.var_Re = 'VRe'
- self.key_tr = setup['tide_trans']
- nc = Dataset(self.fileID) # pass variable ids to nc
- lon = np.ravel(nc.variables['lon_v'][:,:])
- lat = np.ravel(nc.variables['lat_v'][:,:])
- bat = np.ravel(nb.variables['hv'][:,:])
- msk = np.ravel(nb.variables['mv'][:,:])
-
- # Pull out the constituents that are avaibable
- self.cons = []
- for ncon in range(nc.variables['con'].shape[0]):
- self.cons.append(nc.variables['con'][ncon,:].tostring().strip())
-
- nc.close() # Close Datafile
- nb.close() # Close Bathymetry file
-
- # Find nearest neighbours on the source grid to each dst bdy point
- source_tree = sp.cKDTree(list(zip(lon, lat)))
- dst_pts = list(zip(dst_lon, dst_lat))
- nn_dist, self.nn_id = source_tree.query(dst_pts, k=4, eps=0, p=2,
- distance_upper_bound=0.5)
-
- # Create a weighting index for interpolation onto dst bdy point
- # need to check for missing values
-
- ind = nn_dist == np.inf
-
- self.nn_id[ind] = 0 # better way of carrying None in the indices?
- dx = (lon[self.nn_id] - np.repeat(np.reshape(dst_lon,[dst_lon.size, 1]),4,axis=1) ) * np.cos(np.repeat(np.reshape(dst_lat,[dst_lat.size, 1]),4,axis=1) * np.pi / 180.)
- dy = lat[self.nn_id] - np.repeat(np.reshape(dst_lat,[dst_lat.size, 1]),4,axis=1)
-
- dist_tot = np.power((np.power(dx, 2) + np.power(dy, 2)), 0.5)
-
- self.msk = msk[self.nn_id]
- self.bat = bat[self.nn_id]
-
- dist_tot[ind | self.msk] = np.nan
-
- dist_wei = 1/( np.divide(dist_tot,(np.repeat(np.reshape(np.nansum(dist_tot,axis=1),[dst_lat.size, 1]),4,axis=1)) ) )
-
- self.nn_wei = dist_wei/np.repeat(np.reshape(np.nansum(dist_wei, axis=1),[dst_lat.size, 1]),4,axis=1)
- self.nn_wei[ind | self.msk] = 0.
-
- # Need to identify missing points and throw a warning and set values to zero
-
- mv = np.sum(self.wei,axis=1) == 0
- print('##WARNING## There are', np.sum(mv), 'missing values, these will be set to ZERO')
-
- def extract_con(self, con):
-
- if con in self.cons:
- con_ind = self.cons.index(con)
-
- # Extract the complex amplitude components
-
- nc = Dataset(self.fileID) # pass variable ids to nc
-
- vIm = np.ravel(nc.variables[self.var_Im][con_ind,:,:])
- vRe = np.ravel(nc.variables[self.var_Re][con_ind,:,:])
-
- nc.close()
-
- if self.g_type != 't':
-
- self.harm_Im[con] = np.sum(vIm[self.nn_id]*self.nn_wei,axis=1)
- self.harm_Re[con] = np.sum(vRe[self.nn_id]*self.nn_wei,axis=1)
-
- else: # Convert transports to velocities
-
- if self.key_tr == True: # We convert to velocity using tidal model bathymetry
-
- self.harm_Im[con] = np.sum(vIm[self.nn_id]*self.nn_wei,axis=1)/np.sum(self.bat[self.nn_id]*self.nn_wei,axis=1)
- self.harm_Re[con] = np.sum(vRe[self.nn_id]*self.nn_wei,axis=1)/np.sum(self.bat[self.nn_id]*self.nn_wei,axis=1)
-
- else: # We convert to velocity using the regional model bathymetry
-
- self.harm_Im[con] = np.sum(vIm[self.nn_id]*self.nn_wei,axis=1)/self.dst_dep
- self.harm_Re[con] = np.sum(vRe[self.nn_id]*self.nn_wei,axis=1)/self.dst_dep
-
-
- # Rotate vectors
-
- self.harm_Im_rot[con] = self.rot_rep(self.harm_Im[con], self.harm_Im[con], self.rot_str,
- 'en to %s' %self.rot_dir, self.dst_gcos, self.dst_gsin)
- self.harm_Re_rot[con] = self.rot_rep(self.harm_Re[con], self.harm_Re[con], self.rot_str,
- 'en to %s' %self.rot_dir, self.dst_gcos, self.dst_gsin)
-
- else:
-
- # throw some warning
- print('##WARNING## Missing constituent values will be set to ZERO')
-
- self.harm_Im[con] = np.zeros(self.nn_id[:,0].size)
- self.harm_Re[con] = np.zeros(self.nn_id[:,0].size)
-
-
-
-
-
diff --git a/pynemo/tide/nemo_bdy_tide2.py b/pynemo/tide/nemo_bdy_tide2.py
deleted file mode 100644
index 581a012..0000000
--- a/pynemo/tide/nemo_bdy_tide2.py
+++ /dev/null
@@ -1 +0,0 @@
-Non
\ No newline at end of file
diff --git a/pynemo/tide/nemo_bdy_tide3.py b/pynemo/tide/nemo_bdy_tide3.py
index 61adddf..ac89e54 100644
--- a/pynemo/tide/nemo_bdy_tide3.py
+++ b/pynemo/tide/nemo_bdy_tide3.py
@@ -31,13 +31,17 @@ def nemo_bdy_tide_rot(setup, DstCoord, Grid_T, Grid_U, Grid_V, comp,tide_model):
nbdyu = len(Grid_U.bdy_i)
nbdyv = len(Grid_V.bdy_i)
+ # TODO: change from if statement defining HC extract to string passed that defines HC extract script
+ # e.g. pass 'tpxo' for tpxo_extract_HC.py or 'fes' fro fes_extract_HC.py. This will make easier to add new
+ # databases of HC
+
#convert the dst_lon into TMD Conventions (0E/360E)
dst_lon[dst_lon < 0.0] = dst_lon[dst_lon < 0.0]+360.0
#extract the surface elevation at each z-point
if tide_model == 'tpxo':
- tpxo_z = tpxo_extract_HC.TpxoExtract(setup.settings, dst_lat, dst_lon, g_type)
+ tpxo_z = tpxo_extract_HC.HcExtract(setup.settings, dst_lat, dst_lon, g_type)
if tide_model == 'fes':
- fes_z = fes_extract_HC.FesExtract(setup.settings,dst_lat,dst_lon,g_type)
+ fes_z = fes_extract_HC.HcExtract(setup.settings,dst_lat,dst_lon,g_type)
#convert back the z-longitudes into the usual conventions (-180E/+180E)
dst_lon[dst_lon > 180.0] = dst_lon[dst_lon > 180.0]-360.0
@@ -65,8 +69,8 @@ def nemo_bdy_tide_rot(setup, DstCoord, Grid_T, Grid_U, Grid_V, comp,tide_model):
#convert the U-longitudes into the TMD conventions (0/360E)
dst_lon[dst_lon < 0.0] = dst_lon[dst_lon < 0.0]+360.0
if tide_model == 'tpxo':
- tpxo_ux = tpxo_extract_HC.TpxoExtract(setup.settings, dst_lat, dst_lon, Grid_U.grid_type)
- tpxo_vx = tpxo_extract_HC.TpxoExtract(setup.settings, dst_lat, dst_lon, Grid_V.grid_type)
+ tpxo_ux = tpxo_extract_HC.HcExtract(setup.settings, dst_lat, dst_lon, Grid_U.grid_type)
+ tpxo_vx = tpxo_extract_HC.HcExtract(setup.settings, dst_lat, dst_lon, Grid_V.grid_type)
ampuX = tpxo_ux.amp
phauX = tpxo_ux.gph
@@ -74,13 +78,13 @@ def nemo_bdy_tide_rot(setup, DstCoord, Grid_T, Grid_U, Grid_V, comp,tide_model):
phavX = tpxo_vx.gph
if tide_model == 'fes':
- fes_uy = fes_extract_HC.FesExtract(setup.settings, dst_lat, dst_lon, Grid_U.grid_type)
- fes_vy = fes_extract_HC.FesExtract(setup.settings, dst_lat, dst_lon, Grid_V.grid_type)
+ fes_ux = fes_extract_HC.HcExtract(setup.settings, dst_lat, dst_lon, Grid_U.grid_type)
+ fes_vx = fes_extract_HC.HcExtract(setup.settings, dst_lat, dst_lon, Grid_V.grid_type)
- ampuY = fes_uy.amp
- phauY = fes_uy.gph
- ampvY = fes_vy.amp
- phavY = fes_vy.gph
+ ampuX = fes_ux.amp
+ phauX = fes_ux.gph
+ ampvX = fes_vx.amp
+ phavX = fes_vx.gph
#check if ux data are missing
ind = np.where((np.isnan(ampuX)) | (np.isnan(phauX)))
@@ -105,8 +109,8 @@ def nemo_bdy_tide_rot(setup, DstCoord, Grid_T, Grid_U, Grid_V, comp,tide_model):
#convert the U-longitudes into the TMD conventions (0/360E)
dst_lon[dst_lon < 0.0] = dst_lon[dst_lon < 0.0]+360.0
if tide_model == 'tpxo':
- tpxo_uy = tpxo_extract_HC.TpxoExtract(setup.settings, dst_lat, dst_lon, Grid_U.grid_type)
- tpxo_vy = tpxo_extract_HC.TpxoExtract(setup.settings, dst_lat, dst_lon, Grid_V.grid_type)
+ tpxo_uy = tpxo_extract_HC.HcExtract(setup.settings, dst_lat, dst_lon, Grid_U.grid_type)
+ tpxo_vy = tpxo_extract_HC.HcExtract(setup.settings, dst_lat, dst_lon, Grid_V.grid_type)
ampuY = tpxo_uy.amp
phauY = tpxo_uy.gph
@@ -114,8 +118,8 @@ def nemo_bdy_tide_rot(setup, DstCoord, Grid_T, Grid_U, Grid_V, comp,tide_model):
phavY = tpxo_vy.gph
if tide_model == 'fes':
- fes_uy = fes_extract_HC.FesExtract(setup.settings, dst_lat, dst_lon, Grid_U.grid_type)
- fes_vy = fes_extract_HC.FesExtract(setup.settings, dst_lat, dst_lon, Grid_V.grid_type)
+ fes_uy = fes_extract_HC.HcExtract(setup.settings, dst_lat, dst_lon, Grid_U.grid_type)
+ fes_vy = fes_extract_HC.HcExtract(setup.settings, dst_lat, dst_lon, Grid_V.grid_type)
ampuY = fes_uy.amp
phauY = fes_uy.gph
@@ -208,7 +212,12 @@ def nemo_bdy_tide_rot(setup, DstCoord, Grid_T, Grid_U, Grid_V, comp,tide_model):
cosvY = np.zeros((numharm, nbdyv))
sinvY = np.zeros((numharm, nbdyv))
- compindx = constituents_index(tpxo_z.cons, comp)
+ if tide_model == 'tpxo':
+ compindx = constituents_index(tpxo_z.cons, comp)
+
+ if tide_model == 'fes':
+ compindx = constituents_index(fes_z.cons, comp)
+
for h in range(0, numharm):
c = int(compindx[h])
if c != -1:
diff --git a/pynemo/tide/tpxo_extract_HC.py b/pynemo/tide/tpxo_extract_HC.py
index 10cf627..1943c53 100644
--- a/pynemo/tide/tpxo_extract_HC.py
+++ b/pynemo/tide/tpxo_extract_HC.py
@@ -12,7 +12,7 @@ from netCDF4 import Dataset
from scipy import interpolate
import numpy as np
-class TpxoExtract(object):
+class HcExtract(object):
""" This is TPXO model extract_hc.c implementation in python"""
def __init__(self, settings, lat, lon, grid_type):
"""initialises the Extract of tide information from the netcdf
@@ -51,7 +51,7 @@ class TpxoExtract(object):
# Pull out the constituents that are avaibable
self.cons = []
for ncon in range(self.height_dataset.variables['con'].shape[0]):
- self.cons.append(self.height_dataset.variables['con'][ncon, :].tostring().strip())
+ self.cons.append(self.height_dataset.variables['con'][ncon, :].tostring().strip().decode())
elif tide_model == 'FES':
constituents = ['2N2','EPS2','J1','K1','K2','L2','LA2','M2','M3','M4','M6','M8','MF','MKS2','MM','MN4','MS4','MSF','MSQM','MTM','MU2','N2','N4','NU2','O1','P1','Q1','R2','S1','S2','S4','SA','SSA','T2']
print('did not actually code stuff for FES in this routine. Though that would be ideal. Instead put it in fes_extract_HC.py')
--
GitLab