add nearest-neighbor extrapolation (#35)

update documentation for extrapolation

update get_hash to accept bytesio objects

use kd-trees for nearest neighbors

reduce nearest to close points

only extrapolate in tests

add more test functions for utilities

remove 3.5 from tests

reduce codecov notifications to after 6

remove some instances of trailing whitespace

.github/workflows/python-package.yml

@@ -15,7 +15,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-20.04, macos-latest]
-        python-version: [3.5, 3.6, 3.7, 3.8]
+        python-version: [3.6, 3.7, 3.8]
     env:
       OS: ${{ matrix.os }}
       PYTHON: ${{ matrix.python-version }}
@@ -73,7 +73,7 @@ jobs:
         ref: inferminor
     - name: Test with pytest
       run: |
-        pytest --cov=./ --cov-report=xml \
+        pytest --verbose --cov=./ --cov-report=xml \
             --username=${{ secrets.EARTHDATA_USERNAME }} \
             --password=${{ secrets.EARTHDATA_PASSWORD }} \
             --aws-access=${{ secrets.AWS_ACCESS_KEY_ID }} \

.github/workflows/python-request.yml

@@ -13,7 +13,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-20.04, macos-latest]
-        python-version: [3.5, 3.6, 3.7, 3.8]
+        python-version: [3.6, 3.7, 3.8]
     env:
       OS: ${{ matrix.os }}
       PYTHON: ${{ matrix.python-version }}
@@ -71,7 +71,7 @@ jobs:
         ref: inferminor
     - name: Test with pytest
       run: |
-        pytest --cov=./ --cov-report=xml \
+        pytest --verbose --cov=./ --cov-report=xml \
             --username=${{ secrets.EARTHDATA_USERNAME }} \
             --password=${{ secrets.EARTHDATA_PASSWORD }} \
             --aws-access=${{ secrets.AWS_ACCESS_KEY_ID }} \

codecov.yml

@@ -1,6 +1,6 @@
 codecov:
     notify:
-        after_n_builds: 8
+        after_n_builds: 6
 
 coverage:
   status:

doc/source/getting_started/Getting-Started.md

@@ -101,3 +101,9 @@ The default coordinate system in pyTMD is WGS84 geodetic coordinates in latitude
 pyTMD uses [pyproj](https://pypi.org/project/pyproj/) to convert from different coordinate systems and datums.
 Some regional tide models are projected in a different coordinate system.
 For these cases, pyTMD will [convert from latitude and longitude to the model coordinate system](https://github.com/tsutterley/pyTMD/blob/main/pyTMD/convert_ll_xy.py).
+
+#### Interpolation
+For converting from model coordinates, pyTMD uses spatial interpolation routines from [scipy](https://docs.scipy.org/doc/scipy/reference/interpolate.html) along with a built-in [`'bilinear'`](https://github.com/tsutterley/pyTMD/blob/main/pyTMD/bilinear_interp.py) interpolation routine.
+The default interpolator uses a [biharmonic `'spline'`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RectBivariateSpline.html) function to interpolate from the model coordinate system to the output coordinates.
+There are options to use `'nearest'` and `'linear'` interpolators with the [regular grid](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RegularGridInterpolator.html) function.
+For coastal or near-grounded points, the model can be extrapolated using a [nearest-neighbor](https://github.com/tsutterley/pyTMD/blob/main/pyTMD/nearest_extrap.py) routine.

doc/source/index.rst

@@ -56,6 +56,7 @@ conventions for calculating radial pole tide displacements.
     user_guide/infer_minor_corrections.md
     user_guide/load_constituent.md
     user_guide/load_nodal_corrections.md
+    user_guide/nearest_extrap.md
     user_guide/output_otis_tides.md
     user_guide/predict_tidal_ts.md
     user_guide/predict_tide_drift.md

doc/source/user_guide/bilinear_interp.md

@@ -20,7 +20,6 @@ data = bilinear_interp(ilon,ilat,idata,lon,lat)
 #### Options
  - `fill_value`: invalid value
  - `dtype`: output data type
- - `extrapolate`: extrapolate points
 
 #### Outputs
  - `data`: interpolated data

doc/source/user_guide/compute_tidal_currents.md

@@ -51,5 +51,6 @@ python compute_tidal_currents.py --directory <path_to_directory> --tide <model>
      * `'linear'`
      * `'nearest'`
      * `'bilinear'`
+ - `-E X`, `--extrapolate X`: Extrapolate with nearest-neighbors
  - `-V`, `--verbose`: Verbose output of processing run
  - `-M X`, `--mode X`: Permission mode of output file

doc/source/user_guide/compute_tidal_elevations.md

@@ -63,5 +63,6 @@ python compute_tidal_elevations.py --directory <path_to_directory> --tide <model
      * `'linear'`
      * `'nearest'`
      * `'bilinear'`
+ - `-E X`, `--extrapolate X`: Extrapolate with nearest-neighbors
  - `-V`, `--verbose`: Verbose output of processing run
  - `-M X`, `--mode X`: Permission mode of output file

doc/source/user_guide/compute_tide_corrections.md

@@ -37,6 +37,7 @@ tide = compute_tide_corrections(x, y, delta_time, DIRECTORY=DIRECTORY,
      * `bilinear`: quick bilinear interpolation
      * `spline`: scipy bivariate spline interpolation (default)
      * `linear`, `nearest`: scipy regular grid interpolations
+ - `EXTRAPOLATE`: Extrapolate with nearest-neighbors
  - `FILL_VALUE`: output invalid value
 
 #### Outputs

doc/source/user_guide/compute_tides_ICESat2_ATL03.md

@@ -46,5 +46,6 @@ python compute_tides_ICESat2_ATL03.py --directory <path_to_directory> --tide <mo
      * `'linear'`
      * `'nearest'`
      * `'bilinear'`
+ - `-E X`, `--extrapolate X`: Extrapolate with nearest-neighbors
  - `-M X`, `--mode X`: Permission mode of output file
  - `-V`, `--verbose`: Output information about each created file

doc/source/user_guide/compute_tides_ICESat2_ATL06.md

@@ -44,5 +44,6 @@ python compute_tides_ICESat2_ATL06.py --directory <path_to_directory> --tide <mo
      * `'linear'`
      * `'nearest'`
      * `'bilinear'`
+ - `-E X`, `--extrapolate X`: Extrapolate with nearest-neighbors
  - `-M X`, `--mode X`: Permission mode of output file
  - `-V`, `--verbose`: Output information about each created file

doc/source/user_guide/compute_tides_ICESat2_ATL07.md

@@ -44,5 +44,6 @@ python compute_tides_ICESat2_ATL07.py --directory <path_to_directory> --tide <mo
      * `'linear'`
      * `'nearest'`
      * `'bilinear'`
+ - `-E X`, `--extrapolate X`: Extrapolate with nearest-neighbors
  - `-M X`, `--mode X`: Permission mode of output file
  - `-V`, `--verbose`: Output information about each created file

doc/source/user_guide/compute_tides_ICESat2_ATL11.md

@@ -44,5 +44,6 @@ python compute_tides_ICESat2_ATL11.py --directory <path_to_directory> --tide <mo
      * `'linear'`
      * `'nearest'`
      * `'bilinear'`
+ - `-E X`, `--extrapolate X`: Extrapolate with nearest-neighbors
  - `-M X`, `--mode X`: Permission mode of output file
  - `-V`, `--verbose`: Output information about each created file

doc/source/user_guide/compute_tides_ICESat2_ATL12.md

@@ -44,5 +44,6 @@ python compute_tides_ICESat2_ATL12.py --directory <path_to_directory> --tide <mo
      * `'linear'`
      * `'nearest'`
      * `'bilinear'`
+ - `-E X`, `--extrapolate X`: Extrapolate with nearest-neighbors
  - `-M X`, `--mode X`: Permission mode of output file
  - `-V`, `--verbose`: Output information about each created file

doc/source/user_guide/compute_tides_ICESat_GLA12.md

@@ -44,5 +44,6 @@ python compute_tides_ICESat_GLA12.py --directory <path_to_directory> --tide <mod
      * `'linear'`
      * `'nearest'`
      * `'bilinear'`
+ - `-E X`, `--extrapolate X`: Extrapolate with nearest-neighbors
  - `-M X`, `--mode X`: Permission mode of output file
  - `-V`, `--verbose`: Output information about each created file

doc/source/user_guide/compute_tides_icebridge_data.md

@@ -44,5 +44,6 @@ python compute_tides_icebridge_data.py --directory <path_to_directory> --tide <m
      * `'linear'`
      * `'nearest'`
      * `'bilinear'`
+ - `-E X`, `--extrapolate X`: Extrapolate with nearest-neighbors
  - `-M X`, `--mode X`: Permission mode of output file
  - `-V`, `--verbose`: Output information about each created file

doc/source/user_guide/nearest_extrap.md

@@ -0,0 +1,27 @@
+nearest_extrap.py
+=================
+
+ - Uses kd-trees for nearest-neighbor extrapolation of valid model data
+
+#### Calling Sequence
+```python
+from pyTMD.nearest_extrap import nearest_extrap
+data = nearest_extrap(ilon,ilat,idata,lon,lat)
+```
+[Source code](https://github.com/tsutterley/pyTMD/blob/main/pyTMD/nearest_extrap.py)
+
+#### Inputs
+ 1. `ilon`: longitude of tidal model
+ 2. `ilat`: latitude of tidal model
+ 3. `idata`: tide model data
+ 4. `lon`: output longitude
+ 5. `lat`: output latitude
+
+#### Options
+ - `fill_value`: invalid value
+ - `dtype`: output data type
+ - `cutoff`: return only neighbors within distance in kilometers
+ - `EPSG`: projection of tide model data
+
+#### Outputs
+ - `data`: interpolated data

doc/source/user_guide/utilities.rst

@@ -29,11 +29,11 @@ General Methods
 
 .. method:: pyTMD.utilities.get_hash(local)
 
-    Get the MD5 hash value from a local file
+    Get the MD5 hash value from a local file or BytesIO object
 
     Arguments:
 
-        `local`: path to file
+        `local`: BytesIO object or path to file
 
 
 .. method:: pyTMD.utilities.url_split(s)

pyTMD/__init__.py

@@ -26,6 +26,7 @@
 from pyTMD.infer_minor_corrections import infer_minor_corrections
 from pyTMD.load_constituent import load_constituent
 from pyTMD.load_nodal_corrections import load_nodal_corrections
+from pyTMD.nearest_extrap import nearest_extrap
 from pyTMD.output_otis_tides import output_otis_grid
 from pyTMD.read_tide_model import extract_tidal_constants
 from pyTMD.read_netcdf_model import extract_netcdf_constants

pyTMD/bilinear_interp.py

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 u"""
-bilinear_interp.py (08/2020)
+bilinear_interp.py (12/2020)
 Bilinear interpolation of input data to output coordinates
 
 CALLING SEQUENCE:
@@ -16,7 +16,6 @@
 OPTIONS:
     fill_value: invalid value
     dtype: output data type
-    extrapolate: extrapolate points
 
 OUTPUT:
     data: interpolated data
@@ -27,6 +26,7 @@
         https://numpy.org/doc/stable/user/numpy-for-matlab-users.html
 
 UPDATE HISTORY:
+    Updated 12/2020: using numpy isclose to check corner points
     Updated 08/2020: check that output coordinates are within bounds
         allow small extrapolations if individual grid cells are invalid
     Updated 07/2020: split into separate function
@@ -37,7 +37,7 @@
 
 #-- PURPOSE: bilinear interpolation of input data to output data
 def bilinear_interp(ilon,ilat,idata,lon,lat,fill_value=np.nan,
-    dtype=np.float,extrapolate=False):
+    dtype=np.float):
     """
     Bilinear interpolation of input data to output coordinates
 
@@ -53,15 +53,11 @@ def bilinear_interp(ilon,ilat,idata,lon,lat,fill_value=np.nan,
     -----------------
     fill_value: invalid value
     dtype: output data type
-    extrapolate: extrapolate points
 
     Returns
     -------
     data: interpolated data
     """
-    #-- grid step size of tide model
-    dlon = np.abs(ilon[1] - ilon[0])
-    dlat = np.abs(ilat[1] - ilat[0])
     #-- find valid points (within bounds)
     valid, = np.nonzero((lon >= ilon.min()) & (lon <= ilon.max()) &
         (lat > ilat.min()) & (lat < ilat.max()))
@@ -88,19 +84,19 @@ def bilinear_interp(ilon,ilat,idata,lon,lat,fill_value=np.nan,
             IM.mask[j], = idata.mask[YI,XI]
             WM[j], = np.abs(lon[i]-ilon[XI])*np.abs(lat[i]-ilat[YI])
         #-- if on corner value: use exact
-        if ((lat[i] == ilat[iy]) & (lon[i] == ilon[ix])):
+        if (np.isclose(lat[i],ilat[iy]) & np.isclose(lon[i],ilon[ix])):
             data.data[i] = idata.data[iy,ix]
             data.mask[i] = idata.mask[iy,ix]
-        elif ((lat[i] == ilat[iy+1]) & (lon[i] == ilon[ix])):
+        elif (np.isclose(lat[i],ilat[iy+1]) & np.isclose(lon[i],ilon[ix])):
             data.data[i] = idata.data[iy+1,ix]
             data.mask[i] = idata.mask[iy+1,ix]
-        elif ((lat[i] == ilat[iy]) & (lon[i] == ilon[ix+1])):
+        elif (np.isclose(lat[i],ilat[iy]) & np.isclose(lon[i],ilon[ix+1])):
             data.data[i] = idata.data[iy,ix+1]
             data.mask[i] = idata.mask[iy,ix+1]
-        elif ((lat[i] == ilat[iy+1]) & (lon[i] == ilon[ix+1])):
+        elif (np.isclose(lat[i],ilat[iy+1]) & np.isclose(lon[i],ilon[ix+1])):
             data.data[i] = idata.data[iy+1,ix+1]
             data.mask[i] = idata.mask[iy+1,ix+1]
-        elif np.all(np.isfinite(IM) & (~IM.mask)) or extrapolate:
+        elif np.any(np.isfinite(IM) & (~IM.mask)):
             #-- find valid indices for data summation and weight matrix
             ii, = np.nonzero(np.isfinite(IM) & (~IM.mask))
             #-- calculate interpolated value for i

pyTMD/compute_tide_corrections.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 compute_tide_corrections.py
-Written by Tyler Sutterley (11/2020)
+Written by Tyler Sutterley (12/2020)
 Calculates tidal elevations for correcting elevation or imagery data
 
 Uses OTIS format tidal solutions provided by Ohio State University and ESR
@@ -62,8 +62,11 @@
     read_netcdf_model.py: extract tidal harmonic constants from netcdf models
     read_GOT_model.py: extract tidal harmonic constants from GSFC GOT models
     read_FES_model.py: extract tidal harmonic constants from FES tide models
+    bilinear_interp.py: bilinear interpolation of data to coordinates
+    nearest_extrap.py: nearest-neighbor extrapolation of data to coordinates
 
 UPDATE HISTORY:
+    Updated 12/2020: added valid data extrapolation with nearest_extrap
     Updated 11/2020: added model constituents from TPXO9-atlas-v3
     Updated 08/2020: using builtin time operations.
         calculate difference in leap seconds from start of epoch
@@ -92,7 +95,7 @@
 #-- PURPOSE: compute tides at points and times using tide model algorithms
 def compute_tide_corrections(x, y, delta_time, DIRECTORY=None, MODEL=None,
     EPSG=3031, EPOCH=(2000,1,1,0,0,0), TYPE='drift', TIME='UTC',
-    METHOD='spline', FILL_VALUE=np.nan):
+    METHOD='spline', EXTRAPOLATE=False, FILL_VALUE=np.nan):
     """
     Compute tides at points and times using tidal harmonics
 
@@ -340,21 +343,24 @@ def compute_tide_corrections(x, y, delta_time, DIRECTORY=None, MODEL=None,
     #-- read tidal constants and interpolate to grid points
     if model_format in ('OTIS','ATLAS'):
         amp,ph,D,c = extract_tidal_constants(lon, lat, grid_file, model_file,
-            model_EPSG, TYPE=model_type, METHOD=METHOD, GRID=model_format)
+            model_EPSG, TYPE=model_type, METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE,
+            GRID=model_format)
         deltat = np.zeros_like(t)
     elif (model_format == 'netcdf'):
         amp,ph,D,c = extract_netcdf_constants(lon, lat, model_directory,
-            grid_file, model_files, TYPE=model_type, METHOD=METHOD, SCALE=SCALE)
+            grid_file, model_files, TYPE=model_type, METHOD=METHOD, 
+            EXTRAPOLATE=EXTRAPOLATE, SCALE=SCALE)
         deltat = np.zeros_like(t)
     elif (model_format == 'GOT'):
         amp,ph = extract_GOT_constants(lon, lat, model_directory, model_files,
-            METHOD=METHOD, SCALE=SCALE)
+            METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE, SCALE=SCALE)
         #-- interpolate delta times from calendar dates to tide time
         delta_file = pyTMD.utilities.get_data_path(['data','merged_deltat.data'])
         deltat = calc_delta_time(delta_file, t)
     elif (model_format == 'FES'):
         amp,ph = extract_FES_constants(lon, lat, model_directory, model_files,
-            TYPE=model_type, VERSION=MODEL, METHOD=METHOD, SCALE=SCALE)
+            TYPE=model_type, VERSION=MODEL, METHOD=METHOD,
+            EXTRAPOLATE=EXTRAPOLATE, SCALE=SCALE)
         #-- interpolate delta times from calendar dates to tide time
         delta_file = pyTMD.utilities.get_data_path(['data','merged_deltat.data'])
         deltat = calc_delta_time(delta_file, t)