Skip to content
GitLab
Projects
Groups
Snippets
Help
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Open sidebar
thopri
PyNEMO
Commits
9bb4558e
Commit
9bb4558e
authored
5 years ago
by
thopri
Browse files
Options
Download
Email Patches
Plain Diff
Implemented FES hc extraction
parent
6cf0fc2e
Changes
6
Hide whitespace changes
Inline
Side-by-side
Showing
6 changed files
with
65 additions
and
214 deletions
+65
-214
inputs/namelist_remote.bdy
inputs/namelist_remote.bdy
+1
-2
pynemo/tide/fes_extract_HC.py
pynemo/tide/fes_extract_HC.py
+38
-33
pynemo/tide/nemo_bdy_tide.py
pynemo/tide/nemo_bdy_tide.py
+0
-161
pynemo/tide/nemo_bdy_tide2.py
pynemo/tide/nemo_bdy_tide2.py
+0
-1
pynemo/tide/nemo_bdy_tide3.py
pynemo/tide/nemo_bdy_tide3.py
+24
-15
pynemo/tide/tpxo_extract_HC.py
pynemo/tide/tpxo_extract_HC.py
+2
-2
No files found.
inputs/namelist_remote.bdy
View file @
9bb4558e
...
...
@@ -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
!------------------------------------------------------------------------------
...
...
This diff is collapsed.
Click to expand it.
pynemo/tide/fes_extract_HC.py
View file @
9bb4558e
...
...
@@ -12,7 +12,7 @@ from netCDF4 import Dataset
from
scipy
import
interpolate
import
numpy
as
np
class
Fes
Extract
(
object
):
class
Hc
Extract
(
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
.
U
velocity_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
.
V
velocity_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 sc
ript.
#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
,
...
...
This diff is collapsed.
Click to expand it.
pynemo/tide/nemo_bdy_tide.py
deleted
100644 → 0
View file @
6cf0fc2e
'''
'''
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
)
This diff is collapsed.
Click to expand it.
pynemo/tide/nemo_bdy_tide2.py
deleted
100644 → 0
View file @
6cf0fc2e
Non
\ No newline at end of file
This diff is collapsed.
Click to expand it.
pynemo/tide/nemo_bdy_tide3.py
View file @
9bb4558e
...
...
@@ -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
.
Tpxo
Extract
(
setup
.
settings
,
dst_lat
,
dst_lon
,
g_type
)
tpxo_z
=
tpxo_extract_HC
.
Hc
Extract
(
setup
.
settings
,
dst_lat
,
dst_lon
,
g_type
)
if
tide_model
==
'fes'
:
fes_z
=
fes_extract_HC
.
Fes
Extract
(
setup
.
settings
,
dst_lat
,
dst_lon
,
g_type
)
fes_z
=
fes_extract_HC
.
Hc
Extract
(
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
.
Tpxo
Extract
(
setup
.
settings
,
dst_lat
,
dst_lon
,
Grid_U
.
grid_type
)
tpxo_vx
=
tpxo_extract_HC
.
Tpxo
Extract
(
setup
.
settings
,
dst_lat
,
dst_lon
,
Grid_V
.
grid_type
)
tpxo_ux
=
tpxo_extract_HC
.
Hc
Extract
(
setup
.
settings
,
dst_lat
,
dst_lon
,
Grid_U
.
grid_type
)
tpxo_vx
=
tpxo_extract_HC
.
Hc
Extract
(
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_u
y
=
fes_extract_HC
.
Fes
Extract
(
setup
.
settings
,
dst_lat
,
dst_lon
,
Grid_U
.
grid_type
)
fes_v
y
=
fes_extract_HC
.
Fes
Extract
(
setup
.
settings
,
dst_lat
,
dst_lon
,
Grid_V
.
grid_type
)
fes_u
x
=
fes_extract_HC
.
Hc
Extract
(
setup
.
settings
,
dst_lat
,
dst_lon
,
Grid_U
.
grid_type
)
fes_v
x
=
fes_extract_HC
.
Hc
Extract
(
setup
.
settings
,
dst_lat
,
dst_lon
,
Grid_V
.
grid_type
)
ampu
Y
=
fes_u
y
.
amp
phau
Y
=
fes_u
y
.
gph
ampv
Y
=
fes_v
y
.
amp
phav
Y
=
fes_v
y
.
gph
ampu
X
=
fes_u
x
.
amp
phau
X
=
fes_u
x
.
gph
ampv
X
=
fes_v
x
.
amp
phav
X
=
fes_v
x
.
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
.
Tpxo
Extract
(
setup
.
settings
,
dst_lat
,
dst_lon
,
Grid_U
.
grid_type
)
tpxo_vy
=
tpxo_extract_HC
.
Tpxo
Extract
(
setup
.
settings
,
dst_lat
,
dst_lon
,
Grid_V
.
grid_type
)
tpxo_uy
=
tpxo_extract_HC
.
Hc
Extract
(
setup
.
settings
,
dst_lat
,
dst_lon
,
Grid_U
.
grid_type
)
tpxo_vy
=
tpxo_extract_HC
.
Hc
Extract
(
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
.
Fes
Extract
(
setup
.
settings
,
dst_lat
,
dst_lon
,
Grid_U
.
grid_type
)
fes_vy
=
fes_extract_HC
.
Fes
Extract
(
setup
.
settings
,
dst_lat
,
dst_lon
,
Grid_V
.
grid_type
)
fes_uy
=
fes_extract_HC
.
Hc
Extract
(
setup
.
settings
,
dst_lat
,
dst_lon
,
Grid_U
.
grid_type
)
fes_vy
=
fes_extract_HC
.
Hc
Extract
(
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
:
...
...
This diff is collapsed.
Click to expand it.
pynemo/tide/tpxo_extract_HC.py
View file @
9bb4558e
...
...
@@ -12,7 +12,7 @@ from netCDF4 import Dataset
from
scipy
import
interpolate
import
numpy
as
np
class
Tpxo
Extract
(
object
):
class
Hc
Extract
(
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'
)
...
...
This diff is collapsed.
Click to expand it.
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment