feat: added new Gr1km-v2 1km Greenland model from ESR (#59)

feat: added 3413 projection for new 1km Greenland model from ESR

feat: update notebooks for Gr1km-v2

fix: tidal currents for bilinear interpolation to address part of #57

fix: add warning for tide models being entered as string to address #58

docs: update documentation for new model

docs: add notes to interpolation with warnings

fix: check for nan points when reading OTIS models

test: update list of invalid stations in Arctic

Dockerfile

@@ -37,6 +37,7 @@ RUN wget -q http://zlib.net/zlib-${ZLIB_VERSION}.tar.gz && \
     ./configure --prefix=/usr/local && \
     make --quiet --jobs=${JOBS} && \
     make --quiet install && \
+    make clean && \
     cd contrib/minizip && \
     autoreconf -fvi && \
     ./configure --quiet --prefix=/usr/local/ && \

doc/source/getting_started/Getting-Started.rst

@@ -42,10 +42,11 @@ This structure was chosen based on the different formats of each tide model.
     * `CATS2008 <https://www.usap-dc.org/view/dataset/601235>`_: ``<path_to_tide_models>/CATS2008/``
     * CATS2008_load: ``<path_to_tide_models>/CATS2008a_SPOTL_Load/``
 
-- Arctic Ocean Tidal Simulations [Padman2004]_
+- Arctic Ocean and Greenland Coast Tidal Simulations [Padman2004]_
     * `AODTM-5 <https://arcticdata.io/catalog/view/doi:10.18739/A2901ZG3N>`_: ``<path_to_tide_models>/aodtm5_tmd/``
     * `AOTIM-5 <https://arcticdata.io/catalog/view/doi:10.18739/A2S17SS80>`_: ``<path_to_tide_models>/aotim5_tmd/``
     * `AOTIM-5-2018 <https://arcticdata.io/catalog/view/doi:10.18739/A21R6N14K>`_: ``<path_to_tide_models>/Arc5km2018/``
+    * Gr1km-v2: ``<path_to_tide_models>/greenlandTMD_v2/``
 
 - TOPEX/POSEIDON global tide models [Egbert2002]_
     * `TPXO9-atlas <https://www.tpxo.net/tpxo-products-and-registration>`_: ``<path_to_tide_models>/TPXO9_atlas/``
@@ -148,7 +149,8 @@ There are options to use nearest and linear interpolators with the
 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.
 The default maximum extrapolation distance is 10 kilometers.
-This extrapolation cutoff can be set to any distance in kilometers.
+This default distance may not be a large enough extrapolation for some applications and models.
+The extrapolation cutoff can be set to any distance in kilometers, but should be used with caution in cases such as narrow fjords or ice sheet grounding zones [Padman2018]_.
 
 References
 ##########
@@ -161,6 +163,8 @@ References
 
 .. [Padman2008] L. Padman, S. Y. Erofeeva, and H. A. Fricker, "Improving Antarctic tide models by assimilation of ICESat laser altimetry over ice shelves", *Geophysical Research Letters*, 35, L22504, (2008). `doi: 10.1029/2008GL035592 <https://doi.org/10.1029/2008GL035592>`_
 
+.. [Padman2018] L. Padman, M. R. Siegfried, and H. A. Fricker, "Ocean Tide Influences on the Antarctic and Greenland Ice Sheets", *Reviews of Geophysics*, 56, (2018). `doi: 10.1002/2016RG000546 <https://doi.org/10.1002/2016RG000546>`_
+
 .. [Ray1999] R. D. Ray, "A Global Ocean Tide Model From TOPEX/POSEIDON Altimetry: GOT99.2", *NASA Technical Memorandum*, `NASA/TM--1999-209478 <https://ntrs.nasa.gov/search.jsp?R=19990089548>`_.
 
 .. [Stammer2014] D. Stammer et al., "Accuracy assessment of global barotropic ocean tide models", *Reviews of Geophysics*, 52, 243--282, (2014). `doi: 10.1002/2014RG000450 <https://doi.org/10.1002/2014RG000450>`_

doc/source/user_guide/compute_tidal_currents.md

@@ -31,6 +31,7 @@ python compute_tidal_currents.py --directory <path_to_directory> --tide <model>
     * `'AODTM-5'`
     * `'AOTIM-5'`
     * `'AOTIM-5-2018'`
+    * `'Gr1km-v2'`
     * `'FES2014'`
 - `--format X`: input and output data format
     * `'csv'` (default)

doc/source/user_guide/compute_tidal_elevations.md

@@ -34,6 +34,7 @@ python compute_tidal_elevations.py --directory <path_to_directory> --tide <model
     * `'AODTM-5'`
     * `'AOTIM-5'`
     * `'AOTIM-5-2018'`
+    * `'Gr1km-v2'`
     * `'GOT4.7'`
     * `'GOT4.7_load'`
     * `'GOT4.8'`

doc/source/user_guide/compute_tides_ICESat2_ATL03.md

@@ -34,6 +34,7 @@ python compute_tides_ICESat2_ATL03.py --directory <path_to_directory> --tide <mo
     * `'AODTM-5'`
     * `'AOTIM-5'`
     * `'AOTIM-5-2018'`
+    * `'Gr1km-v2'`
     * `'GOT4.7'`
     * `'GOT4.7_load'`
     * `'GOT4.8'`

doc/source/user_guide/compute_tides_ICESat2_ATL06.md

@@ -32,6 +32,7 @@ python compute_tides_ICESat2_ATL06.py --directory <path_to_directory> --tide <mo
     * `'AODTM-5'`
     * `'AOTIM-5'`
     * `'AOTIM-5-2018'`
+    * `'Gr1km-v2'`
     * `'GOT4.7'`
     * `'GOT4.7_load'`
     * `'GOT4.8'`

doc/source/user_guide/compute_tides_ICESat2_ATL07.md

@@ -32,6 +32,7 @@ python compute_tides_ICESat2_ATL07.py --directory <path_to_directory> --tide <mo
     * `'AODTM-5'`
     * `'AOTIM-5'`
     * `'AOTIM-5-2018'`
+    * `'Gr1km-v2'`
     * `'GOT4.7'`
     * `'GOT4.7_load'`
     * `'GOT4.8'`

doc/source/user_guide/compute_tides_ICESat2_ATL11.md

@@ -33,6 +33,7 @@ python compute_tides_ICESat2_ATL11.py --directory <path_to_directory> --tide <mo
     * `'AODTM-5'`
     * `'AOTIM-5'`
     * `'AOTIM-5-2018'`
+    * `'Gr1km-v2'`
     * `'GOT4.7'`
     * `'GOT4.7_load'`
     * `'GOT4.8'`

doc/source/user_guide/compute_tides_ICESat2_ATL12.md

@@ -32,6 +32,7 @@ python compute_tides_ICESat2_ATL12.py --directory <path_to_directory> --tide <mo
     * `'AODTM-5'`
     * `'AOTIM-5'`
     * `'AOTIM-5-2018'`
+    * `'Gr1km-v2'`
     * `'GOT4.7'`
     * `'GOT4.7_load'`
     * `'GOT4.8'`

doc/source/user_guide/compute_tides_ICESat_GLA12.md

@@ -32,6 +32,7 @@ python compute_tides_ICESat_GLA12.py --directory <path_to_directory> --tide <mod
     * `'AODTM-5'`
     * `'AOTIM-5'`
     * `'AOTIM-5-2018'`
+    * `'Gr1km-v2'`
     * `'GOT4.7'`
     * `'GOT4.7_load'`
     * `'GOT4.8'`

doc/source/user_guide/compute_tides_icebridge_data.md

@@ -32,6 +32,7 @@ python compute_tides_icebridge_data.py --directory <path_to_directory> --tide <m
     * `'AODTM-5'`
     * `'AOTIM-5'`
     * `'AOTIM-5-2018'`
+    * `'Gr1km-v2'`
     * `'GOT4.7'`
     * `'GOT4.7_load'`
     * `'GOT4.8'`

notebooks/Check Tide Map.ipynb

@@ -92,7 +92,7 @@
     "# dropdown menu for setting tide model\n",
     "model_list = ['CATS0201','CATS2008','TPXO9-atlas','TPXO9-atlas-v2',\n",
     "    'TPXO9-atlas-v3','TPXO9-atlas-v4','TPXO9.1','TPXO8-atlas','TPXO7.2',\n",
-    "    'AODTM-5','AOTIM-5','AOTIM-5-2018',\n",
+    "    'AODTM-5','AOTIM-5','AOTIM-5-2018','Gr1km-v2',\n",
     "    'GOT4.7','GOT4.8','GOT4.10','FES2014']\n",
     "modelDropdown = widgets.Dropdown(\n",
     "    options=model_list,\n",
@@ -211,6 +211,10 @@
     "    grid_file = os.path.join(tide_dir,'Arc5km2018','grid_Arc5km2018')\n",
     "    model_format = 'OTIS'\n",
     "    EPSG = 'PSNorth'\n",
+    "elif (TIDE_MODEL == 'Gr1km-v2'):\n",
+    "    grid_file = os.path.join(tide_dir,'greenlandTMD_v2','grid_Greenland8.v2')\n",
+    "    model_format = 'OTIS'\n",
+    "    EPSG = '3413'\n",
     "elif (TIDE_MODEL == 'GOT4.7'):\n",
     "    model_directory = os.path.join(tide_dir,'GOT4.7','grids_oceantide')\n",
     "    model_files = ['q1.d.gz','o1.d.gz','p1.d.gz','k1.d.gz','n2.d.gz',\n",
@@ -336,4 +340,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 4
-}
+}

notebooks/Plot Arctic Ocean Map.ipynb

@@ -117,8 +117,8 @@
     "#-- dropdown menu for setting tide model\n",
     "model_list = ['TPXO9-atlas','TPXO9-atlas-v2','TPXO9-atlas-v3',\n",
     "    'TPXO9-atlas-v4','TPXO9.1','TPXO8-atlas','TPXO7.2',\n",
-    "    'AODTM-5','AOTIM-5','AOTIM-5-2018','GOT4.7','GOT4.8',\n",
-    "    'GOT4.10','FES2014']\n",
+    "    'AODTM-5','AOTIM-5','AOTIM-5-2018','Gr1km-v2',\n",
+    "    'GOT4.7','GOT4.8','GOT4.10','FES2014']\n",
     "modelDropdown = widgets.Dropdown(\n",
     "    options=model_list,\n",
     "    value='AOTIM-5-2018',\n",
@@ -257,6 +257,13 @@
     "    model_format = 'OTIS'\n",
     "    EPSG = 'PSNorth'\n",
     "    TYPE = 'z'\n",
+    "elif (MODEL == 'Gr1km-v2'):\n",
+    "    grid_file = os.path.join(tide_dir,'greenlandTMD_v2','grid_Greenland8.v2')\n",
+    "    model_file = os.path.join(tide_dir,'greenlandTMD_v2','h_Greenland8.v2')\n",
+    "    reference = 'https://www.esr.org/research/polar-tide-models/'\n",
+    "    model_format = 'OTIS'\n",
+    "    EPSG = '3413'\n",
+    "    TYPE = 'z'\n",
     "elif (MODEL == 'GOT4.7'):\n",
     "    model_directory = os.path.join(tide_dir,'GOT4.7','grids_oceantide')\n",
     "    model_files = ['q1.d.gz','o1.d.gz','p1.d.gz','k1.d.gz','n2.d.gz',\n",
@@ -496,4 +503,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 4
-}
+}

notebooks/Plot Tide Forecasts.ipynb

@@ -102,7 +102,7 @@
     "#-- dropdown menu for setting tide model\n",
     "model_list = ['CATS0201','CATS2008','TPXO9-atlas','TPXO9-atlas-v2',\n",
     "    'TPXO9-atlas-v3','TPXO9-atlas-v4','TPXO9.1','TPXO8-atlas','TPXO7.2',\n",
-    "    'AODTM-5','AOTIM-5','AOTIM-5-2018',\n",
+    "    'AODTM-5','AOTIM-5','AOTIM-5-2018','Gr1km-v2',\n",
     "    'GOT4.7','GOT4.8','GOT4.10','FES2014']\n",
     "modelDropdown = widgets.Dropdown(\n",
     "    options=model_list,\n",
@@ -297,6 +297,12 @@
     "    model_format = 'OTIS'\n",
     "    EPSG = 'PSNorth'\n",
     "    TYPE = 'z'\n",
+    "elif (TIDE_MODEL == 'Gr1km-v2'):\n",
+    "    grid_file = os.path.join(tide_dir,'greenlandTMD_v2','grid_Greenland8.v2')\n",
+    "    model_file = os.path.join(tide_dir,'greenlandTMD_v2','h_Greenland8.v2')\n",
+    "    model_format = 'OTIS'\n",
+    "    EPSG = '3413'\n",
+    "    TYPE = 'z'\n",
     "elif (TIDE_MODEL == 'GOT4.7'):\n",
     "    model_directory = os.path.join(tide_dir,'GOT4.7','grids_oceantide')\n",
     "    model_files = ['q1.d.gz','o1.d.gz','p1.d.gz','k1.d.gz','n2.d.gz',\n",
@@ -428,9 +434,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.5-final"
+   "version": "3.8.5"
   }
  },
  "nbformat": 4,
  "nbformat_minor": 4
-}
+}

pyTMD/compute_tide_corrections.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 compute_tide_corrections.py
-Written by Tyler Sutterley (05/2021)
+Written by Tyler Sutterley (06/2021)
 Calculates tidal elevations for correcting elevation or imagery data
 
 Uses OTIS format tidal solutions provided by Ohio State University and ESR
@@ -69,6 +69,8 @@
     nearest_extrap.py: nearest-neighbor extrapolation of data to coordinates
 
 UPDATE HISTORY:
+    Updated 06/2021: added new Gr1km-v2 1km Greenland model from ESR
+        add try/except for input projection strings
     Updated 05/2021: added option for extrapolation cutoff in kilometers
     Updated 03/2021: added TPXO9-atlas-v4 in binary OTIS format
         simplified netcdf inputs to be similar to binary OTIS read program
@@ -259,6 +261,12 @@ def compute_tide_corrections(x, y, delta_time, DIRECTORY=None, MODEL=None,
         model_format = 'OTIS'
         model_EPSG = 'PSNorth'
         model_type = 'z'
+    elif (MODEL == 'Gr1km-v2'):
+        grid_file = os.path.join(DIRECTORY,'greenlandTMD_v2','grid_Greenland8.v2')
+        model_file = os.path.join(DIRECTORY,'greenlandTMD_v2','h_Greenland8.v2')
+        model_format = 'OTIS'
+        model_EPSG = '3413'
+        model_type = 'z'
     elif (MODEL == 'GOT4.7'):
         model_directory = os.path.join(DIRECTORY,'GOT4.7','grids_oceantide')
         model_files = ['q1.d.gz','o1.d.gz','p1.d.gz','k1.d.gz','n2.d.gz',
@@ -348,7 +356,12 @@ def compute_tide_corrections(x, y, delta_time, DIRECTORY=None, MODEL=None,
         raise Exception("Unlisted tide model")
 
     #-- converting x,y from EPSG to latitude/longitude
-    crs1 = pyproj.CRS.from_string("epsg:{0:d}".format(EPSG))
+    try:
+        #-- EPSG projection code string or int
+        crs1 = pyproj.CRS.from_string("epsg:{0:d}".format(int(EPSG)))
+    except (ValueError,pyproj.exceptions.CRSError):
+        #-- Projection SRS string
+        crs1 = pyproj.CRS.from_string(EPSG)
     crs2 = pyproj.CRS.from_string("epsg:{0:d}".format(4326))
     transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True)
     lon,lat = transformer.transform(x.flatten(), y.flatten())

pyTMD/convert_ll_xy.py

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 u"""
-convert_ll_xy.py (08/2020)
+convert_ll_xy.py (06/2021)
 Wrapper function to convert lat/lon points to and from projected coordinates
 
 CALLING SEQUENCE:
@@ -29,6 +29,7 @@
         https://pyproj4.github.io/pyproj/
 
 UPDATE HISTORY:
+    Updated 06/2021: added 3413 for new 1km Greenland model from ESR
     Updated 08/2020: using conversion protocols following pyproj-2 updates
         https://pyproj4.github.io/pyproj/stable/gotchas.html
     Updated 07/2020: added function docstrings. changed function name
@@ -62,6 +63,7 @@ def convert_ll_xy(i1,i2,PROJ,BF,EPSG=4326):
     #-- python dictionary with conversion functions
     conversion_functions = {}
     conversion_functions['3031'] = xy_ll_EPSG3031
+    conversion_functions['3413'] = xy_ll_EPSG3413
     conversion_functions['CATS2008'] = xy_ll_CATS2008
     conversion_functions['3976'] = xy_ll_EPSG3976
     conversion_functions['PSNorth'] = xy_ll_PSNorth
@@ -90,6 +92,23 @@ def xy_ll_EPSG3031(i1,i2,BF,EPSG=4326):
     #-- return the output variables
     return transformer.transform(i1, i2, direction=direction)
 
+#-- wrapper function for models in EPSG 3413 (Sea Ice Polar Stereographic North)
+def xy_ll_EPSG3413(i1,i2,BF,EPSG=4326):
+    #-- projections for converting from input EPSG (default latitude/longitude)
+    crs1 = pyproj.CRS.from_string("epsg:{0:d}".format(EPSG))
+    crs2 = pyproj.CRS.from_user_input({'proj':'stere','lat_0':90,'lat_ts':70,
+        'lon_0':-45,'x_0':0.,'y_0':0.,'ellps': 'WGS84','datum': 'WGS84',
+        'units':'km'})
+    transformer = pyproj.Transformer.from_crs(crs1, crs2, always_xy=True)
+    #-- convert lat/lon to Polar-Stereographic x/y
+    if (BF.upper() == 'F'):
+        direction = pyproj.enums.TransformDirection.FORWARD
+    #-- convert Polar-Stereographic x/y to lat/lon
+    elif (BF.upper() == 'B'):
+        direction = pyproj.enums.TransformDirection.INVERSE
+    #-- return the output variables
+    return transformer.transform(i1, i2, direction=direction)
+
 #-- wrapper function for CATS2008 tide models
 def xy_ll_CATS2008(i1,i2,BF,EPSG=4326):
     #-- projections for converting from input EPSG (default latitude/longitude)

pyTMD/read_FES_model.py

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 u"""
-read_FES_model.py (05/2021)
+read_FES_model.py (06/2021)
 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
@@ -53,6 +53,7 @@
     nearest_extrap.py: nearest-neighbor extrapolation of data to coordinates
 
 UPDATE HISTORY:
+    Updated 06/2021: add warning for tide models being entered as string
     Updated 05/2021: added option for extrapolation cutoff in kilometers
     Updated 03/2021: add extrapolation check where there are no invalid points
         prevent ComplexWarning for fill values when calculating amplitudes
@@ -71,6 +72,7 @@
 import gzip
 import uuid
 import netCDF4
+import warnings
 import numpy as np
 import scipy.interpolate
 from pyTMD.bilinear_interp import bilinear_interp
@@ -118,6 +120,11 @@ def extract_FES_constants(ilon, ilat, model_files, TYPE='z', VERSION=None,
     phase: phases of tidal constituents
     """
 
+    #-- raise warning if model files are entered as a string
+    if isinstance(model_files,str):
+        warnings.warn("Tide model is entered as a string")
+        model_files = [model_files]
+
     #-- adjust dimensions of input coordinates to be iterable
     ilon = np.atleast_1d(ilon)
     ilat = np.atleast_1d(ilat)

pyTMD/read_GOT_model.py

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 u"""
-read_GOT_model.py (05/2021)
+read_GOT_model.py (06/2021)
 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
@@ -39,6 +39,7 @@
     nearest_extrap.py: nearest-neighbor extrapolation of data to coordinates
 
 UPDATE HISTORY:
+    Updated 06/2021: add warning for tide models being entered as string
     Updated 05/2021: added option for extrapolation cutoff in kilometers
     Updated 03/2021: add extrapolation check where there are no invalid points
         prevent ComplexWarning for fill values when calculating amplitudes
@@ -66,6 +67,7 @@
 import os
 import re
 import gzip
+import warnings
 import numpy as np
 import scipy.interpolate
 from pyTMD.bilinear_interp import bilinear_interp
@@ -104,6 +106,11 @@ def extract_GOT_constants(ilon, ilat, model_files, METHOD=None,
     constituents: list of model constituents
     """
 
+    #-- raise warning if model files are entered as a string
+    if isinstance(model_files,str):
+        warnings.warn("Tide model is entered as a string")
+        model_files = [model_files]
+
     #-- adjust dimensions of input coordinates to be iterable
     ilon = np.atleast_1d(ilon)
     ilat = np.atleast_1d(ilat)