1

2# -*- coding: utf-8 -*-

3

4 u'''Web-services-based elevations and geoid heights.

5

6 Functions to obtain elevations and geoid heights thru web services,

7 for (lat, lon) locations, currently limited to the U{Conterminous

8 US (CONUS)<https://WikiPedia.org/wiki/Contiguous_United_States>},

9 see also modules L{pygeodesy.geoids} and L{pygeodesy.heights} and

10 U{USGS10mElev.py<https://Gist.GitHub.com/pyRobShrk>}.

11

12 B{macOS}: If an C{SSLCertVerificationError} occurs, especially

13 this I{"[SSL: CERTIFICATE_VERIFY_FAILED] certificate verify

14 failed: self "signed certificate in certificate chain ..."},

15 review U{this post<https://StackOverflow.com/questions/27835619/

16 urllib-and-ssl-certificate-verify-failed-error>} for a remedy.

17 From a C{Terminal} window run:

18 C{"/Applications/Python\\ X.Y/Install\\ Certificates.command"}

19'''

20

21 from pygeodesy.basics import clips, ub2str

22 from pygeodesy.errors import ParseError, _xkwds_get

23 from pygeodesy.interns import NN, _AMPERSAND_, _COLONSPACE_, \

24 _elevation_, _height_, _LCURLY_, \

25 _n_a_, _no_, _RCURLY_, _SPACE_

26 from pygeodesy.lazily import _ALL_LAZY

27 from pygeodesy.named import _NamedTuple

28 from pygeodesy.streprs import fabs, Fmt, fstr, lrstrip

29 from pygeodesy.units import Lat, Lon, Meter, Scalar, Str

30

31# from math import fabs # from .karney

32

33__all__ = _ALL_LAZY.elevations

34__version__ = '23.03.19'

35

36 try:

37 from urllib2 import urlopen # quote, urlcleanup

38 from httplib import HTTPException as HTTPError

39

40 except (ImportError, NameError): # Python 3+

41 from urllib.request import urlopen # urlcleanup

42 # from urllib.parse import quote

43 from urllib.error import HTTPError

44

45_JSON_ = 'JSON'

46_QUESTION_ = '?'

47_XML_ = 'XML'

48

49 try:

50 from json import loads as _json

51 except ImportError:

52

53 from pygeodesy.interns import _COMMA_, _QUOTE2_

54 _QUOTE2COLONSPACE_ = _QUOTE2_ + _COLONSPACE_

55

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

86

87

88 def _error(fun, lat, lon, e):

89 '''(INTERNAL) Format an error

90 '''

91 return _COLONSPACE_(Fmt.PAREN(fun.__name__, fstr((lat, lon))), e)

92

93

94 def _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

102

103 s = u.getcode()

104 if s != 200: # http.HTTPStatus.OK or http.client.OK

105 raise HTTPError('code %d: %s' % (s, u.geturl()))

106

107 r = u.read()

108 u.close()

109 # urlcleanup()

110 return ub2str(r).strip()

111

112

113 def _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

131

132

133 class 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)

138

139

140 def elevation2(lat, lon, timeout=2.0):

141 '''Get the geoid elevation at an C{NAD83} to C{NAVD88} location.

142

143 @arg lat: Latitude (C{degrees}).

144 @arg lon: Longitude (C{degrees}).

145 @kwarg timeout: Optional, query timeout (seconds).

146

147 @return: An L{Elevation2Tuple}C{(elevation, data_source)}

148 or (C{None, "error"}) in case of errors.

149

150 @raise ValueError: Invalid B{C{timeout}}.

151

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}).

156

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)

185

186

187 class 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)

193

194

195 def geoidHeight2(lat, lon, model=0, timeout=2.0):

196 '''Get the C{NAVD88} geoid height at an C{NAD83} location.

197

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).

202

203 @return: An L{GeoidHeight2Tuple}C{(height, model_name)}

204 or C{(None, "error"}) in case of errors.

205

206 @raise ValueError: Invalid B{C{timeout}}.

207

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}).

213

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)

241

242

243 if __name__ == '__main__':

244

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

250

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.

272

273# % python -m pygeodesy.elevations

274# elevation2: (1173.79, '3DEP 1/3 arc-second')

275# geoidHeight2: (-31.703, 'GEOID12B')