From 9465aae56d7ba6724f769f237978c61c7f88bf4c Mon Sep 17 00:00:00 2001 From: Tyler Sutterley Date: Thu, 7 Oct 2021 08:59:53 -0700 Subject: [PATCH] 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 | 1 + doc/source/user_guide/spatial.rst | 10 ++++- doc/source/user_guide/utilities.rst | 31 ++++++++++++- environment.yml | 24 +++++----- notebooks/Check Tide Map.ipynb | 9 ++-- notebooks/Plot Tide Forecasts.ipynb | 9 ++-- pyTMD/spatial.py | 31 +++++++++---- pyTMD/utilities.py | 57 +++++++++++++++++++++++- requirements.txt | 18 ++++---- test/test_spatial.py | 11 +++++ test/test_time.py | 6 +-- 11 files changed, 162 insertions(+), 45 deletions(-) diff --git a/doc/source/getting_started/Citations.rst b/doc/source/getting_started/Citations.rst index 6becdc81..5cf8ad3c 100644 --- a/doc/source/getting_started/Citations.rst +++ b/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 `_ - `ipywidgets: interactive HTML widgets for Jupyter notebooks and IPython `_ - `ipyleaflet: Jupyter / Leaflet bridge enabling interactive maps `_ +- `setuptools_scm: manager for python package versions using scm metadata `_ Credits ####### diff --git a/doc/source/user_guide/spatial.rst b/doc/source/user_guide/spatial.rst index 181480df..42567323 100644 --- a/doc/source/user_guide/spatial.rst +++ b/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 diff --git a/doc/source/user_guide/utilities.rst b/doc/source/user_guide/utilities.rst index f3abb710..9fd18d27 100644 --- a/doc/source/user_guide/utilities.rst +++ b/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 diff --git a/environment.yml b/environment.yml index 7e1b034b..c47851c3 100644 --- a/environment.yml +++ b/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 diff --git a/notebooks/Check Tide Map.ipynb b/notebooks/Check Tide Map.ipynb index 4e8a5e78..7faff83d 100644 --- a/notebooks/Check Tide Map.ipynb +++ b/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", diff --git a/notebooks/Plot Tide Forecasts.ipynb b/notebooks/Plot Tide Forecasts.ipynb index ff523dbe..197db6d4 100644 --- a/notebooks/Plot Tide Forecasts.ipynb +++ b/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", diff --git a/pyTMD/spatial.py b/pyTMD/spatial.py index 1d5978a7..b1ef5c03 100644 --- a/pyTMD/spatial.py +++ b/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 diff --git a/pyTMD/utilities.py b/pyTMD/utilities.py index a1796983..7810c73e 100644 --- a/pyTMD/utilities.py +++ b/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): diff --git a/requirements.txt b/requirements.txt index 213616f8..b2c2de4b 100644 --- a/requirements.txt +++ b/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 diff --git a/test/test_spatial.py b/test/test_spatial.py index bf51395d..601b3ca7 100644 --- a/test/test_spatial.py +++ b/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() diff --git a/test/test_time.py b/test/test_time.py index 88b85604..15d1b265 100644 --- a/test/test_time.py +++ b/test/test_time.py @@ -33,13 +33,12 @@ def test_julian(YEAR,MONTH): YY,MM,DD,HH,MN,SS = pyTMD.time.convert_julian(JD, FORMAT='tuple', ASTYPE=np.float64) #-- assert dates - eps = np.finfo(np.float16).eps assert (YY == YEAR) assert (MM == MONTH) assert (DD == DAY) assert (HH == HOUR) assert (MN == MINUTE) - assert (np.abs(SS - SECOND) < eps) + assert np.isclose(SS, SECOND) #-- parameterize calendar dates @pytest.mark.parametrize("YEAR", np.random.randint(1992,2020,size=2)) @@ -76,13 +75,12 @@ def test_decimal_dates(YEAR,MONTH): minute_temp = np.mod(hour_temp,1)*60.0 second = np.mod(minute_temp,1)*60.0 #-- assert dates - eps = np.finfo(np.float16).eps assert (np.floor(tdec) == YEAR) assert (month == MONTH) assert (day == DAY) assert (np.floor(hour_temp) == HOUR) assert (np.floor(minute_temp) == MINUTE) - assert (np.abs(second - SECOND) < eps) + assert np.isclose(second, SECOND) #-- PURPOSE: test UNIX time def test_unix_time():