From 11120639a8895e5b3cc3067f991ee39047f59194 Mon Sep 17 00:00:00 2001 From: Tyler Sutterley Date: Sun, 20 Dec 2020 20:45:03 -0800 Subject: [PATCH] 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 | 4 +- .github/workflows/python-request.yml | 4 +- codecov.yml | 2 +- doc/source/getting_started/Getting-Started.md | 6 + doc/source/index.rst | 1 + doc/source/user_guide/bilinear_interp.md | 1 - .../user_guide/compute_tidal_currents.md | 1 + .../user_guide/compute_tidal_elevations.md | 1 + .../user_guide/compute_tide_corrections.md | 1 + .../user_guide/compute_tides_ICESat2_ATL03.md | 1 + .../user_guide/compute_tides_ICESat2_ATL06.md | 1 + .../user_guide/compute_tides_ICESat2_ATL07.md | 1 + .../user_guide/compute_tides_ICESat2_ATL11.md | 1 + .../user_guide/compute_tides_ICESat2_ATL12.md | 1 + .../user_guide/compute_tides_ICESat_GLA12.md | 1 + .../compute_tides_icebridge_data.md | 1 + doc/source/user_guide/nearest_extrap.md | 27 ++++ doc/source/user_guide/utilities.rst | 4 +- pyTMD/__init__.py | 1 + pyTMD/bilinear_interp.py | 20 +-- pyTMD/compute_tide_corrections.py | 18 ++- pyTMD/nearest_extrap.py | 143 ++++++++++++++++++ pyTMD/read_FES_model.py | 22 ++- pyTMD/read_GOT_model.py | 21 ++- pyTMD/read_netcdf_model.py | 38 ++++- pyTMD/read_tide_model.py | 41 ++++- pyTMD/utilities.py | 14 +- scripts/compute_LPET_ICESat2_ATL11.py | 14 +- scripts/compute_LPET_ICESat_GLA12.py | 4 +- scripts/compute_LPT_ICESat_GLA12.py | 4 +- scripts/compute_OPT_ICESat_GLA12.py | 4 +- scripts/compute_tidal_currents.py | 23 ++- scripts/compute_tidal_elevations.py | 26 +++- scripts/compute_tides_ICESat2_ATL03.py | 26 +++- scripts/compute_tides_ICESat2_ATL06.py | 25 ++- scripts/compute_tides_ICESat2_ATL07.py | 25 ++- scripts/compute_tides_ICESat2_ATL11.py | 39 +++-- scripts/compute_tides_ICESat2_ATL12.py | 25 ++- scripts/compute_tides_ICESat_GLA12.py | 30 ++-- scripts/compute_tides_icebridge_data.py | 26 +++- setup.py | 1 - test/test_atlas_read.py | 8 +- test/test_perth3_read.py | 8 +- test/test_time.py | 5 + test/test_utilities.py | 33 ++++ 45 files changed, 560 insertions(+), 143 deletions(-) create mode 100644 doc/source/user_guide/nearest_extrap.md create mode 100644 pyTMD/nearest_extrap.py create mode 100644 test/test_utilities.py diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index b9ac9a2b..6051d865 100644 --- a/.github/workflows/python-package.yml +++ b/.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 }} \ diff --git a/.github/workflows/python-request.yml b/.github/workflows/python-request.yml index 70b712d6..5daafc74 100644 --- a/.github/workflows/python-request.yml +++ b/.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 }} \ diff --git a/codecov.yml b/codecov.yml index 320164f4..69abc9ab 100644 --- a/codecov.yml +++ b/codecov.yml @@ -1,6 +1,6 @@ codecov: notify: - after_n_builds: 8 + after_n_builds: 6 coverage: status: diff --git a/doc/source/getting_started/Getting-Started.md b/doc/source/getting_started/Getting-Started.md index 355207c9..fe0c02d6 100644 --- a/doc/source/getting_started/Getting-Started.md +++ b/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. diff --git a/doc/source/index.rst b/doc/source/index.rst index f9313416..6d4e477e 100644 --- a/doc/source/index.rst +++ b/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 diff --git a/doc/source/user_guide/bilinear_interp.md b/doc/source/user_guide/bilinear_interp.md index 131cad73..92e7c429 100644 --- a/doc/source/user_guide/bilinear_interp.md +++ b/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 diff --git a/doc/source/user_guide/compute_tidal_currents.md b/doc/source/user_guide/compute_tidal_currents.md index 7c233465..4f4131b9 100644 --- a/doc/source/user_guide/compute_tidal_currents.md +++ b/doc/source/user_guide/compute_tidal_currents.md @@ -51,5 +51,6 @@ python compute_tidal_currents.py --directory --tide * `'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 diff --git a/doc/source/user_guide/compute_tidal_elevations.md b/doc/source/user_guide/compute_tidal_elevations.md index 81591002..cd429baa 100644 --- a/doc/source/user_guide/compute_tidal_elevations.md +++ b/doc/source/user_guide/compute_tidal_elevations.md @@ -63,5 +63,6 @@ python compute_tidal_elevations.py --directory --tide --tide --tide --tide --tide --tide --tide --tide = 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 diff --git a/pyTMD/compute_tide_corrections.py b/pyTMD/compute_tide_corrections.py index 58e53b44..4f22c51e 100644 --- a/pyTMD/compute_tide_corrections.py +++ b/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) diff --git a/pyTMD/nearest_extrap.py b/pyTMD/nearest_extrap.py new file mode 100644 index 00000000..2e2b7fb1 --- /dev/null +++ b/pyTMD/nearest_extrap.py @@ -0,0 +1,143 @@ +#!/usr/bin/env python +u""" +nearest_extrap.py (12/2020) +Uses kd-trees for nearest-neighbor extrapolation of valid model data + +CALLING SEQUENCE: + data = nearest_extrap(ilon,ilat,idata,lon,lat) + +INPUTS: + ilon: longitude of tidal model + ilat: latitude of tidal model + idata: tide model data + lat: output latitude + lon: output longitude + +OPTIONS: + fill_value: invalid value + dtype: output data type + cutoff: return only neighbors within distance [km] + EPSG: projection of tide model data + +OUTPUT: + data: extrapolated data + +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/ + +UPDATE HISTORY: + Written 12/2020 +""" +import numpy as np +import scipy.spatial + +#-- PURPOSE: Nearest-neighbor extrapolation of valid data to output data +def nearest_extrap(ilon,ilat,idata,lon,lat,fill_value=np.nan, + dtype=np.float,cutoff=np.inf,EPSG='4326'): + """ + Nearest-neighbor extrapolation of valid model data + + Arguments + --------- + ilon: longitude of tidal model + ilat: latitude of tidal model + idata: tide model data + lat: output latitude + lon: output longitude + + Keyword arguments + ----------------- + fill_value: invalid value + dtype: output data type + cutoff: return only neighbors within distance [km] + EPSG: projection of tide model data + + Returns + ------- + data: interpolated data + """ + #-- grid step size of tide model + dlon = np.abs(ilon[1] - ilon[0]) + dlat = np.abs(ilat[1] - ilat[0]) + #-- extrapolate valid data values to data + npts = len(lon) + #-- allocate to output extrapolate data array + data = np.ma.zeros((npts),dtype=dtype,fill_value=fill_value) + data.mask = np.ones((npts),dtype=np.bool) + #-- initially set all data to fill value + data.data[:] = data.fill_value + #-- range of output points + xmin,xmax = (np.min(lon),np.max(lon)) + ymin,ymax = (np.min(lat),np.max(lat)) + + #-- calculate meshgrid of model coordinates + gridlon,gridlat = np.meshgrid(ilon,ilat) + #-- find where input grid is valid and close to output points + indy,indx = np.nonzero((~idata.mask) & np.isfinite(idata.data) & + (gridlon >= (xmin-2.0*dlon)) & (gridlon <= (xmax+2.0*dlon)) & + (gridlat >= (ymin-2.0*dlat)) & (gridlat <= (ymax+2.0*dlat))) + #-- flattened valid data array + iflat = idata.data[indy,indx] + + + #-- calculate coordinates for nearest-neighbors + if (EPSG == '4326'): + #-- valid grid latitude and longitude in radians + iphi = np.pi*gridlon[indy,indx]/180.0 + itheta = np.pi*gridlat[indy,indx]/180.0 + #-- WGS84 ellipsoid parameters + a_axis = 6378.1366 + f = 1.0/298.257223563 + ecc1 = np.sqrt((2.0*f - f**2)*a_axis**2)/a_axis + #-- calculate Cartesian coordinates of input grid + N = a_axis/np.sqrt(1.0-ecc1**2.0*np.sin(itheta)**2.0) + xflat = N*np.cos(itheta)*np.cos(iphi) + yflat = N*np.cos(itheta)*np.sin(iphi) + zflat = (N*(1.0-ecc1**2.0))*np.sin(itheta) + tree = scipy.spatial.cKDTree(np.c_[xflat,yflat,zflat]) + #-- calculate Cartesian coordinates of output coordinates + ns = a_axis/np.sqrt(1.0-ecc1**2.0*np.sin(np.pi*lat/180.0)**2.0) + xs = ns*np.cos(np.pi*lat/180.0)*np.cos(np.pi*lon/180.0) + ys = ns*np.cos(np.pi*lat/180.0)*np.sin(np.pi*lon/180.0) + zs = (ns*(1.0-ecc1**2.0))*np.sin(np.pi*lat/180.0) + points = np.c_[xs,ys,zs] + else: + #-- flattened model coordinates + tree = scipy.spatial.cKDTree(np.c_[gridlon[indy,indx], + gridlat[indy,indx]]) + #-- output coordinates + points = np.c_[lon,lat] + + #-- query output data points and find nearest neighbor within cutoff + dd,ii = tree.query(points,k=1,distance_upper_bound=cutoff) + ind, = np.nonzero(np.isfinite(dd)) + data.data[ind] = iflat[ii[ind]] + data.mask[ind] = False + #-- return extrapolated values + return data + +#-- PURPOSE: calculate Euclidean distances between points +def distance_matrix(c1,c2): + """ + Calculate Euclidean distances between points + + Arguments + --------- + c1: first set of coordinates + c2: second set of coordinates + + Returns + ------- + c: Euclidean distance + """ + #-- decompose Euclidean distance: (x-y)^2 = x^2 - 2xy + y^2 + dx2 = np.sum(c1**2) + dxy = np.dot(c1[np.newaxis,:], c2.T) + dy2 = np.sum(c2**2, axis=1) + #-- calculate Euclidean distance + D, = np.sqrt(dx2 - 2.0*dxy + dy2) + return D diff --git a/pyTMD/read_FES_model.py b/pyTMD/read_FES_model.py index f0eb38a2..48f48a98 100644 --- a/pyTMD/read_FES_model.py +++ b/pyTMD/read_FES_model.py @@ -1,6 +1,6 @@ #!/usr/bin/env python u""" -read_FES_model.py (09/2020) +read_FES_model.py (12/2020) Reads files for a tidal model and makes initial calculations to run tide program Includes functions to extract tidal harmonic constants from the FES (Finite Element Solution) tide models for given locations @@ -30,6 +30,7 @@ bilinear: quick bilinear interpolation spline: scipy bivariate spline interpolation linear, nearest: scipy regular grid interpolations + EXTRAPOLATE: extrapolate model using nearest-neighbors GZIP: input ascii or netCDF4 files are compressed SCALE: scaling factor for converting to output units @@ -47,9 +48,11 @@ https://unidata.github.io/netcdf4-python/netCDF4/index.html PROGRAM DEPENDENCIES: - bilinear_interp.py: bilinear interpolation of data to specified coordinates + 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 nearest-neighbor data extrapolation Updated 09/2020: set bounds error to false for regular grid interpolations adjust dimensions of input coordinates to be iterable Updated 08/2020: replaced griddata with scipy regular grid interpolators @@ -61,10 +64,12 @@ import numpy as np import scipy.interpolate from pyTMD.bilinear_interp import bilinear_interp +from pyTMD.nearest_extrap import nearest_extrap #-- PURPOSE: extract tidal harmonic constants from tide models at coordinates def extract_FES_constants(ilon, ilat, directory, model_files, - TYPE='z', VERSION=None, METHOD='spline', GZIP=True, SCALE=1): + TYPE='z', VERSION=None, METHOD='spline', EXTRAPOLATE=False, + GZIP=True, SCALE=1): """ Reads files for an ascii or netCDF4 tidal model Makes initial calculations to run the tide program @@ -93,6 +98,7 @@ def extract_FES_constants(ilon, ilat, directory, model_files, bilinear: quick bilinear interpolation spline: scipy bivariate spline interpolation linear, nearest: scipy regular grid interpolations + EXTRAPOLATE: extrapolate model using nearest-neighbors GZIP: input files are compressed SCALE: scaling factor for converting to output units @@ -166,6 +172,16 @@ def extract_FES_constants(ilon, ilat, directory, model_files, #-- replace invalid values with fill_value hci.mask[:] |= (hci.data == hci.fill_value) hci.data[hci.mask] = hci.fill_value + #-- extrapolate data using nearest-neighbors + if EXTRAPOLATE: + #-- find invalid data points + inv, = np.nonzero(hci.mask) + #-- extrapolate points within 10km of valid model points + hci.data[inv] = nearest_extrap(lon,lat,hc,ilon[inv],ilat[inv], + dtype=hc.dtype,cutoff=10.0) + #-- replace nan values with fill_value + hci.mask[inv] = np.isnan(hci.data[inv]) + hci.data[hci.mask] = hci.fill_value #-- convert amplitude from input units to meters amplitude.data[:,i] = np.abs(hci)*SCALE amplitude.mask[:,i] = np.copy(hci.mask) diff --git a/pyTMD/read_GOT_model.py b/pyTMD/read_GOT_model.py index ad05183c..49d32faa 100644 --- a/pyTMD/read_GOT_model.py +++ b/pyTMD/read_GOT_model.py @@ -1,6 +1,6 @@ #!/usr/bin/env python u""" -read_GOT_model.py (09/2020) +read_GOT_model.py (12/2020) Reads files for Richard Ray's Global Ocean Tide (GOT) models and makes initial calculations to run the tide program Includes functions to extract tidal harmonic constants out of a tidal model for @@ -17,6 +17,7 @@ bilinear: quick bilinear interpolation spline: scipy bivariate spline interpolation linear, nearest: scipy regular grid interpolations + EXTRAPOLATE: extrapolate model using nearest-neighbors GZIP: input files are compressed SCALE: scaling factor for converting to output units @@ -32,9 +33,11 @@ https://docs.scipy.org/doc/ PROGRAM DEPENDENCIES: - bilinear_interp.py: bilinear interpolation of data to specified coordinates + 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 09/2020: set bounds error to false for regular grid interpolations adjust dimensions of input coordinates to be iterable Updated 08/2020: replaced griddata with scipy regular grid interpolators @@ -56,10 +59,11 @@ import numpy as np import scipy.interpolate from pyTMD.bilinear_interp import bilinear_interp +from pyTMD.nearest_extrap import nearest_extrap #-- PURPOSE: extract tidal harmonic constants out of GOT model at coordinates def extract_GOT_constants(ilon, ilat, directory, model_files, - METHOD=None, GZIP=True, SCALE=1): + METHOD=None, EXTRAPOLATE=False, GZIP=True, SCALE=1): """ Reads files for Richard Ray's Global Ocean Tide (GOT) models Makes initial calculations to run the tide program @@ -78,6 +82,7 @@ def extract_GOT_constants(ilon, ilat, directory, model_files, bilinear: quick bilinear interpolation spline: scipy bivariate spline interpolation linear, nearest: scipy regular grid interpolations + EXTRAPOLATE: extrapolate model using nearest-neighbors GZIP: input files are compressed SCALE: scaling factor for converting to output units @@ -152,6 +157,16 @@ def extract_GOT_constants(ilon, ilat, directory, model_files, #-- replace invalid values with fill_value hci.mask[:] |= (hci.data == hci.fill_value) hci.data[hci.mask] = hci.fill_value + #-- extrapolate data using nearest-neighbors + if EXTRAPOLATE: + #-- find invalid data points + inv, = np.nonzero(hci.mask) + #-- extrapolate points within 10km of valid model points + hci.data[inv] = nearest_extrap(lon,lat,hc,ilon[inv],ilat[inv], + dtype=hc.dtype,cutoff=10.0) + #-- replace nan values with fill_value + hci.mask[inv] = np.isnan(hci.data[inv]) + hci.data[hci.mask] = hci.fill_value #-- convert amplitude from input units to meters amplitude.data[:,i] = np.abs(hci)*SCALE amplitude.mask[:,i] = np.copy(hci.mask) diff --git a/pyTMD/read_netcdf_model.py b/pyTMD/read_netcdf_model.py index 37ff33eb..cae38b11 100644 --- a/pyTMD/read_netcdf_model.py +++ b/pyTMD/read_netcdf_model.py @@ -1,6 +1,6 @@ #!/usr/bin/env python u""" -read_netcdf_model.py (11/2020) +read_netcdf_model.py (12/2020) Reads files for a tidal model and makes initial calculations to run tide program Includes functions to extract tidal harmonic constants from OTIS tide models for given locations @@ -29,6 +29,7 @@ bilinear: quick bilinear interpolation spline: scipy bivariate spline interpolation linear, nearest: scipy regular grid interpolations + EXTRAPOLATE: extrapolate model using nearest-neighbors GZIP: input netCDF4 files are compressed SCALE: scaling factor for converting to output units @@ -48,9 +49,12 @@ https://unidata.github.io/netcdf4-python/netCDF4/index.html PROGRAM DEPENDENCIES: - bilinear_interp.py: bilinear interpolation of data to specified coordinates + 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 + replace tostring with tobytes to fix DeprecationWarning Updated 11/2020: create function to read bathymetry and spatial coordinates Updated 09/2020: set bounds error to false for regular grid interpolations adjust dimensions of input coordinates to be iterable @@ -67,10 +71,11 @@ import numpy as np import scipy.interpolate from pyTMD.bilinear_interp import bilinear_interp +from pyTMD.nearest_extrap import nearest_extrap #-- PURPOSE: extract tidal harmonic constants from tide models at coordinates def extract_netcdf_constants(ilon, ilat, directory, grid_file, model_files, - TYPE='z', METHOD='spline', GZIP=True, SCALE=1): + TYPE='z', METHOD='spline', EXTRAPOLATE=False, GZIP=True, SCALE=1): """ Reads files for a netCDF4 tidal model Makes initial calculations to run the tide program @@ -96,6 +101,7 @@ def extract_netcdf_constants(ilon, ilat, directory, grid_file, model_files, bilinear: quick bilinear interpolation spline: scipy bivariate spline interpolation linear, nearest: scipy regular grid interpolations + EXTRAPOLATE: extrapolate model using nearest-neighbors GZIP: input netCDF4 files are compressed SCALE: scaling factor for converting to output units @@ -118,8 +124,6 @@ def extract_netcdf_constants(ilon, ilat, directory, grid_file, model_files, bathymetry = extend_matrix(bathymetry) #-- create masks bathymetry.mask = (bathymetry.data == 0) - #-- create meshes from latitude and longitude - gridlon,gridlat = np.meshgrid(lon,lat) #-- adjust dimensions of input coordinates to be iterable ilon = np.atleast_1d(ilon) @@ -211,6 +215,16 @@ def extract_netcdf_constants(ilon, ilat, directory, grid_file, model_files, #-- mask invalid values z1.mask[:] |= np.copy(D.mask) z1.data[z1.mask] = z1.fill_value + #-- extrapolate data using nearest-neighbors + if EXTRAPOLATE: + #-- find invalid data points + inv, = np.nonzero(z1.mask) + #-- extrapolate points within 10km of valid model points + z1.data[inv] = nearest_extrap(lon,lat,z,ilon[inv],ilat[inv], + dtype=z.dtype,cutoff=10.0) + #-- replace nan values with fill_value + z1.mask[inv] = np.isnan(z1.data[inv]) + z1.data[z1.mask] = z1.fill_value #-- amplitude and phase of the constituent ampl[:,i] = np.abs(z1) phase[:,i] = np.arctan2(-np.imag(z1),np.real(z1)) @@ -248,6 +262,16 @@ def extract_netcdf_constants(ilon, ilat, directory, grid_file, model_files, #-- mask invalid values tr1.mask[:] |= np.copy(D.mask) tr1.data[tr1.mask] = tr1.fill_value + #-- extrapolate data using nearest-neighbors + if EXTRAPOLATE: + #-- find invalid data points + inv, = np.nonzero(tr1.mask) + #-- extrapolate points within 10km of valid model points + tr1.data[inv] = nearest_extrap(lon,lat,tr,ilon[inv],ilat[inv], + dtype=tr.dtype,cutoff=10.0) + #-- replace nan values with fill_value + tr1.mask[inv] = np.isnan(tr1.data[inv]) + tr1.data[tr1.mask] = tr1.fill_value #-- convert units tr1 = tr1/unit_conv #-- amplitude and phase of the constituent @@ -386,7 +410,7 @@ def read_elevation_file(input_file,GZIP): else: fileID = netCDF4.Dataset(input_file,'r') #-- constituent name - con = fileID.variables['con'][:].tostring().decode('utf-8') + con = fileID.variables['con'][:].tobytes().decode('utf-8') #-- variable dimensions nx = fileID.dimensions['nx'].size ny = fileID.dimensions['ny'].size @@ -433,7 +457,7 @@ def read_transport_file(input_file,TYPE,GZIP): else: fileID = netCDF4.Dataset(input_file,'r') #-- constituent name - con = fileID.variables['con'][:].tostring().decode('utf-8') + con = fileID.variables['con'][:].tobytes().decode('utf-8') #-- variable dimensions nx = fileID.dimensions['nx'].size ny = fileID.dimensions['ny'].size diff --git a/pyTMD/read_tide_model.py b/pyTMD/read_tide_model.py index 27cf2d69..ad42997a 100644 --- a/pyTMD/read_tide_model.py +++ b/pyTMD/read_tide_model.py @@ -1,6 +1,6 @@ #!/usr/bin/env python u""" -read_tide_model.py (09/2020) +read_tide_model.py (12/2020) Reads files for a tidal model and makes initial calculations to run tide program Includes functions to extract tidal harmonic constants from OTIS tide models for given locations @@ -28,6 +28,7 @@ bilinear: quick bilinear interpolation spline: scipy bivariate spline interpolation linear, nearest: scipy regular grid interpolations + EXTRAPOLATE: extrapolate model using nearest-neighbors GRID: binary file type to read ATLAS: reading a global solution with localized solutions OTIS: combined global solution @@ -47,9 +48,11 @@ PROGRAM DEPENDENCIES: convert_ll_xy.py: converts lat/lon points to and from projected coordinates - bilinear_interp.py: bilinear interpolation of data to specified coordinates + 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 09/2020: set bounds error to false for regular grid interpolations adjust dimensions of input coordinates to be iterable use masked arrays with atlas models and grids. make 2' grid with nearest @@ -73,10 +76,11 @@ import scipy.interpolate from pyTMD.convert_ll_xy import convert_ll_xy from pyTMD.bilinear_interp import bilinear_interp +from pyTMD.nearest_extrap import nearest_extrap #-- PURPOSE: extract tidal harmonic constants from tide models at coordinates def extract_tidal_constants(ilon, ilat, grid_file, model_file, EPSG, TYPE='z', - METHOD='spline', GRID='OTIS'): + METHOD='spline', EXTRAPOLATE=False, GRID='OTIS'): """ Reads files for an OTIS-formatted tidal model Makes initial calculations to run the tide program @@ -100,6 +104,7 @@ def extract_tidal_constants(ilon, ilat, grid_file, model_file, EPSG, TYPE='z', bilinear: quick bilinear interpolation spline: scipy bivariate spline interpolation linear, nearest: scipy regular grid interpolations + EXTRAPOLATE: extrapolate model using nearest-neighbors GRID: binary file type to read ATLAS: reading a global solution with localized solutions OTIS: combined global solution @@ -257,6 +262,16 @@ def extract_tidal_constants(ilon, ilat, grid_file, model_file, EPSG, TYPE='z', #-- replace invalid values with fill_value z1.mask = (z1.data == z1.fill_value) | (~mz1.astype(np.bool)) z1.data[z1.mask] = z1.fill_value + #-- extrapolate data using nearest-neighbors + if EXTRAPOLATE: + #-- find invalid data points + inv, = np.nonzero(z1.mask) + #-- extrapolate points within 10km of valid model points + z1.data[inv] = nearest_extrap(xi,yi,z,x[inv],y[inv], + dtype=np.complex128,cutoff=10.0,EPSG=EPSG) + #-- replace nan values with fill_value + z1.mask[inv] = np.isnan(z1.data[inv]) + z1.data[z1.mask] = z1.fill_value #-- amplitude and phase of the constituent amplitude.data[:,i] = np.abs(z1.data) amplitude.mask[:,i] = np.copy(z1.mask) @@ -302,6 +317,16 @@ def extract_tidal_constants(ilon, ilat, grid_file, model_file, EPSG, TYPE='z', #-- replace invalid values with fill_value u1.mask = (u1.data == u1.fill_value) | (~mu1.astype(np.bool)) u1.data[u1.mask] = u1.fill_value + #-- extrapolate data using nearest-neighbors + if EXTRAPOLATE: + #-- find invalid data points + inv, = np.nonzero(u1.mask) + #-- extrapolate points within 10km of valid model points + u1.data[inv] = nearest_extrap(xu,yi,u,x[inv],y[inv], + dtype=np.complex128,cutoff=10.0,EPSG=EPSG) + #-- replace nan values with fill_value + u1.mask[inv] = np.isnan(u1.data[inv]) + u1.data[u1.mask] = u1.fill_value #-- convert units u1 = u1/unit_conv #-- amplitude and phase of the constituent @@ -349,6 +374,16 @@ def extract_tidal_constants(ilon, ilat, grid_file, model_file, EPSG, TYPE='z', #-- replace invalid values with fill_value v1.mask = (v1.data == v1.fill_value) | (~mv1.astype(np.bool)) v1.data[v1.mask] = v1.fill_value + #-- extrapolate data using nearest-neighbors + if EXTRAPOLATE: + #-- find invalid data points + inv, = np.nonzero(v1.mask) + #-- extrapolate points within 10km of valid model points + v1.data[inv] = nearest_extrap(x,yv,v,x[inv],y[inv], + dtype=np.complex128,cutoff=10.0,EPSG=EPSG) + #-- replace nan values with fill_value + v1.mask[inv] = np.isnan(v1.data[inv]) + v1.data[v1.mask] = v1.fill_value #-- convert units v1 = v1/unit_conv #-- amplitude and phase of the constituent diff --git a/pyTMD/utilities.py b/pyTMD/utilities.py index 3fe40ba9..a585f0c9 100644 --- a/pyTMD/utilities.py +++ b/pyTMD/utilities.py @@ -63,14 +63,16 @@ def get_data_path(relpath): #-- PURPOSE: get the MD5 hash value of a file def 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 """ - #-- check if local file exists - if os.access(os.path.expanduser(local),os.F_OK): + #-- check if open file object or if local file exists + if isinstance(local, io.IOBase): + return hashlib.md5(local.getvalue()).hexdigest() + elif os.access(os.path.expanduser(local),os.F_OK): #-- generate checksum hash for local file #-- open the local_file in binary read mode with open(os.path.expanduser(local), 'rb') as local_buffer: @@ -105,10 +107,12 @@ def roman_to_int(roman): """ #-- mapping between Roman and Arabic numerals roman_map = {'i':1, 'v':5, 'x':10, 'l':50, 'c':100, 'd':500, 'm':1000} + #-- verify case + roman = roman.lower() output = 0 #-- iterate through roman numerals in string and calculate total for i,s in enumerate(roman): - if i > 0 and roman_map[s] > roman_map[roman[i-1]]: + if (i > 0) and (roman_map[s] > roman_map[roman[i-1]]): output += roman_map[s] - 2*roman_map[roman[i-1]] else: output += roman_map[s] diff --git a/scripts/compute_LPET_ICESat2_ATL11.py b/scripts/compute_LPET_ICESat2_ATL11.py index 088f1b9f..ecf2e3a8 100644 --- a/scripts/compute_LPET_ICESat2_ATL11.py +++ b/scripts/compute_LPET_ICESat2_ATL11.py @@ -100,10 +100,10 @@ def compute_LPET_ICESat2(FILE, VERBOSE=False, MODE=0o775): #-- for each input beam pair within the file for ptx in sorted(IS2_atl11_pairs): #-- output data dictionaries for beam - IS2_atl11_tide[ptx] = {} - IS2_atl11_fill[ptx] = {} - IS2_atl11_dims[ptx] = {} - IS2_atl11_tide_attrs[ptx] = {} + IS2_atl11_tide[ptx] = dict(cycle_stats={}) + IS2_atl11_fill[ptx] = dict(cycle_stats={}) + IS2_atl11_dims[ptx] = dict(cycle_stats={}) + IS2_atl11_tide_attrs[ptx] = dict(cycle_stats={}) #-- number of average segments and number of included cycles delta_time = fileID[ptx]['delta_time'][:].copy() @@ -227,10 +227,6 @@ def compute_LPET_ICESat2(FILE, VERBOSE=False, MODE=0o775): "delta_time latitude longitude" #-- cycle statistics variables - IS2_atl11_tide[ptx]['cycle_stats'] = {} - IS2_atl11_fill[ptx]['cycle_stats'] = {} - IS2_atl11_dims[ptx]['cycle_stats'] = {} - IS2_atl11_tide_attrs[ptx]['cycle_stats'] = {} IS2_atl11_tide_attrs[ptx]['cycle_stats']['Description'] = ("The cycle_stats subgroup " "contains summary information about segments for each reference point, including " "the uncorrected mean heights for reference surfaces, blowing snow and cloud " @@ -326,7 +322,7 @@ def HDF5_ATL11_tide_write(IS2_atl11_tide, IS2_atl11_attrs, INPUT=None, h5[ptx][k].dims[i].attach_scale(h5[ptx][dim]) else: #-- make dimension - h5[ptx][k].make_scale(k) + h5[ptx][k].make_scale(k) #-- add HDF5 variable attributes for att_name,att_val in attrs.items(): h5[ptx][k].attrs[att_name] = att_val diff --git a/scripts/compute_LPET_ICESat_GLA12.py b/scripts/compute_LPET_ICESat_GLA12.py index 6c7d547d..fe74225a 100644 --- a/scripts/compute_LPET_ICESat_GLA12.py +++ b/scripts/compute_LPET_ICESat_GLA12.py @@ -66,7 +66,7 @@ def compute_LPET_ICESat(FILE, VERBOSE=False, MODE=0o775): #-- new reference orbit ground track file is obtained.) #-- INST: Instance number (increments every time the satellite enters a #-- different reference orbit) - #-- CYCL: Cycle of reference orbit for this phase + #-- CYCL: Cycle of reference orbit for this phase #-- TRK: Track within reference orbit #-- SEG: Segment of orbit #-- GRAN: Granule version number @@ -264,7 +264,7 @@ def HDF5_GLA12_tide_write(IS_gla12_tide, IS_gla12_attrs, #-- add HDF5 group attributes for att_name,att_val in IS_gla12_attrs['Data_40HZ'][group].items(): if not isinstance(att_val,dict): - fileID['Data_40HZ'][group].attrs[att_name] = att_val + fileID['Data_40HZ'][group].attrs[att_name] = att_val #-- for each variable in the group for key,val in IS_gla12_tide['Data_40HZ'][group].items(): fillvalue = FILL_VALUE['Data_40HZ'][group][key] diff --git a/scripts/compute_LPT_ICESat_GLA12.py b/scripts/compute_LPT_ICESat_GLA12.py index 022f323c..445c8e65 100644 --- a/scripts/compute_LPT_ICESat_GLA12.py +++ b/scripts/compute_LPT_ICESat_GLA12.py @@ -70,7 +70,7 @@ def compute_LPT_ICESat(FILE, VERBOSE=False, MODE=0o775): #-- new reference orbit ground track file is obtained.) #-- INST: Instance number (increments every time the satellite enters a #-- different reference orbit) - #-- CYCL: Cycle of reference orbit for this phase + #-- CYCL: Cycle of reference orbit for this phase #-- TRK: Track within reference orbit #-- SEG: Segment of orbit #-- GRAN: Granule version number @@ -341,7 +341,7 @@ def HDF5_GLA12_tide_write(IS_gla12_tide, IS_gla12_attrs, #-- add HDF5 group attributes for att_name,att_val in IS_gla12_attrs['Data_40HZ'][group].items(): if not isinstance(att_val,dict): - fileID['Data_40HZ'][group].attrs[att_name] = att_val + fileID['Data_40HZ'][group].attrs[att_name] = att_val #-- for each variable in the group for key,val in IS_gla12_tide['Data_40HZ'][group].items(): fillvalue = FILL_VALUE['Data_40HZ'][group][key] diff --git a/scripts/compute_OPT_ICESat_GLA12.py b/scripts/compute_OPT_ICESat_GLA12.py index 4bcb5f1f..fa5c95f9 100644 --- a/scripts/compute_OPT_ICESat_GLA12.py +++ b/scripts/compute_OPT_ICESat_GLA12.py @@ -76,7 +76,7 @@ def compute_OPT_ICESat(FILE, METHOD=None, VERBOSE=False, MODE=0o775): #-- new reference orbit ground track file is obtained.) #-- INST: Instance number (increments every time the satellite enters a #-- different reference orbit) - #-- CYCL: Cycle of reference orbit for this phase + #-- CYCL: Cycle of reference orbit for this phase #-- TRK: Track within reference orbit #-- SEG: Segment of orbit #-- GRAN: Granule version number @@ -354,7 +354,7 @@ def HDF5_GLA12_tide_write(IS_gla12_tide, IS_gla12_attrs, #-- add HDF5 group attributes for att_name,att_val in IS_gla12_attrs['Data_40HZ'][group].items(): if not isinstance(att_val,dict): - fileID['Data_40HZ'][group].attrs[att_name] = att_val + fileID['Data_40HZ'][group].attrs[att_name] = att_val #-- for each variable in the group for key,val in IS_gla12_tide['Data_40HZ'][group].items(): fillvalue = FILL_VALUE['Data_40HZ'][group][key] diff --git a/scripts/compute_tidal_currents.py b/scripts/compute_tidal_currents.py index 47e0e8ed..d9072571 100755 --- a/scripts/compute_tidal_currents.py +++ b/scripts/compute_tidal_currents.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_tidal_currents.py -Written by Tyler Sutterley (11/2020) +Written by Tyler Sutterley (12/2020) Calculates zonal and meridional tidal currents for an input file Uses OTIS format tidal solutions provided by Ohio State University and ESR @@ -54,6 +54,7 @@ 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 @@ -87,9 +88,12 @@ read_tide_model.py: extract tidal harmonic constants from OTIS tide models read_netcdf_model.py: extract tidal harmonic constants from netcdf 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 predict_tide_drift.py: predict tidal elevations using harmonic constants UPDATE HISTORY: + Updated 12/2020: added valid data extrapolation with nearest_extrap Updated 11/2020: added options to read from and write to geotiff image files Updated 10/2020: using argparse to set command line parameters Forked 09/2020 from compute_tidal_elevations.py @@ -127,7 +131,8 @@ def compute_tidal_currents(tide_dir, input_file, output_file, TIDE_MODEL=None, FORMAT='csv', VARIABLES=['time','lat','lon','data'], HEADER=0, TYPE='drift', TIME_UNITS='days since 1858-11-17T00:00:00', - TIME=None, PROJECTION='4326', METHOD='spline', VERBOSE=False, MODE=0o775): + TIME=None, PROJECTION='4326', METHOD='spline', EXTRAPOLATE=False, + VERBOSE=False, MODE=0o775): #-- select between tide models if (TIDE_MODEL == 'CATS0201'): @@ -344,17 +349,18 @@ def compute_tidal_currents(tide_dir, input_file, output_file, if model_format in ('OTIS','ATLAS'): amp,ph,D,c = extract_tidal_constants(lon.flatten(), lat.flatten(), grid_file, model_file, EPSG, TYPE=t, METHOD=METHOD, - GRID=model_format) + EXTRAPOLATE=EXTRAPOLATE, GRID=model_format) deltat = np.zeros((nt)) elif (model_format == 'netcdf'): amp,ph,D,c = extract_netcdf_constants(lon.flatten(), lat.flatten(), model_directory, grid_file, model_files[t], TYPE=t, - METHOD=METHOD, SCALE=model_scale, GZIP=GZIP) + METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE, SCALE=model_scale, + GZIP=GZIP) deltat = np.zeros((nt)) elif (model_format == 'FES'): amp,ph = extract_FES_constants(lon.flatten(), lat.flatten(), model_directory[t], model_files, TYPE=t, VERSION=TIDE_MODEL, - METHOD=METHOD, SCALE=model_scale) + METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE, SCALE=model_scale) #-- interpolate delta times from calendar dates to tide time delta_file = get_data_path(['data','merged_deltat.data']) deltat = calc_delta_time(delta_file, tide_time) @@ -471,6 +477,10 @@ def main(): metavar='METHOD', type=str, default='spline', choices=('spline','linear','nearest','bilinear'), help='Spatial interpolation method') + #-- extrapolate with nearest-neighbors + parser.add_argument('--extrapolate','-E', + default=False, action='store_true', + help='Extrapolate with nearest-neighbors') #-- verbose output of processing run #-- print information about each input and output file parser.add_argument('--verbose','-V', @@ -493,7 +503,8 @@ def main(): FORMAT=args.format, TIDE_MODEL=args.tide, VARIABLES=args.variables, HEADER=args.header, TYPE=args.type, TIME_UNITS=args.epoch, TIME=args.deltatime, PROJECTION=args.projection, - METHOD=args.interpolate, VERBOSE=args.verbose, MODE=args.mode) + METHOD=args.interpolate, EXTRAPOLATE=args.extrapolate, + VERBOSE=args.verbose, MODE=args.mode) #-- run main program if __name__ == '__main__': diff --git a/scripts/compute_tidal_elevations.py b/scripts/compute_tidal_elevations.py index 7dc3742e..d0e591bf 100755 --- a/scripts/compute_tidal_elevations.py +++ b/scripts/compute_tidal_elevations.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_tidal_elevations.py -Written by Tyler Sutterley (11/2020) +Written by Tyler Sutterley (12/2020) Calculates tidal elevations for an input file Uses OTIS format tidal solutions provided by Ohio State University and ESR @@ -64,6 +64,7 @@ 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 @@ -98,9 +99,12 @@ 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 predict_tide_drift.py: predict tidal elevations using harmonic constants UPDATE HISTORY: + Updated 12/2020: added valid data extrapolation with nearest_extrap Updated 11/2020: added model constituents from TPXO9-atlas-v3 added options to read from and write to geotiff image files Updated 10/2020: using argparse to set command line parameters @@ -139,7 +143,8 @@ def compute_tidal_elevations(tide_dir, input_file, output_file, TIDE_MODEL=None, FORMAT='csv', VARIABLES=['time','lat','lon','data'], HEADER=0, TYPE='drift', TIME_UNITS='days since 1858-11-17T00:00:00', - TIME=None, PROJECTION='4326', METHOD='spline', VERBOSE=False, MODE=0o775): + TIME=None, PROJECTION='4326', METHOD='spline', EXTRAPOLATE=False, + VERBOSE=False, MODE=0o775): #-- select between tide models if (TIDE_MODEL == 'CATS0201'): @@ -470,16 +475,18 @@ def compute_tidal_elevations(tide_dir, input_file, output_file, #-- read tidal constants and interpolate to grid points if model_format in ('OTIS','ATLAS'): amp,ph,D,c = extract_tidal_constants(lon.flatten(), lat.flatten(), - grid_file, model_file, EPSG, TYPE=model_type, METHOD=METHOD) + grid_file, model_file, EPSG, TYPE=model_type, METHOD=METHOD, + EXTRAPOLATE=EXTRAPOLATE) deltat = np.zeros((nt)) elif (model_format == 'netcdf'): amp,ph,D,c = extract_netcdf_constants(lon.flatten(), lat.flatten(), model_directory, grid_file, model_files, TYPE=model_type, - METHOD=METHOD, SCALE=model_scale) + METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE, SCALE=model_scale) deltat = np.zeros((nt)) elif (model_format == 'GOT'): amp,ph = extract_GOT_constants(lon.flatten(), lat.flatten(), - model_directory, model_files, METHOD=METHOD, SCALE=model_scale) + model_directory, model_files, METHOD=METHOD, + EXTRAPOLATE=EXTRAPOLATE, SCALE=model_scale) #-- convert times from modified julian days to days since 1992-01-01 #-- interpolate delta times from calendar dates to tide time delta_file = get_data_path(['data','merged_deltat.data']) @@ -487,7 +494,7 @@ def compute_tidal_elevations(tide_dir, input_file, output_file, elif (model_format == 'FES'): amp,ph = extract_FES_constants(lon.flatten(), lat.flatten(), model_directory, model_files, TYPE=model_type, VERSION=TIDE_MODEL, - METHOD=METHOD, SCALE=model_scale) + METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE, SCALE=model_scale) #-- convert times from modified julian days to days since 1992-01-01 #-- interpolate delta times from calendar dates to tide time delta_file = get_data_path(['data','merged_deltat.data']) @@ -604,6 +611,10 @@ def main(): metavar='METHOD', type=str, default='spline', choices=('spline','linear','nearest','bilinear'), help='Spatial interpolation method') + #-- extrapolate with nearest-neighbors + parser.add_argument('--extrapolate','-E', + default=False, action='store_true', + help='Extrapolate with nearest-neighbors') #-- verbose output of processing run #-- print information about each input and output file parser.add_argument('--verbose','-V', @@ -626,7 +637,8 @@ def main(): FORMAT=args.format, TIDE_MODEL=args.tide, VARIABLES=args.variables, HEADER=args.header, TYPE=args.type, TIME_UNITS=args.epoch, TIME=args.deltatime, PROJECTION=args.projection, - METHOD=args.interpolate, VERBOSE=args.verbose, MODE=args.mode) + METHOD=args.interpolate, EXTRAPOLATE=args.extrapolate, + VERBOSE=args.verbose, MODE=args.mode) #-- run main program if __name__ == '__main__': diff --git a/scripts/compute_tides_ICESat2_ATL03.py b/scripts/compute_tides_ICESat2_ATL03.py index 117d4408..05beb305 100644 --- a/scripts/compute_tides_ICESat2_ATL03.py +++ b/scripts/compute_tides_ICESat2_ATL03.py @@ -42,6 +42,7 @@ linear nearest bilinear + -E X, --extrapolate X: Extrapolate with nearest-neighbors -M X, --mode X: Permission mode of directories and files created -V, --verbose: Output information about each created file @@ -67,14 +68,17 @@ infer_minor_corrections.py: return corrections for minor constituents load_constituent.py: loads parameters for a given tidal constituent load_nodal_corrections.py: load the nodal corrections for tidal constituents - predict_tide_drift.py: predict tidal elevations using harmonic constants read_tide_model.py: extract tidal harmonic constants from OTIS tide models 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 + predict_tide_drift.py: predict tidal elevations using harmonic constants UPDATE HISTORY: Updated 12/2020: H5py deprecation warning change to use make_scale + added valid data extrapolation with nearest_extrap Updated 11/2020: added model constituents from TPXO9-atlas-v3 Updated 10/2020: using argparse to set command line parameters Updated 08/2020: using builtin time operations. python3 regular expressions @@ -116,7 +120,7 @@ #-- PURPOSE: read ICESat-2 geolocated photon data (ATL03) from NSIDC #-- compute tides at points and times using tidal model driver algorithms def compute_tides_ICESat2(tide_dir, FILE, TIDE_MODEL=None, METHOD='spline', - VERBOSE=False, MODE=0o775): + EXTRAPOLATE=False, VERBOSE=False, MODE=0o775): #-- select between tide models if (TIDE_MODEL == 'CATS0201'): grid_file = os.path.join(tide_dir,'cats0201_tmd','grid_CATS') @@ -478,21 +482,24 @@ def compute_tides_ICESat2(tide_dir, FILE, TIDE_MODEL=None, METHOD='spline', #-- 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, EPSG, TYPE=TYPE, METHOD=METHOD, GRID=model_format) + model_file, EPSG, TYPE=TYPE, METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE, + GRID=model_format) deltat = np.zeros_like(tide_time) elif (model_format == 'netcdf'): amp,ph,D,c = extract_netcdf_constants(lon, lat, model_directory, - grid_file, model_files, TYPE=TYPE, METHOD=METHOD, SCALE=SCALE) + grid_file, model_files, TYPE=TYPE, METHOD=METHOD, + EXTRAPOLATE=EXTRAPOLATE, SCALE=SCALE) deltat = np.zeros_like(tide_time) elif (model_format == 'GOT'): amp,ph = extract_GOT_constants(lon, lat, model_directory, - model_files, METHOD=METHOD, SCALE=SCALE) + model_files, 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, tide_time) elif (model_format == 'FES'): amp,ph = extract_FES_constants(lon,lat,model_directory,model_files, - TYPE=TYPE, VERSION=TIDE_MODEL, METHOD=METHOD, SCALE=SCALE) + TYPE=TYPE, VERSION=TIDE_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, tide_time) @@ -821,6 +828,10 @@ def main(): metavar='METHOD', type=str, default='spline', choices=('spline','linear','nearest','bilinear'), help='Spatial interpolation method') + #-- extrapolate with nearest-neighbors + parser.add_argument('--extrapolate','-E', + default=False, action='store_true', + help='Extrapolate with nearest-neighbors') #-- verbosity settings #-- verbose will output information about each output file parser.add_argument('--verbose','-V', @@ -835,7 +846,8 @@ def main(): #-- run for each input ATL03 file for FILE in args.infile: compute_tides_ICESat2(args.directory, FILE, TIDE_MODEL=args.tide, - METHOD=args.interpolate, VERBOSE=args.verbose, MODE=args.mode) + METHOD=args.interpolate, EXTRAPOLATE=args.extrapolate, + VERBOSE=args.verbose, MODE=args.mode) #-- run main program if __name__ == '__main__': diff --git a/scripts/compute_tides_ICESat2_ATL06.py b/scripts/compute_tides_ICESat2_ATL06.py index aaff9daa..26cd5932 100644 --- a/scripts/compute_tides_ICESat2_ATL06.py +++ b/scripts/compute_tides_ICESat2_ATL06.py @@ -40,6 +40,7 @@ linear nearest bilinear + -E X, --extrapolate X: Extrapolate with nearest-neighbors -M X, --mode X: Permission mode of directories and files created -V, --verbose: Output information about each created file @@ -65,14 +66,17 @@ infer_minor_corrections.py: return corrections for minor constituents load_constituent.py: loads parameters for a given tidal constituent load_nodal_corrections.py: load the nodal corrections for tidal constituents - predict_tide_drift.py: predict tidal elevations using harmonic constants read_tide_model.py: extract tidal harmonic constants from OTIS tide models 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 + predict_tide_drift.py: predict tidal elevations using harmonic constants UPDATE HISTORY: Updated 12/2020: H5py deprecation warning change to use make_scale + added valid data extrapolation with nearest_extrap Updated 11/2020: added model constituents from TPXO9-atlas-v3 Updated 10/2020: using argparse to set command line parameters Updated 08/2020: using builtin time operations. python3 regular expressions @@ -116,7 +120,7 @@ #-- PURPOSE: read ICESat-2 land ice data (ATL06) from NSIDC #-- compute tides at points and times using tidal model driver algorithms def compute_tides_ICESat2(tide_dir, FILE, TIDE_MODEL=None, METHOD='spline', - VERBOSE=False, MODE=0o775): + EXTRAPOLATE=False, VERBOSE=False, MODE=0o775): #-- select between tide models if (TIDE_MODEL == 'CATS0201'): grid_file = os.path.join(tide_dir,'cats0201_tmd','grid_CATS') @@ -473,23 +477,25 @@ def compute_tides_ICESat2(tide_dir, FILE, TIDE_MODEL=None, METHOD='spline', if model_format in ('OTIS','ATLAS'): amp,ph,D,c = extract_tidal_constants(val['longitude'], val['latitude'], grid_file, model_file, EPSG, TYPE=TYPE, - METHOD=METHOD, GRID=model_format) + METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE, GRID=model_format) deltat = np.zeros_like(tide_time) elif (model_format == 'netcdf'): amp,ph,D,c = extract_netcdf_constants(val['longitude'], val['latitude'], model_directory, grid_file, - model_files, TYPE=TYPE, METHOD=METHOD, SCALE=SCALE) + model_files, TYPE=TYPE, METHOD=METHOD, + EXTRAPOLATE=EXTRAPOLATE, SCALE=SCALE) deltat = np.zeros_like(tide_time) elif (model_format == 'GOT'): amp,ph = extract_GOT_constants(val['longitude'], val['latitude'], - model_directory, model_files, METHOD=METHOD, SCALE=SCALE) + model_directory, model_files, 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, tide_time) elif (model_format == 'FES'): amp,ph = extract_FES_constants(val['longitude'], val['latitude'], model_directory, model_files, TYPE=TYPE, VERSION=TIDE_MODEL, - 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, tide_time) @@ -834,6 +840,10 @@ def main(): metavar='METHOD', type=str, default='spline', choices=('spline','linear','nearest','bilinear'), help='Spatial interpolation method') + #-- extrapolate with nearest-neighbors + parser.add_argument('--extrapolate','-E', + default=False, action='store_true', + help='Extrapolate with nearest-neighbors') #-- verbosity settings #-- verbose will output information about each output file parser.add_argument('--verbose','-V', @@ -848,7 +858,8 @@ def main(): #-- run for each input ATL06 file for FILE in args.infile: compute_tides_ICESat2(args.directory, FILE, TIDE_MODEL=args.tide, - METHOD=args.interpolate, VERBOSE=args.verbose, MODE=args.mode) + METHOD=args.interpolate, EXTRAPOLATE=args.extrapolate, + VERBOSE=args.verbose, MODE=args.mode) #-- run main program if __name__ == '__main__': diff --git a/scripts/compute_tides_ICESat2_ATL07.py b/scripts/compute_tides_ICESat2_ATL07.py index ac27c08f..82e196f6 100644 --- a/scripts/compute_tides_ICESat2_ATL07.py +++ b/scripts/compute_tides_ICESat2_ATL07.py @@ -40,6 +40,7 @@ linear nearest bilinear + -E X, --extrapolate X: Extrapolate with nearest-neighbors -M X, --mode X: Permission mode of directories and files created -V, --verbose: Output information about each created file @@ -65,14 +66,17 @@ infer_minor_corrections.py: return corrections for minor constituents load_constituent.py: loads parameters for a given tidal constituent load_nodal_corrections.py: load the nodal corrections for tidal constituents - predict_tide_drift.py: predict tidal elevations using harmonic constants read_tide_model.py: extract tidal harmonic constants from OTIS tide models 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 + predict_tide_drift.py: predict tidal elevations using harmonic constants UPDATE HISTORY: Updated 12/2020: H5py deprecation warning change to use make_scale + added valid data extrapolation with nearest_extrap Updated 11/2020: added model constituents from TPXO9-atlas-v3 Updated 10/2020: using argparse to set command line parameters Updated 08/2020: using builtin time operations. python3 regular expressions @@ -115,7 +119,7 @@ #-- PURPOSE: read ICESat-2 sea ice height (ATL07) from NSIDC #-- compute tides at points and times using tidal model driver algorithms def compute_tides_ICESat2(tide_dir, FILE, TIDE_MODEL=None, METHOD='spline', - VERBOSE=False, MODE=0o775): + EXTRAPOLATE=False, VERBOSE=False, MODE=0o775): #-- select between tide models if (TIDE_MODEL == 'CATS0201'): grid_file = os.path.join(tide_dir,'cats0201_tmd','grid_CATS') @@ -470,23 +474,25 @@ def compute_tides_ICESat2(tide_dir, FILE, TIDE_MODEL=None, METHOD='spline', if model_format in ('OTIS','ATLAS'): amp,ph,D,c = extract_tidal_constants(val['longitude'], val['latitude'], grid_file, model_file, EPSG, TYPE=TYPE, - METHOD=METHOD, GRID=model_format) + METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE, GRID=model_format) deltat = np.zeros_like(tide_time) elif (model_format == 'netcdf'): amp,ph,D,c = extract_netcdf_constants(val['longitude'], val['latitude'], model_directory, grid_file, - model_files, TYPE=TYPE, METHOD=METHOD, SCALE=SCALE) + model_files, TYPE=TYPE, METHOD=METHOD, + EXTRAPOLATE=EXTRAPOLATE, SCALE=SCALE) deltat = np.zeros_like(tide_time) elif (model_format == 'GOT'): amp,ph = extract_GOT_constants(val['longitude'], val['latitude'], - model_directory, model_files, METHOD=METHOD, SCALE=SCALE) + model_directory, model_files, 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, tide_time) elif (model_format == 'FES'): amp,ph = extract_FES_constants(val['longitude'], val['latitude'], model_directory, model_files, TYPE=TYPE, VERSION=TIDE_MODEL, - 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, tide_time) @@ -861,6 +867,10 @@ def main(): metavar='METHOD', type=str, default='spline', choices=('spline','linear','nearest','bilinear'), help='Spatial interpolation method') + #-- extrapolate with nearest-neighbors + parser.add_argument('--extrapolate','-E', + default=False, action='store_true', + help='Extrapolate with nearest-neighbors') #-- verbosity settings #-- verbose will output information about each output file parser.add_argument('--verbose','-V', @@ -875,7 +885,8 @@ def main(): #-- run for each input ATL07 file for FILE in args.infile: compute_tides_ICESat2(args.directory, FILE, TIDE_MODEL=args.tide, - METHOD=args.interpolate, VERBOSE=args.verbose, MODE=args.mode) + METHOD=args.interpolate, EXTRAPOLATE=args.extrapolate, + VERBOSE=args.verbose, MODE=args.mode) #-- run main program if __name__ == '__main__': diff --git a/scripts/compute_tides_ICESat2_ATL11.py b/scripts/compute_tides_ICESat2_ATL11.py index 806feab2..1ca207d0 100644 --- a/scripts/compute_tides_ICESat2_ATL11.py +++ b/scripts/compute_tides_ICESat2_ATL11.py @@ -40,6 +40,7 @@ linear nearest bilinear + -E X, --extrapolate X: Extrapolate with nearest-neighbors -M X, --mode X: Permission mode of directories and files created -V, --verbose: Output information about each created file @@ -64,11 +65,13 @@ infer_minor_corrections.py: return corrections for minor constituents load_constituent.py: loads parameters for a given tidal constituent load_nodal_corrections.py: load the nodal corrections for tidal constituents - predict_tide_drift.py: predict tidal elevations using harmonic constants read_tide_model.py: extract tidal harmonic constants from OTIS tide models 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 + predict_tide_drift.py: predict tidal elevations using harmonic constants UPDATE HISTORY: Written 12/2020 @@ -96,7 +99,7 @@ #-- PURPOSE: read ICESat-2 annual land ice height data (ATL11) from NSIDC #-- compute tides at points and times using tidal model driver algorithms def compute_tides_ICESat2(tide_dir, FILE, TIDE_MODEL=None, METHOD='spline', - VERBOSE=False, MODE=0o775): + EXTRAPOLATE=False, VERBOSE=False, MODE=0o775): #-- select between tide models if (TIDE_MODEL == 'CATS0201'): grid_file = os.path.join(tide_dir,'cats0201_tmd','grid_CATS') @@ -446,10 +449,10 @@ def compute_tides_ICESat2(tide_dir, FILE, TIDE_MODEL=None, METHOD='spline', #-- for each input beam pair within the file for ptx in sorted(IS2_atl11_pairs): #-- output data dictionaries for beam - IS2_atl11_tide[ptx] = {} - IS2_atl11_fill[ptx] = {} - IS2_atl11_dims[ptx] = {} - IS2_atl11_tide_attrs[ptx] = {} + IS2_atl11_tide[ptx] = dict(cycle_stats={}) + IS2_atl11_fill[ptx] = dict(cycle_stats={}) + IS2_atl11_dims[ptx] = dict(cycle_stats={}) + IS2_atl11_tide_attrs[ptx] = dict(cycle_stats={}) #-- number of average segments and number of included cycles delta_time = fileID[ptx]['delta_time'][:].copy() @@ -467,24 +470,27 @@ def compute_tides_ICESat2(tide_dir, FILE, TIDE_MODEL=None, METHOD='spline', if model_format in ('OTIS','ATLAS'): amp,ph,D,c = extract_tidal_constants(fileID[ptx]['longitude'][:], fileID[ptx]['latitude'][:], grid_file, model_file, EPSG, - TYPE=TYPE, METHOD=METHOD, GRID=model_format) + TYPE=TYPE, METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE, + GRID=model_format) deltat = np.zeros_like(tide_time) elif (model_format == 'netcdf'): amp,ph,D,c = extract_netcdf_constants(fileID[ptx]['longitude'][:], fileID[ptx]['latitude'][:], model_directory, grid_file, - model_files, TYPE=TYPE, METHOD=METHOD, SCALE=SCALE) + model_files, TYPE=TYPE, METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE, + SCALE=SCALE) deltat = np.zeros_like(tide_time) elif (model_format == 'GOT'): amp,ph = extract_GOT_constants(fileID[ptx]['longitude'][:], fileID[ptx]['latitude'][:], 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, tide_time) elif (model_format == 'FES'): amp,ph = extract_FES_constants(fileID[ptx]['longitude'][:], fileID[ptx]['latitude'][:], model_directory, model_files, - TYPE=TYPE, VERSION=TIDE_MODEL, METHOD=METHOD, SCALE=SCALE) + TYPE=TYPE, VERSION=TIDE_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, tide_time) @@ -603,10 +609,6 @@ def compute_tides_ICESat2(tide_dir, FILE, TIDE_MODEL=None, METHOD='spline', "delta_time latitude longitude" #-- cycle statistics variables - IS2_atl11_tide[ptx]['cycle_stats'] = {} - IS2_atl11_fill[ptx]['cycle_stats'] = {} - IS2_atl11_dims[ptx]['cycle_stats'] = {} - IS2_atl11_tide_attrs[ptx]['cycle_stats'] = {} IS2_atl11_tide_attrs[ptx]['cycle_stats']['Description'] = ("The cycle_stats subgroup " "contains summary information about segments for each reference point, including " "the uncorrected mean heights for reference surfaces, blowing snow and cloud " @@ -700,7 +702,7 @@ def HDF5_ATL11_tide_write(IS2_atl11_tide, IS2_atl11_attrs, INPUT=None, h5[ptx][k].dims[i].attach_scale(h5[ptx][dim]) else: #-- make dimension - h5[ptx][k].make_scale(k) + h5[ptx][k].make_scale(k) #-- add HDF5 variable attributes for att_name,att_val in attrs.items(): h5[ptx][k].attrs[att_name] = att_val @@ -833,6 +835,10 @@ def main(): metavar='METHOD', type=str, default='spline', choices=('spline','linear','nearest','bilinear'), help='Spatial interpolation method') + #-- extrapolate with nearest-neighbors + parser.add_argument('--extrapolate','-E', + default=False, action='store_true', + help='Extrapolate with nearest-neighbors') #-- verbosity settings #-- verbose will output information about each output file parser.add_argument('--verbose','-V', @@ -847,7 +853,8 @@ def main(): #-- run for each input ATL11 file for FILE in args.infile: compute_tides_ICESat2(args.directory, FILE, TIDE_MODEL=args.tide, - METHOD=args.interpolate, VERBOSE=args.verbose, MODE=args.mode) + METHOD=args.interpolate, EXTRAPOLATE=args.extrapolate, + VERBOSE=args.verbose, MODE=args.mode) #-- run main program if __name__ == '__main__': diff --git a/scripts/compute_tides_ICESat2_ATL12.py b/scripts/compute_tides_ICESat2_ATL12.py index 388776c9..8e7c0050 100644 --- a/scripts/compute_tides_ICESat2_ATL12.py +++ b/scripts/compute_tides_ICESat2_ATL12.py @@ -40,6 +40,7 @@ linear nearest bilinear + -E X, --extrapolate X: Extrapolate with nearest-neighbors -M X, --mode X: Permission mode of directories and files created -V, --verbose: Output information about each created file @@ -65,14 +66,17 @@ infer_minor_corrections.py: return corrections for minor constituents load_constituent.py: loads parameters for a given tidal constituent load_nodal_corrections.py: load the nodal corrections for tidal constituents - predict_tide_drift.py: predict tidal elevations using harmonic constants read_tide_model.py: extract tidal harmonic constants from OTIS tide models 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 + predict_tide_drift.py: predict tidal elevations using harmonic constants UPDATE HISTORY: Updated 12/2020: H5py deprecation warning change to use make_scale + added valid data extrapolation with nearest_extrap Updated 11/2020: added model constituents from TPXO9-atlas-v3 Updated 10/2020: using argparse to set command line parameters Updated 08/2020: using builtin time operations. python3 regular expressions @@ -116,7 +120,7 @@ #-- PURPOSE: read ICESat-2 ocean surface height (ATL12) from NSIDC #-- compute tides at points and times using tidal model driver algorithms def compute_tides_ICESat2(tide_dir, FILE, TIDE_MODEL=None, METHOD='spline', - VERBOSE=False, MODE=0o775): + EXTRAPOLATE=False, VERBOSE=False, MODE=0o775): #-- select between tide models if (TIDE_MODEL == 'CATS0201'): grid_file = os.path.join(tide_dir,'cats0201_tmd','grid_CATS') @@ -471,23 +475,25 @@ def compute_tides_ICESat2(tide_dir, FILE, TIDE_MODEL=None, METHOD='spline', if model_format in ('OTIS','ATLAS'): amp,ph,D,c = extract_tidal_constants(val['longitude'], val['latitude'], grid_file, model_file, EPSG, TYPE=TYPE, - METHOD=METHOD, GRID=model_format) + METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE, GRID=model_format) deltat = np.zeros_like(tide_time) elif (model_format == 'netcdf'): amp,ph,D,c = extract_netcdf_constants(val['longitude'], val['latitude'], model_directory, grid_file, - model_files, TYPE=TYPE, METHOD=METHOD, SCALE=SCALE) + model_files, TYPE=TYPE, METHOD=METHOD, + EXTRAPOLATE=EXTRAPOLATE, SCALE=SCALE) deltat = np.zeros_like(tide_time) elif (model_format == 'GOT'): amp,ph = extract_GOT_constants(val['longitude'], val['latitude'], - model_directory, model_files, METHOD=METHOD, SCALE=SCALE) + model_directory, model_files, 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, tide_time) elif (model_format == 'FES'): amp,ph = extract_FES_constants(val['longitude'], val['latitude'], model_directory, model_files, TYPE=TYPE, VERSION=TIDE_MODEL, - 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, tide_time) @@ -825,6 +831,10 @@ def main(): metavar='METHOD', type=str, default='spline', choices=('spline','linear','nearest','bilinear'), help='Spatial interpolation method') + #-- extrapolate with nearest-neighbors + parser.add_argument('--extrapolate','-E', + default=False, action='store_true', + help='Extrapolate with nearest-neighbors') #-- verbosity settings #-- verbose will output information about each output file parser.add_argument('--verbose','-V', @@ -839,7 +849,8 @@ def main(): #-- run for each input ATL12 file for FILE in args.infile: compute_tides_ICESat2(args.directory, FILE, TIDE_MODEL=args.tide, - METHOD=args.interpolate, VERBOSE=args.verbose, MODE=args.mode) + METHOD=args.interpolate, EXTRAPOLATE=args.extrapolate, + VERBOSE=args.verbose, MODE=args.mode) #-- run main program if __name__ == '__main__': diff --git a/scripts/compute_tides_ICESat_GLA12.py b/scripts/compute_tides_ICESat_GLA12.py index 8901e0dd..3b2d8b28 100644 --- a/scripts/compute_tides_ICESat_GLA12.py +++ b/scripts/compute_tides_ICESat_GLA12.py @@ -41,6 +41,7 @@ linear nearest bilinear + -E X, --extrapolate X: Extrapolate with nearest-neighbors -M X, --mode X: Permission mode of directories and files created -V, --verbose: Output information about each created file @@ -65,15 +66,18 @@ infer_minor_corrections.py: return corrections for minor constituents load_constituent.py: loads parameters for a given tidal constituent load_nodal_corrections.py: load the nodal corrections for tidal constituents - predict_tide_drift.py: predict tidal elevations using harmonic constants read_tide_model.py: extract tidal harmonic constants from OTIS tide models 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 + predict_tide_drift.py: predict tidal elevations using harmonic constants UPDATE HISTORY: Updated 12/2020: updated for public release H5py deprecation warning change to use make_scale and not create_scale + added valid data extrapolation with nearest_extrap Updated 11/2020: added model constituents from TPXO9-atlas-v3 Updated 10/2020: using argparse to set command line parameters Updated 08/2020: using builtin time operations. python3 regular expressions @@ -112,7 +116,7 @@ #-- PURPOSE: read ICESat ice sheet HDF5 elevation data (GLAH12) from NSIDC #-- compute tides at points and times using tidal model driver algorithms def compute_tides_ICESat(tide_dir, FILE, TIDE_MODEL=None, METHOD='spline', - VERBOSE=False, MODE=0o775): + EXTRAPOLATE=False, VERBOSE=False, MODE=0o775): #-- select between tide models if (TIDE_MODEL == 'CATS0201'): grid_file = os.path.join(tide_dir,'cats0201_tmd','grid_CATS') @@ -429,7 +433,7 @@ def compute_tides_ICESat(tide_dir, FILE, TIDE_MODEL=None, METHOD='spline', #-- new reference orbit ground track file is obtained.) #-- INST: Instance number (increments every time the satellite enters a #-- different reference orbit) - #-- CYCL: Cycle of reference orbit for this phase + #-- CYCL: Cycle of reference orbit for this phase #-- TRK: Track within reference orbit #-- SEG: Segment of orbit #-- GRAN: Granule version number @@ -465,21 +469,24 @@ def compute_tides_ICESat(tide_dir, FILE, TIDE_MODEL=None, METHOD='spline', #-- read tidal constants and interpolate to grid points if model_format in ('OTIS','ATLAS'): amp,ph,D,c = extract_tidal_constants(lon_40HZ, lat_40HZ, grid_file, - model_file, EPSG, TYPE=TYPE, METHOD=METHOD, GRID=model_format) + model_file, EPSG, TYPE=TYPE, METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE, + GRID=model_format) deltat = np.zeros_like(tide_time) elif (model_format == 'netcdf'): amp,ph,D,c = extract_netcdf_constants(lon_40HZ, lat_40HZ, model_directory, - grid_file, model_files, TYPE=TYPE, METHOD=METHOD, SCALE=SCALE) + grid_file, model_files, TYPE=TYPE, METHOD=METHOD, + EXTRAPOLATE=EXTRAPOLATE, SCALE=SCALE) deltat = np.zeros_like(tide_time) elif (model_format == 'GOT'): amp,ph = extract_GOT_constants(lon_40HZ, lat_40HZ, model_directory, - model_files, METHOD=METHOD, SCALE=SCALE) + model_files, 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, tide_time) elif (model_format == 'FES'): amp,ph = extract_FES_constants(lon_40HZ, lat_40HZ, model_directory, - model_files, TYPE=TYPE, VERSION=TIDE_MODEL, METHOD=METHOD, SCALE=SCALE) + model_files, TYPE=TYPE, VERSION=TIDE_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, tide_time) @@ -664,7 +671,7 @@ def HDF5_GLA12_tide_write(IS_gla12_tide, IS_gla12_attrs, #-- add HDF5 group attributes for att_name,att_val in IS_gla12_attrs['Data_40HZ'][group].items(): if not isinstance(att_val,dict): - fileID['Data_40HZ'][group].attrs[att_name] = att_val + fileID['Data_40HZ'][group].attrs[att_name] = att_val #-- for each variable in the group for key,val in IS_gla12_tide['Data_40HZ'][group].items(): fillvalue = FILL_VALUE['Data_40HZ'][group][key] @@ -724,6 +731,10 @@ def main(): metavar='METHOD', type=str, default='spline', choices=('spline','linear','nearest','bilinear'), help='Spatial interpolation method') + #-- extrapolate with nearest-neighbors + parser.add_argument('--extrapolate','-E', + default=False, action='store_true', + help='Extrapolate with nearest-neighbors') #-- verbosity settings #-- verbose will output information about each output file parser.add_argument('--verbose','-V', @@ -738,7 +749,8 @@ def main(): #-- run for each input GLA12 file for FILE in args.infile: compute_tides_ICESat(args.directory, FILE, TIDE_MODEL=args.tide, - METHOD=args.interpolate, VERBOSE=args.verbose, MODE=args.mode) + METHOD=args.interpolate, EXTRAPOLATE=args.extrapolate, + VERBOSE=args.verbose, MODE=args.mode) #-- run main program if __name__ == '__main__': diff --git a/scripts/compute_tides_icebridge_data.py b/scripts/compute_tides_icebridge_data.py index 0d16e36b..6b1e79a5 100644 --- a/scripts/compute_tides_icebridge_data.py +++ b/scripts/compute_tides_icebridge_data.py @@ -1,7 +1,7 @@ #!/usr/bin/env python u""" compute_tides_icebridge_data.py -Written by Tyler Sutterley (11/2020) +Written by Tyler Sutterley (12/2020) Calculates tidal elevations for correcting Operation IceBridge elevation data Uses OTIS format tidal solutions provided by Ohio State University and ESR @@ -43,6 +43,7 @@ linear nearest bilinear + -E X, --extrapolate X: Extrapolate with nearest-neighbors -M X, --mode X: Permission mode of directories and files created -V, --verbose: Output information about each created file @@ -69,14 +70,17 @@ infer_minor_corrections.py: return corrections for minor constituents load_constituent.py: loads parameters for a given tidal constituent load_nodal_corrections.py: load the nodal corrections for tidal constituents - predict_tide_drift.py: predict tidal elevations using harmonic constants read_tide_model.py: extract tidal harmonic constants from OTIS tide models 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 + predict_tide_drift.py: predict tidal elevations using harmonic constants read_ATM1b_QFIT_binary.py: read ATM1b QFIT binary files (NSIDC version 1) UPDATE HISTORY: + Updated 12/2020: added valid data extrapolation with nearest_extrap Updated 11/2020: added model constituents from TPXO9-atlas-v3 Updated 10/2020: using argparse to set command line parameters Updated 09/2020: output ocean and load tide as tide_ocean and tide_load @@ -411,7 +415,7 @@ def read_LVIS_HDF5_file(input_file, input_subsetter): #-- PURPOSE: read Operation IceBridge data from NSIDC #-- compute tides at points and times using tidal model driver algorithms def compute_tides_icebridge_data(tide_dir, arg, TIDE_MODEL, - METHOD='spline', VERBOSE=False, MODE=0o775): + METHOD='spline', EXTRAPOLATE=False, VERBOSE=False, MODE=0o775): #-- extract file name and subsetter indices lists match_object = re.match(r'(.*?)(\[(.*?)\])?$',arg) @@ -759,23 +763,24 @@ def compute_tides_icebridge_data(tide_dir, arg, TIDE_MODEL, if model_format in ('OTIS','ATLAS'): amp,ph,D,c = extract_tidal_constants(dinput['lon'], dinput['lat'], grid_file, model_file, EPSG, TYPE=TYPE, METHOD=METHOD, - GRID=model_format) + EXTRAPOLATE=EXTRAPOLATE, GRID=model_format) deltat = np.zeros_like(t) elif model_format in ('netcdf'): amp,ph,D,c = extract_netcdf_constants(dinput['lon'], dinput['lat'], model_directory, grid_file, model_files, TYPE=TYPE, METHOD=METHOD, - SCALE=SCALE) + EXTRAPOLATE=EXTRAPOLATE, SCALE=SCALE) deltat = np.zeros_like(t) elif (model_format == 'GOT'): amp,ph = extract_GOT_constants(dinput['lon'], dinput['lat'], - model_directory, model_files, METHOD=METHOD, SCALE=SCALE) + model_directory, model_files, METHOD=METHOD, + EXTRAPOLATE=EXTRAPOLATE, SCALE=SCALE) #-- interpolate delta times from calendar dates to tide time delta_file = get_data_path(['data','merged_deltat.data']) deltat = calc_delta_time(delta_file, t) elif (model_format == 'FES'): amp,ph = extract_FES_constants(dinput['lon'], dinput['lat'], model_directory, model_files, TYPE=TYPE, VERSION=TIDE_MODEL, - METHOD=METHOD, SCALE=SCALE) + METHOD=METHOD, EXTRAPOLATE=EXTRAPOLATE, SCALE=SCALE) #-- interpolate delta times from calendar dates to tide time delta_file = get_data_path(['data','merged_deltat.data']) deltat = calc_delta_time(delta_file, t) @@ -923,6 +928,10 @@ def main(): metavar='METHOD', type=str, default='spline', choices=('spline','linear','nearest','bilinear'), help='Spatial interpolation method') + #-- extrapolate with nearest-neighbors + parser.add_argument('--extrapolate','-E', + default=False, action='store_true', + help='Extrapolate with nearest-neighbors') #-- verbosity settings #-- verbose will output information about each output file parser.add_argument('--verbose','-V', @@ -937,7 +946,8 @@ def main(): #-- run for each input Operation IceBridge file for arg in args.infile: compute_tides_icebridge_data(args.directory, arg, TIDE_MODEL=args.tide, - METHOD=args.interpolate, VERBOSE=args.verbose, MODE=args.mode) + METHOD=args.interpolate, EXTRAPOLATE=args.extrapolate, + VERBOSE=args.verbose, MODE=args.mode) #-- run main program if __name__ == '__main__': diff --git a/setup.py b/setup.py index bfc0ea0c..fb153bca 100644 --- a/setup.py +++ b/setup.py @@ -53,7 +53,6 @@ def check_output(cmd): 'Topic :: Scientific/Engineering :: Physics', 'License :: OSI Approved :: MIT License', 'Programming Language :: Python :: 3', - 'Programming Language :: Python :: 3.5', 'Programming Language :: Python :: 3.6', 'Programming Language :: Python :: 3.7', 'Programming Language :: Python :: 3.8', diff --git a/test/test_atlas_read.py b/test/test_atlas_read.py index 29a36d55..1cba07e8 100644 --- a/test/test_atlas_read.py +++ b/test/test_atlas_read.py @@ -97,12 +97,13 @@ def test_download_TPXO9_v2(aws_access_key_id,aws_secret_access_key,aws_region_na @pytest.mark.parametrize("MODEL", ['TPXO8-atlas','TPXO9-atlas-v2']) #-- parameterize interpolation method @pytest.mark.parametrize("METHOD", ['spline','nearest','bilinear']) +@pytest.mark.parametrize("EXTRAPOLATE", [True]) #-- PURPOSE: test the tide correction wrapper function -def test_Ross_Ice_Shelf(MODEL, METHOD): +def test_Ross_Ice_Shelf(MODEL, METHOD, EXTRAPOLATE): #-- create an image around the Ross Ice Shelf xlimits = np.array([-740000,520000]) ylimits = np.array([-1430000,-300000]) - spacing = np.array([5e3,-5e3]) + spacing = np.array([10e3,-10e3]) #-- x and y coordinates x = np.arange(xlimits[0],xlimits[1]+spacing[0],spacing[0]) y = np.arange(ylimits[1],ylimits[0]+spacing[1],spacing[1]) @@ -115,5 +116,6 @@ def test_Ross_Ice_Shelf(MODEL, METHOD): #-- calculate tide map tide = pyTMD.compute_tide_corrections(xgrid, ygrid, delta_time, DIRECTORY=filepath, MODEL=MODEL, EPOCH=(2000,1,1,0,0,0), - TYPE='grid', TIME='TAI', EPSG=3031, METHOD=METHOD) + TYPE='grid', TIME='TAI', EPSG=3031, METHOD=METHOD, + EXTRAPOLATE=EXTRAPOLATE) assert np.any(tide) diff --git a/test/test_perth3_read.py b/test/test_perth3_read.py index 73549aeb..bf96d35f 100644 --- a/test/test_perth3_read.py +++ b/test/test_perth3_read.py @@ -133,12 +133,13 @@ def test_verify_GOT47(METHOD): #-- parameterize interpolation method @pytest.mark.parametrize("METHOD", ['spline','nearest','bilinear']) +@pytest.mark.parametrize("EXTRAPOLATE", [True]) #-- PURPOSE: test the tide correction wrapper function -def test_Ross_Ice_Shelf(METHOD): +def test_Ross_Ice_Shelf(METHOD, EXTRAPOLATE): #-- create an image around the Ross Ice Shelf xlimits = np.array([-740000,520000]) ylimits = np.array([-1430000,-300000]) - spacing = np.array([5e3,-5e3]) + spacing = np.array([10e3,-10e3]) #-- x and y coordinates x = np.arange(xlimits[0],xlimits[1]+spacing[0],spacing[0]) y = np.arange(ylimits[1],ylimits[0]+spacing[1],spacing[1]) @@ -151,5 +152,6 @@ def test_Ross_Ice_Shelf(METHOD): #-- calculate tide map tide = pyTMD.compute_tide_corrections(xgrid, ygrid, delta_time, DIRECTORY=filepath, MODEL='GOT4.7', EPOCH=(2018,1,1,0,0,0), - TYPE='grid', TIME='GPS', EPSG=3031, METHOD=METHOD) + TYPE='grid', TIME='GPS', EPSG=3031, METHOD=METHOD, + EXTRAPOLATE=EXTRAPOLATE) assert np.any(tide) diff --git a/test/test_time.py b/test/test_time.py index c0a4393a..7f4e5479 100644 --- a/test/test_time.py +++ b/test/test_time.py @@ -83,6 +83,11 @@ def test_decimal_dates(YEAR,MONTH): assert (np.floor(minute_temp) == MINUTE) assert (np.abs(second - SECOND) < eps) +#-- PURPOSE: test UNIX time +def test_unix_time(): + UNIX = pyTMD.utilities.get_unix_time('2018-01-01 00:00:00') + assert (UNIX == 1514764800) + #-- PURPOSE: test parsing time strings def test_parse_date_string(): #-- time string for Modified Julian Days diff --git a/test/test_utilities.py b/test/test_utilities.py new file mode 100644 index 00000000..f56a8370 --- /dev/null +++ b/test/test_utilities.py @@ -0,0 +1,33 @@ +#!/usr/bin/env python +u""" +test_utilities.py (12/2020) +Verify file utility functions +""" +import io +import gzip +import pytest +import warnings +import posixpath +import pyTMD.utilities + +def test_hash(): + # get hash of compressed file + ocean_pole_tide_file = pyTMD.utilities.get_data_path(['data', + 'opoleloadcoefcmcor.txt.gz']) + TEST = pyTMD.utilities.get_hash(ocean_pole_tide_file) + assert (TEST == '9c66edc2d0fbf627e7ae1cb923a9f0e5') + # get hash of uncompressed file + with gzip.open(ocean_pole_tide_file) as fid: + TEST = pyTMD.utilities.get_hash(io.BytesIO(fid.read())) + assert (TEST == 'cea08f83d613ed8e1a81f3b3a9453721') + +def test_urlsplit(): + HOST = ['https://cddis.nasa.gov','archive','products','iers','deltat.data'] + url = posixpath.join(*HOST) + TEST = pyTMD.utilities.url_split(url) + assert (HOST == list(TEST)) + +def test_roman(): + for s,i in zip(['MMXVIII','MCMXVI','DXCIV'],[2018,1916,594]): + TEST = pyTMD.utilities.roman_to_int(s) + assert (TEST == i)