fix: add pole case in stereographic area scale calculation (#70)

feat: added generic list from Apache http server
fix: use getncattr to get attributes from netCDF4 files to prevent deprecation errors
tests: add wrap longitude and update time tests to use np.isclose

doc/source/getting_started/Citations.rst

@@ -52,6 +52,7 @@ This software is also dependent on other commonly used Python packages:
 - `cartopy: Python package designed for geospatial data processing <https://scitools.org.uk/cartopy/docs/latest/>`_
 - `ipywidgets: interactive HTML widgets for Jupyter notebooks and IPython <https://ipywidgets.readthedocs.io/en/latest/>`_
 - `ipyleaflet: Jupyter / Leaflet bridge enabling interactive maps <https://github.com/jupyter-widgets/ipyleaflet>`_
+- `setuptools_scm: manager for python package versions using scm metadata <https://pypi.org/project/setuptools-scm/1.9.0/>`_
 
 Credits
 #######

doc/source/user_guide/spatial.rst

@@ -259,9 +259,17 @@ General Methods
         ``lat``: array of latitudes in degrees
 
     Returns:
-    
+
         ``delta_h``: difference in elevation for two ellipsoids
 
+.. method:: pyTMD.spatial.wrap_longitudes(lon):
+
+    Wraps longitudes to range from -180 to +180
+
+    Arguments:
+
+        ``lon``: longitude
+
 .. method:: pyTMD.spatial.to_cartesian(lon,lat,a_axis=6378137.0,flat=1.0/298.257223563)
 
     Converts geodetic coordinates to Cartesian coordinates

doc/source/user_guide/utilities.rst

@@ -32,7 +32,7 @@ General Methods
     Platform independent file opener
 
     Arguments:
-    
+
         ``filename``: path to file
 
 
@@ -192,6 +192,35 @@ General Methods
         ``HOST``: remote http host
 
 
+.. method:: pyTMD.utilities.http_list(HOST,timeout=None,context=ssl.SSLContext(),parser=lxml.etree.HTMLParser(),format='%Y-%m-%d %H:%M',pattern='',sort=False)
+
+    List a directory on an Apache http Server
+
+    Arguments:
+
+        ``HOST``: remote http host path split as list
+
+    Keyword arguments:
+
+        ``timeout``: timeout in seconds for blocking operations
+
+        ``context``: SSL context for url opener object
+
+        ``parser``: HTML parser for lxml
+
+        ``format``: format for input time string
+
+        ``pattern``: regular expression pattern for reducing list
+
+        ``sort``: sort output list
+
+    Returns:
+
+        ``output``: list of items in a directory
+
+        ``mtimes``: list of last modification times for items in the directory
+
+
 .. method:: pyTMD.utilities.from_http(HOST,timeout=None,context=ssl.SSLContext(),local=None,hash='',chunk=16384,verbose=False,fid=sys.stdout,mode=0o775)
 
     Download a file from a http host

environment.yml

@@ -2,23 +2,23 @@ name: pyTMD
 channels:
   - conda-forge
 dependencies:
-  - python>=3.6
-  - setuptools_scm
-  - notebook
-  - ipywidgets
+  - cartopy
+  - gdal
+  - h5py
   - ipyleaflet
+  - ipywidgets
+  - lxml
+  - matplotlib
+  - netCDF4
+  - notebook
   - numpy
-  - scipy
+  - pip
   - pyproj
+  - python>=3.6
   - python-dateutil
-  - matplotlib
-  - netCDF4
-  - h5py
-  - gdal
-  - lxml
   - pyyaml
-  - cartopy
-  - pip
+  - scipy
+  - setuptools_scm
   - pip:
      - git+https://github.com/tsutterley/read-ICESat-2.git
      - git+https://github.com/tsutterley/read-ATM1b-QFIT-binary.git

notebooks/Check Tide Map.ipynb

@@ -71,6 +71,7 @@
     "import pyTMD.read_netcdf_model\n",
     "import pyTMD.read_GOT_model\n",
     "import pyTMD.read_FES_model\n",
+    "from pyTMD.spatial import wrap_longitudes\n",
     "# autoreload\n",
     "%load_ext autoreload\n",
     "%autoreload 2"
@@ -133,7 +134,7 @@
    "execution_count": null,
    "source": [
     "# default coordinates to use\n",
-    "LAT,LON = (32.93301304,242.7294513)\n",
+    "LAT,LON = (32.86710263,-117.25750387)\n",
     "m = leaflet.Map(center=(LAT,LON), zoom=12, basemap=leaflet.basemaps.Esri.WorldTopoMap)\n",
     "# add control for zoom\n",
     "zoom_slider = widgets.IntSlider(description='Zoom level:', min=0, max=15, value=7)\n",
@@ -153,14 +154,14 @@
     "# add function for setting marker text if location changed\n",
     "def set_marker_text(sender):\n",
     "    LAT,LON = marker.location\n",
-    "    markerText.value = '{0:0.8f},{1:0.8f}'.format(LAT,LON % 360)\n",
+    "    markerText.value = '{0:0.8f},{1:0.8f}'.format(LAT,wrap_longitudes(LON))\n",
     "# add function for setting map center if location changed\n",
     "def set_map_center(sender):\n",
     "    m.center = marker.location\n",
     "# add function for setting marker location if text changed\n",
     "def set_marker_location(sender):\n",
     "    LAT,LON = [float(i) for i in markerText.value.split(',')]\n",
-    "    marker.location = (LAT,LON % 360)\n",
+    "    marker.location = (LAT,LON)\n",
     "    \n",
     "# watch marker widgets for changes\n",
     "marker.observe(set_marker_text)\n",
@@ -181,7 +182,7 @@
     "# leaflet location\n",
     "LAT,LON = marker.location\n",
     "# verify longitudes\n",
-    "LON %= 360\n",
+    "LON = wrap_longitudes(LON)\n",
     "\n",
     "# get model parameters\n",
     "model = pyTMD.model(dirText.value, format=atlasDropdown.value,\n",

notebooks/Plot Tide Forecasts.ipynb

@@ -81,6 +81,7 @@
     "from pyTMD.read_netcdf_model import extract_netcdf_constants\n",
     "from pyTMD.read_GOT_model import extract_GOT_constants\n",
     "from pyTMD.read_FES_model import extract_FES_constants\n",
+    "from pyTMD.spatial import wrap_longitudes\n",
     "# autoreload\n",
     "%load_ext autoreload\n",
     "%autoreload 2"
@@ -151,7 +152,7 @@
    "execution_count": null,
    "source": [
     "# default coordinates to use\n",
-    "LAT,LON = (32.93301304,242.7294513)\n",
+    "LAT,LON = (32.86710263,-117.25750387)\n",
     "m = leaflet.Map(center=(LAT,LON), zoom=12, basemap=leaflet.basemaps.Esri.WorldTopoMap)\n",
     "# add control for zoom\n",
     "zoom_slider = widgets.IntSlider(description='Zoom level:', min=0, max=15, value=7)\n",
@@ -171,14 +172,14 @@
     "# add function for setting marker text if location changed\n",
     "def set_marker_text(sender):\n",
     "    LAT,LON = marker.location\n",
-    "    markerText.value = '{0:0.8f},{1:0.8f}'.format(LAT,LON % 360)\n",
+    "    markerText.value = '{0:0.8f},{1:0.8f}'.format(LAT,wrap_longitudes(LON))\n",
     "# add function for setting map center if location changed\n",
     "def set_map_center(sender):\n",
     "    m.center = marker.location\n",
     "# add function for setting marker location if text changed\n",
     "def set_marker_location(sender):\n",
     "    LAT,LON = [float(i) for i in markerText.value.split(',')]\n",
-    "    marker.location = (LAT,LON % 360)\n",
+    "    marker.location = (LAT,LON)\n",
     "    \n",
     "# watch marker widgets for changes\n",
     "marker.observe(set_marker_text)\n",
@@ -199,7 +200,7 @@
     "# leaflet location\n",
     "LAT,LON = marker.location\n",
     "# verify longitudes\n",
-    "LON %= 360\n",
+    "LON = wrap_longitudes(LON)\n",
     "\n",
     "# convert from calendar date to days relative to Jan 1, 1992 (48622 MJD)\n",
     "YMD = datepick.value\n",

pyTMD/spatial.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 u"""
 spatial.py
-Written by Tyler Sutterley (09/2021)
+Written by Tyler Sutterley (10/2021)
 
 Utilities for reading, writing and operating on spatial data
 
@@ -19,6 +19,7 @@
         https://github.com/yaml/pyyaml
 
 UPDATE HISTORY:
+    Updated 10/2021: add pole case in stereographic area scale calculation
     Updated 09/2021: can calculate height differences between ellipsoids
     Updated 07/2021: added function for determining input variable type
     Updated 03/2021: added polar stereographic area scale calculation
@@ -188,8 +189,8 @@ def from_netCDF4(filename, compression=None, verbose=False,
     for attr in ['title','description','projection']:
         #-- try getting the attribute
         try:
-            ncattr, = [s for s in dir(fileID) if re.match(attr,s,re.I)]
-            dinput['attributes'][attr] = getattr(fileID,ncattr)
+            ncattr, = [s for s in fileID.ncattrs() if re.match(attr,s,re.I)]
+            dinput['attributes'][attr] = fileID.getncattr(ncattr)
         except (ValueError,AttributeError):
             pass
     #-- list of attributes to attempt to retrieve from included variables
@@ -206,9 +207,10 @@ def from_netCDF4(filename, compression=None, verbose=False,
         for attr in attributes_list:
             #-- try getting the attribute
             try:
-                ncattr, = [s for s in dir(fileID) if re.match(attr,s,re.I)]
+                ncattr, = [s for s in fileID.variables[nc].ncattrs()
+                    if re.match(attr,s,re.I)]
                 dinput['attributes'][key][attr] = \
-                    getattr(fileID.variables[nc],ncattr)
+                    fileID.variables[nc].getncattr(ncattr)
             except (ValueError,AttributeError):
                 pass
     #-- convert to masked array if fill values
@@ -711,6 +713,17 @@ def compute_delta_h(a1, f1, a2, f2, lat):
     delta_h = -(delta_a*np.cos(phi)**2 + delta_b*np.sin(phi)**2)
     return delta_h
 
+def wrap_longitudes(lon):
+    """
+    Wraps longitudes to range from -180 to +180
+
+    Inputs:
+        lon: longitude (degrees east)
+    """
+    phi = np.arctan2(np.sin(lon*np.pi/180.0),np.cos(lon*np.pi/180.0))
+    #-- convert phi from radians to degrees
+    return phi*180.0/np.pi
+
 def to_cartesian(lon,lat,h=0.0,a_axis=6378137.0,flat=1.0/298.257223563):
     """
     Converts geodetic coordinates to Cartesian coordinates
@@ -835,6 +848,7 @@ def to_geodetic(x,y,z,a_axis=6378137.0,flat=1.0/298.257223563):
 def scale_areas(lat, flat=1.0/298.257223563, ref=70.0):
     """
     Calculates area scaling factors for a polar stereographic projection
+    including special case of at the exact pole
 
     Inputs:
         lat: latitude (degrees north)
@@ -869,8 +883,9 @@ def scale_areas(lat, flat=1.0/298.257223563, ref=70.0):
     mref = np.cos(theta_ref)/np.sqrt(1.0 - ecc2*np.sin(theta_ref)**2)
     tref = np.tan(np.pi/4.0 - theta_ref/2.0)/((1.0 - ecc*np.sin(theta_ref)) / \
         (1.0 + ecc*np.sin(theta_ref)))**(ecc/2.0)
-    #-- area scaling
+    #-- distance scaling
     k = (mref/m)*(t/tref)
-    #-- return the area scaling factors
-    scale = 1.0/(k**2)
+    kp = 0.5*mref*np.sqrt(((1.0+ecc)**(1.0+ecc))*((1.0-ecc)**(1.0-ecc)))/tref
+    #-- area scaling
+    scale = np.where(np.isclose(theta,np.pi/2.0),1.0/(kp**2),1.0/(k**2))
     return scale

pyTMD/utilities.py

@@ -1,13 +1,14 @@
 #!/usr/bin/env python
 u"""
 utilities.py
-Written by Tyler Sutterley (07/2021)
+Written by Tyler Sutterley (09/2021)
 Download and management utilities for syncing time and auxiliary files
 
 PYTHON DEPENDENCIES:
     lxml: processing XML and HTML in Python (https://pypi.python.org/pypi/lxml)
 
 UPDATE HISTORY:
+    Updated 09/2021: added generic list from Apache http server
     Updated 08/2021: added function to open a file path
     Updated 07/2021: add parser for converting file files to arguments
     Updated 03/2021: added sha1 option for retrieving file hashes
@@ -407,6 +408,60 @@ def check_connection(HOST):
     else:
         return True
 
+#-- PURPOSE: list a directory on an Apache http Server
+def http_list(HOST,timeout=None,context=ssl.SSLContext(),
+    parser=lxml.etree.HTMLParser(),format='%Y-%m-%d %H:%M',
+    pattern='',sort=False):
+    """
+    List a directory on an Apache http Server
+
+    Arguments
+    ---------
+    HOST: remote http host path split as list
+
+    Keyword arguments
+    -----------------
+    timeout: timeout in seconds for blocking operations
+    context: SSL context for url opener object
+    parser: HTML parser for lxml
+    format: format for input time string
+    pattern: regular expression pattern for reducing list
+    sort: sort output list
+
+    Returns
+    -------
+    output: list of items in a directory
+    mtimes: list of last modification times for items in the directory
+    """
+    #-- try listing from http
+    try:
+        #-- Create and submit request.
+        request=urllib2.Request(posixpath.join(*HOST))
+        response=urllib2.urlopen(request,timeout=timeout,context=context)
+    except (urllib2.HTTPError, urllib2.URLError):
+        raise Exception('List error from {0}'.format(posixpath.join(*HOST)))
+    else:
+        #-- read and parse request for files (column names and modified times)
+        tree = lxml.etree.parse(response,parser)
+        colnames = tree.xpath('//tr/td[not(@*)]//a/@href')
+        #-- get the Unix timestamp value for a modification time
+        lastmod = [get_unix_time(i,format=format)
+            for i in tree.xpath('//tr/td[@align="right"][1]/text()')]
+        #-- reduce using regular expression pattern
+        if pattern:
+            i = [i for i,f in enumerate(colnames) if re.search(pattern,f)]
+            #-- reduce list of column names and last modified times
+            colnames = [colnames[indice] for indice in i]
+            lastmod = [lastmod[indice] for indice in i]
+        #-- sort the list
+        if sort:
+            i = [i for i,j in sorted(enumerate(colnames), key=lambda i: i[1])]
+            #-- sort list of column names and last modified times
+            colnames = [colnames[indice] for indice in i]
+            lastmod = [lastmod[indice] for indice in i]
+        #-- return the list of column names and last modified times
+        return (colnames,lastmod)
+
 #-- PURPOSE: download a file from a http host
 def from_http(HOST,timeout=None,context=ssl.SSLContext(),local=None,hash='',
     chunk=16384,verbose=False,fid=sys.stdout,mode=0o775):

requirements.txt

@@ -1,14 +1,12 @@
-setuptools_scm
-numpy
-scipy
-pyproj
-python-dateutil
-matplotlib
 cartopy --no-binary=cartopy
-netCDF4
-h5py
 gdal
+h5py
 lxml
+matplotlib
+netCDF4
+numpy
+pyproj
+python-dateutil
 pyyaml
-ipywidgets
-ipyleaflet
+scipy
+setuptools_scm

test/test_spatial.py

@@ -282,3 +282,14 @@ def test_convert_ellipsoid():
     maxelevdel = 100.0*np.max(h_WGS84-elev_WGS84)
     assert np.isclose([minlatdel,maxlatdel],[-2.1316282e-14,2.1316282e-14]).all()
     assert np.isclose([minelevdel,maxelevdel],[-1.3287718e-7,1.6830199e-7]).all()
+
+#-- PURPOSE: test wrap longitudes
+def test_wrap_longitudes():
+    #-- number of data points
+    lon = np.arange(360)
+    obs = pyTMD.spatial.wrap_longitudes(lon)
+    #-- expected wrapped longitudes
+    exp = np.zeros((360))
+    exp[:181] = np.arange(181)
+    exp[181:] = np.arange(-179,0)
+    assert np.isclose(obs,exp).all()