Interpolation METHOD="spline" fails for TPXO8

[robbibt]

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

[robbibt]

"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

[tsutterley]

That's not good! I'll work on this today.

[tsutterley]

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!

[robbibt]

Amazing, thanks so much @tsutterley for the quick fix (as always). I'll test it out and let you know.

[robbibt]

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

[tsutterley]

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).

[robbibt]

Neat! That's working nicely for me - I had to make two small changes:

• Update grid_file to grid_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 scale to 0.001 to get outputs in metre units

Am now able to model tides for all pyTMD supported models in a fraction of the time - very exciting!

[tsutterley]

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! 🚀