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():