Coverage for pygeodesy/elevations.py: 80%

70 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-04-12 11:45 -0400

1 

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

3 

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

5 

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

11 

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

20 

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 

30 

31# from math import fabs # from .karney 

32 

33__all__ = _ALL_LAZY.elevations 

34__version__ = '23.03.19' 

35 

36try: 

37 from urllib2 import urlopen # quote, urlcleanup 

38 from httplib import HTTPException as HTTPError 

39 

40except (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 

49try: 

50 from json import loads as _json 

51except 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 

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

89 '''(INTERNAL) Format an error 

90 ''' 

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

92 

93 

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 

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 

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 

131 

132 

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) 

138 

139 

140def 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 

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) 

193 

194 

195def 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 

243if __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')