/
check_points.py
223 lines (204 loc) · 8.65 KB
/
check_points.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
#!/usr/bin/env python
u"""
check_points.py
Written by Tyler Sutterley (04/2024)
Check if points are within a tide model domain
OTIS format tidal solutions provided by Ohio State University and ESR
http://volkov.oce.orst.edu/tides/region.html
https://www.esr.org/research/polar-tide-models/list-of-polar-tide-models/
ftp://ftp.esr.org/pub/datasets/tmd/
Global Tide Model (GOT) solutions provided by Richard Ray at GSFC
or Finite Element Solution (FES) models provided by AVISO
INPUTS:
x: x-coordinates in projection EPSG
y: y-coordinates in projection EPSG
OPTIONS:
DIRECTORY: working data directory for tide models
MODEL: Tide model to use
ATLAS_FORMAT: ATLAS tide model format (OTIS, netcdf)
GZIP: Tide model files are gzip compressed
DEFINITION_FILE: Tide model definition file for use
EPSG: input coordinate system
default: 3031 Polar Stereographic South, WGS84
METHOD: interpolation method
bilinear: quick bilinear interpolation
spline: scipy bivariate spline interpolation
linear, nearest: scipy regular grid interpolations
OUTPUTS:
valid: array describing if input coordinate is within model domain
PYTHON DEPENDENCIES:
numpy: Scientific Computing Tools For Python
https://numpy.org
https://numpy.org/doc/stable/user/numpy-for-matlab-users.html
scipy: Scientific Tools for Python
https://docs.scipy.org/doc/
netCDF4: Python interface to the netCDF C library
https://unidata.github.io/netcdf4-python/netCDF4/index.html
pyproj: Python interface to PROJ library
https://pypi.org/project/pyproj/
PROGRAM DEPENDENCIES:
crs.py: Coordinate Reference System (CRS) routines
io/model.py: retrieves tide model parameters for named tide models
io/OTIS.py: extract tidal harmonic constants from OTIS tide models
io/ATLAS.py: extract tidal harmonic constants from ATLAS netcdf models
io/GOT.py: extract tidal harmonic constants from GSFC GOT models
io/FES.py: extract tidal harmonic constants from FES tide models
interpolate.py: interpolation routines for spatial data
UPDATE HISTORY:
Updated 04/2024: use wrapper to importlib for optional dependencies
Updated 12/2023: use new crs class for coordinate reprojection
Updated 08/2023: changed ESR netCDF4 format to TMD3 format
Updated 04/2023: using pathlib to define and expand paths
Updated 03/2023: add basic variable typing to function inputs
Updated 12/2022: refactored tide read programs under io
refactored bilinear interpolation routine
Updated 11/2022: place some imports within try/except statements
use f-strings for formatting verbose or ascii output
Updated 05/2022: added ESR netCDF4 formats to list of model types
updated keyword arguments to read tide model programs
Updated 04/2022: updated docstrings to numpy documentation format
Updated 09/2021: refactor to use model class for files and attributes
Updated 07/2021: added check that tide model directory is accessible
Updated 06/2021: add try/except for input projection strings
Written 05/2021
"""
from __future__ import print_function, annotations
import logging
import pathlib
import numpy as np
import scipy.interpolate
import pyTMD.crs
import pyTMD.io
import pyTMD.io.model
import pyTMD.interpolate
import pyTMD.utilities
# attempt imports
pyproj = pyTMD.utilities.import_dependency('pyproj')
# PURPOSE: compute tides at points and times using tide model algorithms
def check_points(x: np.ndarray, y: np.ndarray,
DIRECTORY: str | pathlib.Path | None = None,
MODEL: str | None = None,
ATLAS_FORMAT: str = 'netcdf',
GZIP: bool = False,
DEFINITION_FILE: str | pathlib.Path | None = None,
EPSG: str | int = 3031,
METHOD: str = 'spline'
):
"""
Check if points are within a tide model domain
Parameters
----------
x: np.ndarray
x-coordinates in projection EPSG
y: np.ndarray
y-coordinates in projection EPSG
DIRECTORY: str or NoneType, default None
working data directory for tide models
MODEL: str or NoneType, default None
Tide model to use
ATLAS_FORMAT: str, default 'netcdf'
ATLAS tide model format
- ``'OTIS'``
- ``'netcdf'``
GZIP: bool, default False
Tide model files are gzip compressed
DEFINITION_FILE: str or NoneType, default None
Tide model definition file for use
EPSG: str or int, default: 3031 (Polar Stereographic South, WGS84)
Input coordinate system
METHOD: str, default 'spline'
interpolation method
- ```bilinear```: quick bilinear interpolation
- ```spline```: scipy bivariate spline interpolation
- ```linear```, ```nearest```: scipy regular grid interpolations
Returns
-------
valid: bool
array describing if input coordinate is within model domain
"""
# check that tide directory is accessible
if DIRECTORY is not None:
DIRECTORY = pathlib.Path(DIRECTORY).expanduser()
if not DIRECTORY.exists():
raise FileNotFoundError("Invalid tide directory")
# get parameters for tide model
if DEFINITION_FILE is not None:
model = pyTMD.io.model(DIRECTORY).from_file(
pathlib.Path(DEFINITION_FILE).expanduser())
else:
model = pyTMD.io.model(DIRECTORY, format=ATLAS_FORMAT,
compressed=GZIP).elevation(MODEL)
# input shape of data
idim = np.shape(x)
# converting x,y from input coordinate reference system
crs1 = pyTMD.crs().from_input(EPSG)
crs2 = pyproj.CRS.from_epsg(4326)
transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True)
lon, lat = transformer.transform(
np.atleast_1d(x).flatten(), np.atleast_1d(y).flatten()
)
# read tidal constants and interpolate to grid points
if model.format in ('OTIS','ATLAS','TMD3'):
# if reading a single OTIS solution
xi, yi, hz, mz, iob, dt = pyTMD.io.OTIS.read_otis_grid(
pathlib.Path(model.grid_file).expanduser())
# invert model mask
mz = np.logical_not(mz)
# adjust dimensions of input coordinates to be iterable
# run wrapper function to convert coordinate systems of input lat/lon
X, Y = pyTMD.crs().convert(lon, lat, model.projection, 'F')
elif (model.format == 'netcdf'):
# if reading a netCDF OTIS atlas solution
xi, yi, hz = pyTMD.io.ATLAS.read_netcdf_grid(
pathlib.Path(model.grid_file).expanduser(),
compressed=model.compressed, type=model.type)
# copy bathymetry mask
mz = np.copy(hz.mask)
# copy latitude and longitude and adjust longitudes
X,Y = np.copy([lon,lat]).astype(np.float64)
lt0, = np.nonzero(X < 0)
X[lt0] += 360.0
elif (model.format == 'GOT'):
# if reading a NASA GOT solution
hc, xi, yi, c = pyTMD.io.GOT.read_ascii_file(
pathlib.Path(model.model_file[0]).expanduser(),
compressed=model.compressed)
# copy tidal constituent mask
mz = np.copy(hc.mask)
# copy latitude and longitude and adjust longitudes
X, Y = np.copy([lon,lat]).astype(np.float64)
lt0, = np.nonzero(X < 0)
X[lt0] += 360.0
elif (model.format == 'FES'):
# if reading a FES netCDF solution
hc, xi, yi = pyTMD.io.FES.read_netcdf_file(
pathlib.Path(model.model_file[0]).expanduser(),
compressed=model.compressed, type=model.type,
version=model.version)
# copy tidal constituent mask
mz = np.copy(hc.mask)
# copy latitude and longitude and adjust longitudes
X, Y = np.copy([lon,lat]).astype(np.float64)
lt0, = np.nonzero(X < 0)
X[lt0] += 360.0
# interpolate masks
if (METHOD == 'bilinear'):
# replace invalid values with nan
mz1 = pyTMD.interpolate.bilinear(xi, yi, mz, X, Y)
mask = np.floor(mz1).astype(mz.dtype)
elif (METHOD == 'spline'):
f1 = scipy.interpolate.RectBivariateSpline(xi, yi, mz.T,
kx=1, ky=1)
mask = np.floor(f1.ev(X, Y)).astype(mz.dtype)
else:
# use scipy regular grid to interpolate values
r1 = scipy.interpolate.RegularGridInterpolator((yi, xi), mz,
method=METHOD, bounds_error=False, fill_value=1)
mask = np.floor(r1.__call__(np.c_[y, x])).astype(mz.dtype)
# reshape to original dimensions
valid = np.logical_not(mask).reshape(idim).astype(mz.dtype)
# replace points outside model domain with invalid
valid &= (X >= xi.min()) & (X <= xi.max())
valid &= (Y >= yi.min()) & (Y <= yi.max())
# return the valid mask
return valid