Coverage for pygeodesy/elevations.py: 80%
70 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-12 11:45 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-12 11:45 -0400
2# -*- coding: utf-8 -*-
4u'''Web-services-based elevations and geoid heights.
6Functions to obtain elevations and geoid heights thru web services,
7for (lat, lon) locations, currently limited to the U{Conterminous
8US (CONUS)<https://WikiPedia.org/wiki/Contiguous_United_States>},
9see also modules L{pygeodesy.geoids} and L{pygeodesy.heights} and
10U{USGS10mElev.py<https://Gist.GitHub.com/pyRobShrk>}.
12B{macOS}: If an C{SSLCertVerificationError} occurs, especially
13this I{"[SSL: CERTIFICATE_VERIFY_FAILED] certificate verify
14failed: self "signed certificate in certificate chain ..."},
15review U{this post<https://StackOverflow.com/questions/27835619/
16urllib-and-ssl-certificate-verify-failed-error>} for a remedy.
17From a C{Terminal} window run:
18C{"/Applications/Python\\ X.Y/Install\\ Certificates.command"}
19'''
21from pygeodesy.basics import clips, ub2str
22from pygeodesy.errors import ParseError, _xkwds_get
23from pygeodesy.interns import NN, _AMPERSAND_, _COLONSPACE_, \
24 _elevation_, _height_, _LCURLY_, \
25 _n_a_, _no_, _RCURLY_, _SPACE_
26from pygeodesy.lazily import _ALL_LAZY
27from pygeodesy.named import _NamedTuple
28from pygeodesy.streprs import fabs, Fmt, fstr, lrstrip
29from pygeodesy.units import Lat, Lon, Meter, Scalar, Str
31# from math import fabs # from .karney
33__all__ = _ALL_LAZY.elevations
34__version__ = '23.03.19'
36try:
37 from urllib2 import urlopen # quote, urlcleanup
38 from httplib import HTTPException as HTTPError
40except (ImportError, NameError): # Python 3+
41 from urllib.request import urlopen # urlcleanup
42 # from urllib.parse import quote
43 from urllib.error import HTTPError
45_JSON_ = 'JSON'
46_QUESTION_ = '?'
47_XML_ = 'XML'
49try:
50 from json import loads as _json
51except ImportError:
53 from pygeodesy.interns import _COMMA_, _QUOTE2_
54 _QUOTE2COLONSPACE_ = _QUOTE2_ + _COLONSPACE_
56 def _json(ngs):
57 '''(INTERNAL) Convert an NGS response in JSON to a C{dict}.
58 '''
59 # b'{"geoidModel": "GEOID12A",
60 # "station": "UserStation",
61 # "lat": 37.8816,
62 # "latDms": "N375253.76000",
63 # "lon": -121.9142,
64 # "lonDms": "W1215451.12000",
65 # "geoidHeight": -31.703,
66 # "error": 0.064
67 # }'
68 #
69 # or in case of errors:
70 #
71 # b'{"error": "No suitable Geoid model found for model 15"
72 # }'
73 d = {}
74 for t in lrstrip(ngs.strip(), lrpairs={_LCURLY_: _RCURLY_}).split(_COMMA_):
75 t = t.strip()
76 j = t.strip(_QUOTE2_).split(_QUOTE2COLONSPACE_)
77 if len(j) != 2:
78 raise ParseError(json=t)
79 k, v = j
80 try:
81 v = float(v)
82 except (TypeError, ValueError):
83 v = Str(ub2str(v.lstrip().lstrip(_QUOTE2_)), name=k)
84 d[k] = v
85 return d
88def _error(fun, lat, lon, e):
89 '''(INTERNAL) Format an error
90 '''
91 return _COLONSPACE_(Fmt.PAREN(fun.__name__, fstr((lat, lon))), e)
94def _qURL(url, timeout=2, **params):
95 '''(INTERNAL) Build B{C{url}} query, get and verify response.
96 '''
97 if params: # build url query, don't map(quote, params)!
98 p = _AMPERSAND_(*(Fmt.EQUAL(p, v) for p, v in params.items() if v))
99 if p:
100 url = NN(url, _QUESTION_, p)
101 u = urlopen(url, timeout=timeout) # secs
103 s = u.getcode()
104 if s != 200: # http.HTTPStatus.OK or http.client.OK
105 raise HTTPError('code %d: %s' % (s, u.geturl()))
107 r = u.read()
108 u.close()
109 # urlcleanup()
110 return ub2str(r).strip()
113def _xml(tag, xml):
114 '''(INTERNAL) Get a <tag>value</tag> from XML.
115 '''
116 # b'<?xml version="1.0" encoding="utf-8" ?>
117 # <USGS_Elevation_Point_Query_Service>
118 # <Elevation_Query x="-121.914200" y="37.881600">
119 # <Data_Source>3DEP 1/3 arc-second</Data_Source>
120 # <Elevation>3851.03</Elevation>
121 # <Units>Feet</Units>
122 # </Elevation_Query>
123 # </USGS_Elevation_Point_Query_Service>'
124 i = xml.find(Fmt.TAG(tag))
125 if i > 0:
126 i += len(tag) + 2
127 j = xml.find(Fmt.TAGEND(tag), i)
128 if j > i:
129 return Str(xml[i:j].strip(), name=tag)
130 return _no_(_XML_, Fmt.TAG(tag)) # PYCHOK no cover
133class Elevation2Tuple(_NamedTuple): # .elevations.py
134 '''2-Tuple C{(elevation, data_source)} in C{meter} and C{str}.
135 '''
136 _Names_ = (_elevation_, 'data_source')
137 _Units_ = ( Meter, Str)
140def elevation2(lat, lon, timeout=2.0):
141 '''Get the geoid elevation at an C{NAD83} to C{NAVD88} location.
143 @arg lat: Latitude (C{degrees}).
144 @arg lon: Longitude (C{degrees}).
145 @kwarg timeout: Optional, query timeout (seconds).
147 @return: An L{Elevation2Tuple}C{(elevation, data_source)}
148 or (C{None, "error"}) in case of errors.
150 @raise ValueError: Invalid B{C{timeout}}.
152 @note: The returned C{elevation} is C{None} if B{C{lat}} or B{C{lon}} is
153 invalid or outside the C{Conterminous US (CONUS)}, if conversion
154 failed or if the query timed out. The C{"error"} is the C{HTTP-,
155 IO-, SSL-} or other C{-Error} as a string (C{str}).
157 @see: U{USGS Elevation Point Query Service<https://NationalMap.gov/epqs>}, the
158 U{FAQ<https://www.USGS.gov/faqs/what-are-projection-horizontal-and-vertical-
159 datum-units-and-resolution-3dep-standard-dems>}, U{geoid.py<https://Gist.GitHub.com/
160 pyRobShrk>}, module L{geoids}, classes L{GeoidG2012B}, L{GeoidKarney} and
161 L{GeoidPGM}.
162 '''
163 try: # alt 'https://NED.USGS.gov/epqs/pqs.php'
164 x = _qURL('https://NationalMap.USGS.gov/epqs/pqs.php',
165 x=Lon(lon).toStr(prec=6),
166 y=Lat(lat).toStr(prec=6),
167 units='Meters', # 'Feet', capitalized
168 output=_XML_.lower(), # _JSON_, lowercase only
169 timeout=Scalar(timeout=timeout))
170 if x[:6] == '<?xml ':
171 e = _xml('Elevation', x)
172 try:
173 e = float(e)
174 if fabs(e) < 1e6:
175 return Elevation2Tuple(e, _xml('Data_Source', x))
176 e = 'non-CONUS %.2F' % (e,)
177 except (TypeError, ValueError):
178 pass
179 else: # PYCHOK no cover
180 e = _no_(_XML_, Fmt.QUOTE2(clips(x, limit=128, white=_SPACE_)))
181 except Exception as x: # (HTTPError, IOError, TypeError, ValueError)
182 e = repr(x)
183 e = _error(elevation2, lat, lon, e)
184 return Elevation2Tuple(None, e)
187class GeoidHeight2Tuple(_NamedTuple): # .elevations.py
188 '''2-Tuple C{(height, model_name)}, geoid C{height} in C{meter}
189 and C{model_name} as C{str}.
190 '''
191 _Names_ = (_height_, 'model_name')
192 _Units_ = ( Meter, Str)
195def geoidHeight2(lat, lon, model=0, timeout=2.0):
196 '''Get the C{NAVD88} geoid height at an C{NAD83} location.
198 @arg lat: Latitude (C{degrees}).
199 @arg lon: Longitude (C{degrees}).
200 @kwarg model: Optional, geoid model ID (C{int}).
201 @kwarg timeout: Optional, query timeout (seconds).
203 @return: An L{GeoidHeight2Tuple}C{(height, model_name)}
204 or C{(None, "error"}) in case of errors.
206 @raise ValueError: Invalid B{C{timeout}}.
208 @note: The returned C{height} is C{None} if B{C{lat}} or B{C{lon}} is
209 invalid or outside the C{Conterminous US (CONUS)}, if the
210 B{C{model}} was invalid, if conversion failed or if the query
211 timed out. The C{"error"} is the C{HTTP-, IO-, SSL-, URL-} or
212 other C{-Error} as a string (C{str}).
214 @see: U{NOAA National Geodetic Survey
215 <https://www.NGS.NOAA.gov/INFO/geodesy.shtml>},
216 U{Geoid<https://www.NGS.NOAA.gov/web_services/geoid.shtml>},
217 U{USGS10mElev.py<https://Gist.GitHub.com/pyRobShrk>}, module
218 L{geoids}, classes L{GeoidG2012B}, L{GeoidKarney} and
219 L{GeoidPGM}.
220 '''
221 try:
222 j = _qURL('https://Geodesy.NOAA.gov/api/geoid/ght',
223 lat=Lat(lat).toStr(prec=6),
224 lon=Lon(lon).toStr(prec=6),
225 model=(model if model else NN),
226 timeout=Scalar(timeout=timeout)) # PYCHOK indent
227 if j[:1] == _LCURLY_ and j[-1:] == _RCURLY_ and j.find('"error":') > 0:
228 d, e = _json(j), 'geoidHeight'
229 if isinstance(_xkwds_get(d, error=_n_a_), float):
230 h = d.get(e, None)
231 if h is not None:
232 m = _xkwds_get(d, geoidModel=_n_a_)
233 return GeoidHeight2Tuple(h, m)
234 else:
235 e = _JSON_
236 e = _no_(e, Fmt.QUOTE2(clips(j, limit=256, white=_SPACE_)))
237 except Exception as x: # (HTTPError, IOError, ParseError, TypeError, ValueError)
238 e = repr(x)
239 e = _error(geoidHeight2, lat, lon, e)
240 return GeoidHeight2Tuple(None, e)
243if __name__ == '__main__':
245 # <https://WikiPedia.org/wiki/Mount_Diablo>
246 for f in (elevation2, # (1173.79, '3DEP 1/3 arc-second')
247 geoidHeight2): # (-31.699, u'GEOID12B')
248 t = f(37.8816, -121.9142)
249 print(_COLONSPACE_(f.__name__, t))
251# **) MIT License
252#
253# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
254#
255# Permission is hereby granted, free of charge, to any person obtaining a
256# copy of this software and associated documentation files (the "Software"),
257# to deal in the Software without restriction, including without limitation
258# the rights to use, copy, modify, merge, publish, distribute, sublicense,
259# and/or sell copies of the Software, and to permit persons to whom the
260# Software is furnished to do so, subject to the following conditions:
261#
262# The above copyright notice and this permission notice shall be included
263# in all copies or substantial portions of the Software.
264#
265# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
266# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
267# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
268# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
269# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
270# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
271# OTHER DEALINGS IN THE SOFTWARE.
273# % python -m pygeodesy.elevations
274# elevation2: (1173.79, '3DEP 1/3 arc-second')
275# geoidHeight2: (-31.703, 'GEOID12B')