1

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

3

4 u'''Web-services-based elevations and C{CONUS} 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@see: Module L{pygeodesy.geoids} to get geoid heights from other

13 sources and for regions other than C{CONUS}.

14

15@note: If on B{macOS} an C{SSLCertVerificationError} occurs, like

16 I{"[SSL: CERTIFICATE_VERIFY_FAILED] certificate verify failed:

17 self "signed certificate in certificate chain ..."}, review

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

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

20 From a C{Terminal} window run:

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

22'''

23

24 from pygeodesy.basics import clips, ub2str

25 from pygeodesy.errors import ParseError, _xkwds_get

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

27 _elevation_, _height_, _LCURLY_, \

28 _n_a_, _no_, _RCURLY_, _SPACE_

29 from pygeodesy.lazily import _ALL_LAZY

30 from pygeodesy.named import _NamedTuple

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

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

33

34# from math import fabs # from .karney

35

36__all__ = _ALL_LAZY.elevations

37__version__ = '23.05.04'

38

39 try:

40 from urllib2 import urlopen # quote, urlcleanup

41 from httplib import HTTPException as HTTPError

42

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

44 from urllib.request import urlopen # urlcleanup

45 # from urllib.parse import quote

46 from urllib.error import HTTPError

47

48_JSON_ = 'JSON'

49_QUESTION_ = '?'

50_XML_ = 'XML'

51

52 try:

53 from json import loads as _json

54 except ImportError:

55

56 from pygeodesy.interns import _COMMA_, _QUOTE2_

57 _QUOTE2COLONSPACE_ = _QUOTE2_ + _COLONSPACE_

58

59 def _json(ngs):

60 '''(INTERNAL) Convert an NGS response in JSON to a C{dict}.

61 '''

62 # b'{"geoidModel": "GEOID12A",

63 # "station": "UserStation",

64 # "lat": 37.8816,

65 # "latDms": "N375253.76000",

66 # "lon": -121.9142,

67 # "lonDms": "W1215451.12000",

68 # "geoidHeight": -31.703,

69 # "error": 0.064

70 # }'

71 #

72 # or in case of errors:

73 #

74 # b'{"error": "No suitable Geoid model found for model 15"

75 # }'

76 d = {}

77 for t in lrstrip(ngs.strip(), lrpairs={_LCURLY_: _RCURLY_}).split(_COMMA_):

78 t = t.strip()

79 j = t.strip(_QUOTE2_).split(_QUOTE2COLONSPACE_)

80 if len(j) != 2:

81 raise ParseError(json=t)

82 k, v = j

83 try:

84 v = float(v)

85 except (TypeError, ValueError):

86 v = Str(ub2str(v.lstrip().lstrip(_QUOTE2_)), name=k)

87 d[k] = v

88 return d

89

90

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

92 '''(INTERNAL) Format an error

93 '''

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

95

96

97 def _qURL(url, timeout=2, **params):

98 '''(INTERNAL) Build B{C{url}} query, get and verify response.

99 '''

100 if params: # build url query, don't map(quote, params)!

101 p = _AMPERSAND_(*(Fmt.EQUAL(p, v) for p, v in params.items() if v))

102 if p:

103 url = NN(url, _QUESTION_, p)

104 u = urlopen(url, timeout=timeout) # secs

105

106 s = u.getcode()

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

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

109

110 r = u.read()

111 u.close()

112 # urlcleanup()

113 return ub2str(r).strip()

114

115

116 def _xml(tag, xml):

117 '''(INTERNAL) Get a <tag>value</tag> from XML.

118 '''

119 # b'<?xml version="1.0" encoding="utf-8" ?>

120 # <USGS_Elevation_Point_Query_Service>

121 # <Elevation_Query x="-121.914200" y="37.881600">

122 # <Data_Source>3DEP 1/3 arc-second</Data_Source>

123 # <Elevation>3851.03</Elevation>

124 # <Units>Feet</Units>

125 # </Elevation_Query>

126 # </USGS_Elevation_Point_Query_Service>'

127 i = xml.find(Fmt.TAG(tag))

128 if i > 0:

129 i += len(tag) + 2

130 j = xml.find(Fmt.TAGEND(tag), i)

131 if j > i:

132 return Str(xml[i:j].strip(), name=tag)

133 return _no_(_XML_, Fmt.TAG(tag)) # PYCHOK no cover

134

135

136 class Elevation2Tuple(_NamedTuple): # .elevations.py

137 '''2-Tuple C{(elevation, data_source)} in C{meter} and C{str}.

138 '''

139 _Names_ = (_elevation_, 'data_source')

140 _Units_ = ( Meter, Str)

141

142

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

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

145

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

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

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

149

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

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

152

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

154

155 @note: The returned C{elevation} is C{None} if B{C{lat}} or B{C{lon}} is

156 invalid or outside the C{Conterminous US (CONUS)}, if conversion

157 failed or if the query timed out. The C{"error"} is the C{HTTP-,

158 IO-, SSL-} or other C{-Error} as a string (C{str}).

159

160 @see: U{USGS Elevation Point Query Service<https://apps.NationalMap.gov/epqs/>}, the

161 U{FAQ<https://www.USGS.gov/faqs/what-are-projection-horizontal-and-vertical-

162 datum-units-and-resolution-3dep-standard-dems>}, U{geoid.py<https://Gist.GitHub.com/

163 pyRobShrk>}, module L{geoids}, classes L{GeoidG2012B}, L{GeoidKarney} and

164 L{GeoidPGM}.

165 '''

166 try: # alt 'https://NED.USGS.gov/epqs/pqs.php', 'https://epqs.NationalMap.gov/v1'

167 x = _qURL('https://NationalMap.USGS.gov/epqs/pqs.php',

168 x=Lon(lon).toStr(prec=6),

169 y=Lat(lat).toStr(prec=6),

170 units='Meters', # 'Feet', capitalized

171 output=_XML_.lower(), # _JSON_, lowercase only

172 timeout=Scalar(timeout=timeout))

173 if x[:6] == '<?xml ':

174 e = _xml('Elevation', x)

175 try:

176 e = float(e)

177 if fabs(e) < 1e6:

178 return Elevation2Tuple(e, _xml('Data_Source', x))

179 e = 'non-CONUS %.2F' % (e,)

180 except (TypeError, ValueError):

181 pass

182 else: # PYCHOK no cover

183 e = _no_(_XML_, Fmt.QUOTE2(clips(x, limit=128, white=_SPACE_)))

184 except Exception as x: # (HTTPError, IOError, TypeError, ValueError)

185 e = repr(x)

186 e = _error(elevation2, lat, lon, e)

187 return Elevation2Tuple(None, e)

188

189

190 class GeoidHeight2Tuple(_NamedTuple): # .elevations.py

191 '''2-Tuple C{(height, model_name)}, geoid C{height} in C{meter}

192 and C{model_name} as C{str}.

193 '''

194 _Names_ = (_height_, 'model_name')

195 _Units_ = ( Meter, Str)

196

197

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

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

200

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

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

203 @kwarg model: Optional, geoid model ID (C{int}).

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

205

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

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

208

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

210

211 @note: The returned C{height} is C{None} if B{C{lat}} or B{C{lon}}

212 is invalid or outside the C{Conterminous US (CONUS)}, if the

213 B{C{model}} was invalid, if conversion failed or if the query

214 timed out. The C{"error"} is the C{HTTP-, IO-, SSL-, URL-}

215 or other C{-Error} as a string (C{str}).

216

217 @see: U{NOAA National Geodetic Survey

218 <https://www.NGS.NOAA.gov/INFO/geodesy.shtml>},

219 U{Geoid<https://www.NGS.NOAA.gov/web_services/geoid.shtml>},

220 U{USGS10mElev.py<https://Gist.GitHub.com/pyRobShrk>}, module

221 L{geoids}, classes L{GeoidG2012B}, L{GeoidKarney} and

222 L{GeoidPGM}.

223 '''

224 try:

225 j = _qURL('https://Geodesy.NOAA.gov/api/geoid/ght',

226 lat=Lat(lat).toStr(prec=6),

227 lon=Lon(lon).toStr(prec=6),

228 model=(model if model else NN),

229 timeout=Scalar(timeout=timeout)) # PYCHOK indent

230 if j[:1] == _LCURLY_ and j[-1:] == _RCURLY_ and j.find('"error":') > 0:

231 d, e = _json(j), 'geoidHeight'

232 if isinstance(_xkwds_get(d, error=_n_a_), float):

233 h = d.get(e, None)

234 if h is not None:

235 m = _xkwds_get(d, geoidModel=_n_a_)

236 return GeoidHeight2Tuple(h, m)

237 else:

238 e = _JSON_

239 e = _no_(e, Fmt.QUOTE2(clips(j, limit=256, white=_SPACE_)))

240 except Exception as x: # (HTTPError, IOError, ParseError, TypeError, ValueError)

241 e = repr(x)

242 e = _error(geoidHeight2, lat, lon, e)

243 return GeoidHeight2Tuple(None, e)

244

245

246 if __name__ == '__main__':

247

248 # <https://WikiPedia.org/wiki/Mount_Diablo>

249 for f in (elevation2, # (1173.79, '3DEP 1/3 arc-second')

250 geoidHeight2): # (-31.699, u'GEOID12B')

251 t = f(37.8816, -121.9142)

252 print(_COLONSPACE_(f.__name__, t))

253

254# **) MIT License

255#

256# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.

257#

258# Permission is hereby granted, free of charge, to any person obtaining a

259# copy of this software and associated documentation files (the "Software"),

260# to deal in the Software without restriction, including without limitation

261# the rights to use, copy, modify, merge, publish, distribute, sublicense,

262# and/or sell copies of the Software, and to permit persons to whom the

263# Software is furnished to do so, subject to the following conditions:

264#

265# The above copyright notice and this permission notice shall be included

266# in all copies or substantial portions of the Software.

267#

268# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

269# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

270# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

271# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

272# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

273# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

274# OTHER DEALINGS IN THE SOFTWARE.

275

276# % python -m pygeodesy.elevations

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

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