pyTMD
pyTMD copied to clipboard
tsutterley
Interpolation `METHOD="spline"` fails for TPXO8
Issue
I can easily model tides for TPXO8 using "bilinear" interpolation:
import numpy as np
import pandas as pd
from pyTMD import compute_tide_corrections
compute_tide_corrections(
x=np.linspace(155, 160, 10),
y=np.linspace(-30, -40, 10),
delta_time=pd.date_range("2020-01", "2020-02", periods=10),
DIRECTORY="/home/jovyan/tide_models_clipped",
MODEL="TPXO8-atlas",
EPSG=4326,
TYPE="drift",
TIME="datetime",
METHOD="bilinear",
)
However, if I try the same code with "spline" interpolation (the default), this fails with ValueError: x dimension of z must have same number of elements as x:
compute_tide_corrections(
x=np.linspace(155, 160, 10),
y=np.linspace(-30, -40, 10),
delta_time=pd.date_range("2020-01", "2020-02", periods=10),
DIRECTORY="/home/jovyan/tide_models_clipped",
MODEL="TPXO8-atlas",
EPSG=4326,
TYPE="drift",
TIME="datetime",
METHOD="spline",
)
Other models work fine with "spline" interpolation, e.g. FES2014:
compute_tide_corrections(
x=np.linspace(155, 160, 10),
y=np.linspace(-30, -40, 10),
delta_time=pd.date_range("2020-01", "2020-02", periods=10),
DIRECTORY="/home/jovyan/tide_models_clipped",
MODEL="FES2014",
EPSG=4326,
TYPE="drift",
TIME="datetime",
METHOD="spline",
)
Expected outcome
All supported tide models can be run with the default "spline" interpolation method without raising errors.
Details
I'm using pyTMD=2.0.7 and scipy=1.10.1. Full error:
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[23], line 1
----> 1 compute_tide_corrections(
2 x=np.linspace(155, 160, 10),
3 y=np.linspace(-30, -40, 10),
4 delta_time=pd.date_range("2020-01", "2020-02", periods=10),
5 DIRECTORY="/home/jovyan/tide_models_clipped",
6 MODEL="TPXO8-atlas",
7 EPSG=4326,
8 TYPE="drift",
9 TIME="datetime",
10 METHOD="spline",
11 )
File /env/lib/python3.8/site-packages/pyTMD/compute_tide_corrections.py:318, in compute_tide_corrections(x, y, delta_time, DIRECTORY, MODEL, ATLAS_FORMAT, GZIP, DEFINITION_FILE, EPSG, EPOCH, TYPE, TIME, METHOD, EXTRAPOLATE, CUTOFF, APPLY_FLEXURE, FILL_VALUE, **kwargs)
316 # read tidal constants and interpolate to grid points
317 if model.format in ('OTIS','ATLAS','TMD3'):
--> 318 amp,ph,D,c = pyTMD.io.OTIS.extract_constants(lon, lat, model.grid_file,
319 model.model_file, model.projection, type=model.type,
320 method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF,
321 grid=model.format, apply_flexure=APPLY_FLEXURE)
322 # use delta time at 2000.0 to match TMD outputs
323 deltat = np.zeros((nt), dtype=np.float64)
File /env/lib/python3.8/site-packages/pyTMD/io/OTIS.py:410, in extract_constants(ilon, ilat, grid_file, model_file, EPSG, **kwargs)
407 hci.data[hci.mask] = hci.fill_value
408 elif (kwargs['method'] == 'spline'):
409 # use scipy bivariate splines to interpolate values
--> 410 hci = pyTMD.interpolate.spline(xi, yi, hc, x, y,
411 dtype=hc.dtype,
412 reducer=np.ceil,
413 kx=1, ky=1)
414 # replace zero values with fill_value
415 hci.mask |= D.mask
File /env/lib/python3.8/site-packages/pyTMD/interpolate.py:172, in spline(ilon, ilat, idata, lon, lat, fill_value, dtype, reducer, **kwargs)
170 # construct splines for input data and mask
171 if np.iscomplexobj(idata):
--> 172 s1 = scipy.interpolate.RectBivariateSpline(ilon, ilat,
173 idata.data.real.T, **kwargs)
174 s2 = scipy.interpolate.RectBivariateSpline(ilon, ilat,
175 idata.data.imag.T, **kwargs)
176 s3 = scipy.interpolate.RectBivariateSpline(ilon, ilat,
177 idata.mask.T, **kwargs)
File /env/lib/python3.8/site-packages/scipy/interpolate/_fitpack2.py:1494, in RectBivariateSpline.__init__(self, x, y, z, bbox, kx, ky, s)
1492 raise ValueError('y must be strictly increasing')
1493 if not x.size == z.shape[0]:
-> 1494 raise ValueError('x dimension of z must have same number of '
1495 'elements as x')
1496 if not y.size == z.shape[1]:
1497 raise ValueError('y dimension of z must have same number of '
1498 'elements as y')
ValueError: x dimension of z must have same number of elements as x
"linear" also seems to fail with a different error:
compute_tide_corrections(
x=np.linspace(155, 160, 10),
y=np.linspace(-30, -40, 10),
delta_time=pd.date_range("2020-01", "2020-02", periods=10),
DIRECTORY="/home/jovyan/tide_models_clipped",
MODEL="TPXO8-atlas",
EPSG=4326,
TYPE="drift",
TIME="datetime",
METHOD="linear",
)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[26], line 1
----> 1 compute_tide_corrections(
2 x=np.linspace(155, 160, 10),
3 y=np.linspace(-30, -40, 10),
4 delta_time=pd.date_range("2020-01", "2020-02", periods=10),
5 DIRECTORY="/home/jovyan/tide_models_clipped",
6 MODEL="TPXO8-atlas",
7 EPSG=4326,
8 TYPE="drift",
9 TIME="datetime",
10 METHOD="linear",
11 )
File /env/lib/python3.8/site-packages/pyTMD/compute_tide_corrections.py:318, in compute_tide_corrections(x, y, delta_time, DIRECTORY, MODEL, ATLAS_FORMAT, GZIP, DEFINITION_FILE, EPSG, EPOCH, TYPE, TIME, METHOD, EXTRAPOLATE, CUTOFF, APPLY_FLEXURE, FILL_VALUE, **kwargs)
316 # read tidal constants and interpolate to grid points
317 if model.format in ('OTIS','ATLAS','TMD3'):
--> 318 amp,ph,D,c = pyTMD.io.OTIS.extract_constants(lon, lat, model.grid_file,
319 model.model_file, model.projection, type=model.type,
320 method=METHOD, extrapolate=EXTRAPOLATE, cutoff=CUTOFF,
321 grid=model.format, apply_flexure=APPLY_FLEXURE)
322 # use delta time at 2000.0 to match TMD outputs
323 deltat = np.zeros((nt), dtype=np.float64)
File /env/lib/python3.8/site-packages/pyTMD/io/OTIS.py:419, in extract_constants(ilon, ilat, grid_file, model_file, EPSG, **kwargs)
416 hci.data[hci.mask] = hci.fill_value
417 else:
418 # use scipy regular grid to interpolate values
--> 419 hci = pyTMD.interpolate.regulargrid(xi, yi, hc, x, y,
420 fill_value=hc.fill_value,
421 dtype=hc.dtype,
422 method=kwargs['method'],
423 reducer=np.ceil,
424 bounds_error=False)
425 # replace invalid values with fill_value
426 hci.mask = (hci.data == hci.fill_value) | D.mask
File /env/lib/python3.8/site-packages/pyTMD/interpolate.py:259, in regulargrid(ilon, ilat, idata, lon, lat, fill_value, dtype, reducer, **kwargs)
257 data.mask = np.ones((npts), dtype=bool)
258 # use scipy regular grid to interpolate values for a given method
--> 259 r1 = scipy.interpolate.RegularGridInterpolator((ilat, ilon),
260 idata.data, fill_value=fill_value, **kwargs)
261 r2 = scipy.interpolate.RegularGridInterpolator((ilat, ilon),
262 idata.mask, fill_value=1, **kwargs)
263 # evaluate the interpolator at input coordinates
File /env/lib/python3.8/site-packages/scipy/interpolate/_rgi.py:242, in RegularGridInterpolator.__init__(self, points, values, method, bounds_error, fill_value)
240 self.grid, self._descending_dimensions = _check_points(points)
241 self.values = self._check_values(values)
--> 242 self._check_dimensionality(self.grid, self.values)
243 self.fill_value = self._check_fill_value(self.values, fill_value)
244 if self._descending_dimensions:
File /env/lib/python3.8/site-packages/scipy/interpolate/_rgi.py:248, in RegularGridInterpolator._check_dimensionality(self, grid, values)
247 def _check_dimensionality(self, grid, values):
--> 248 _check_dimensionality(grid, values)
File /env/lib/python3.8/site-packages/scipy/interpolate/_rgi.py:45, in _check_dimensionality(points, values)
42 raise ValueError("The points in dimension %d must be "
43 "1-dimensional" % i)
44 if not values.shape[i] == len(p):
---> 45 raise ValueError("There are %d points and %d values in "
46 "dimension %d" % (len(p), values.shape[i], i))
ValueError: There are 10800 points and 10802 values in dimension 1
Thanks @robbibt for finding this bug and pointing it out! Should be fixed in 2.0.8. Let me know if you have any more trouble!
Amazing, thanks so much @tsutterley for the quick fix (as always). I'll test it out and let you know.
It works nicely now on 2.0.8 - awesome! 🚀
A related TPXO8 question: for some of our other models, we've found we can get a really big performance boost by clipping our tide model NetCDF files to Australia before passing them to pyTMD. But for TPXO8, we're currently using what I believe is the "compact" versions of the model, e.g.:
These aren't something I can easily clip, so I was wondering if pyTMD also supports the NetCDF versions of TPXO8's model files referenced on the OSU website here? https://www.tpxo.net/tpxo-products-and-registration
Excellent!
You're right. What's hard coded into the model class for TPXO8-atlas is the compact format. A model definition file should be able to be used with subsetted TPXO8-atlas netCDFs. I wrote the following model definition for someone using the (global) netCDF format. Check that the scaling is right for the units (they said it was decimeters).
format netcdf
name TPXO8-atlas-v1
model_file hf.*_tpxo8_atlas_30c_v1.nc
grid_file grid_tpxo8_atlas_30c_v1.nc
type z
variable tide_ocean
version v1
scale 0.1
compressed False
reference https://www.tpxo.net/global/tpxo8-atlas
p.s. the above definition "file" uses the glob functionality to find files (so set --directory to point to where the files are located).
Neat! That's working nicely for me - I had to make two small changes:
- Update
grid_filetogrid_tpxo8atlas_30_v1.nc(without an underscore and the "c" after 30); the grid file has a different name format to the other NetCDFs supplied by OSU. - Change
scaleto0.001to get outputs in metre units
Am now able to model tides for all pyTMD supported models in a fraction of the time - very exciting!
Thanks @robbibt! I'll make those change to my version of the definition file. I'm still thinking about the best way for users to enter these arguments, while still allowing for ease of use. The definition files are a bit of a kludgy solution.
Awesome! That's super exciting! :rocket: