Coverage for pygeodesy/osgr.py: 97%

303 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-12-02 13:46 -0500

1 

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

3 

4u'''Ordnance Survey Grid References (OSGR) references on the UK U{National Grid 

5<https://www.OrdnanceSurvey.co.UK/documents/resources/guide-to-nationalgrid.pdf>}. 

6 

7Classes L{Osgr} and L{OSGRError} and functions L{parseOSGR} and L{toOsgr}. 

8 

9A pure Python implementation, transcoded from I{Chris Veness}' JavaScript originals U{OS National Grid 

10<https://www.Movable-Type.co.UK/scripts/latlong-os-gridref.html>} and U{Module osgridref 

11<https://www.Movable-Type.co.UK/scripts/geodesy/docs/module-osgridref.html>} and I{Charles Karney}'s 

12C++ class U{OSGB<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1OSGB.html>}. 

13 

14OSGR provides geocoordinate references for UK mapping purposes, converted in 2015 to work with the C{WGS84} 

15or the original C{OSGB36} datum. In addition, this implementation includes both the OS recommended and the 

16Krüger-based method to convert between OSGR and geodetic coordinates (with keyword argument C{kTM} of 

17function L{toOsgr}, method L{Osgr.toLatLon} and method C{toOsgr} of any ellipsoidal C{LatLon} class). 

18 

19See U{Transverse Mercator: Redfearn series<https://WikiPedia.org/wiki/Transverse_Mercator:_Redfearn_series>}, 

20Karney's U{"Transverse Mercator with an accuracy of a few nanometers", 2011<https://ArXiv.org/pdf/1002.1417v3.pdf>} 

21(building on U{"Konforme Abbildung des Erdellipsoids in der Ebene", 1912<https://bib.GFZ-Potsdam.DE/pub/digi/krueger2.pdf>}, 

22U{"Die Mathematik der Gauß-Krueger-Abbildung", 2006<https://DE.WikiPedia.org/wiki/Gauß-Krüger-Koordinatensystem>}, 

23U{A Guide<https://www.OrdnanceSurvey.co.UK/documents/resources/guide-coordinate-systems-great-britain.pdf>} 

24and U{Ordnance Survey National Grid<https://WikiPedia.org/wiki/Ordnance_Survey_National_Grid>}. 

25''' 

26# make sure int/int division yields float quotient, see .basics 

27from __future__ import division as _; del _ # PYCHOK semicolon 

28 

29from pygeodesy.basics import halfs2, isbool, isfloat, map1, \ 

30 _splituple, _xsubclassof 

31from pygeodesy.constants import _1_0, _10_0, _N_2_0 # PYCHOK used! 

32from pygeodesy.datums import Datums, _ellipsoidal_datum, _WGS84 

33# from pygeodesy.dms import parseDMS2 # _MODS 

34from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase as _LLEB 

35from pygeodesy.errors import _parseX, _TypeError, _ValueError, \ 

36 _xkwds, _xkwds_get 

37from pygeodesy.fmath import Fdot, fpowers, Fsum 

38# from pygeodesy.fsums import Fsum # from .fmath 

39from pygeodesy.interns import MISSING, NN, _A_, _COLON_, _COMMA_, \ 

40 _COMMASPACE_, _DOT_, _ellipsoidal_, \ 

41 _latlon_, _not_, _SPACE_ 

42from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS 

43from pygeodesy.named import _NamedBase, nameof, _xnamed 

44from pygeodesy.namedTuples import EasNor2Tuple, LatLon2Tuple, \ 

45 LatLonDatum3Tuple 

46from pygeodesy.props import Property_RO, property_RO 

47from pygeodesy.streprs import _EN_WIDE, enstr2, _enstr2m3, Fmt, \ 

48 _resolution10, unstr, _xzipairs 

49from pygeodesy.units import Easting, Lam_, Lat, Lon, Northing, \ 

50 Phi_, Scalar, _10um, _100km 

51from pygeodesy.utily import degrees90, degrees180, sincostan3, truncate 

52 

53from math import cos, fabs, radians, sin, sqrt 

54 

55__all__ = _ALL_LAZY.osgr 

56__version__ = '23.08.24' 

57 

58_equivalent_ = 'equivalent' 

59_OSGR_ = 'OSGR' 

60_ord_A = ord(_A_) 

61_TRIPS = 33 # .toLatLon convergence 

62 

63 

64class _NG(object): 

65 '''Ordnance Survey National Grid parameters. 

66 ''' 

67 @Property_RO 

68 def a0(self): # equatoradius, scaled 

69 return self.ellipsoid.a * self.k0 

70 

71 @Property_RO 

72 def b0(self): # polaradius, scaled 

73 return self.ellipsoid.b * self.k0 

74 

75 @Property_RO 

76 def datum(self): # datum, Airy130 ellipsoid 

77 return Datums.OSGB36 

78 

79 @Property_RO 

80 def eas0(self): # False origin easting (C{meter}) 

81 return Easting(4 * _100km) 

82 

83 @Property_RO 

84 def easX(self): # easting [0..extent] (C{meter}) 

85 return Easting(7 * _100km) 

86 

87 @Property_RO 

88 def ellipsoid(self): # ellipsoid, Airy130 

89 return self.datum.ellipsoid 

90 

91 def forward2(self, latlon): # convert C{latlon} to (easting, norting), as I{Karney}'s 

92 # U{Forward<https://GeographicLib.SourceForge.io/C++/doc/OSGB_8hpp_source.html>} 

93 t = self.kTM.forward(latlon.lat, latlon.lon, lon0=self.lon0) 

94 e = t.easting + self.eas0 

95 n = t.northing + self.nor0ffset 

96 return e, n 

97 

98 @Property_RO 

99 def k0(self): # central scale (C{float}), like I{Karney}'s CentralScale 

100 # <https://GeographicLib.SourceForge.io/C++/doc/OSGB_8hpp_source.html> 

101 _0_9998268 = (9998268 - 10000000) / 10000000 

102 return Scalar(_10_0**_0_9998268) # 0.9996012717... 

103 

104 @Property_RO 

105 def kTM(self): # the L{KTransverseMercator} instance, like I{Karney}'s OSGBTM 

106 # <https://GeographicLib.SourceForge.io/C++/doc/OSGB_8cpp_source.html> 

107 return _MODS.ktm.KTransverseMercator(self.datum, lon0=0, k0=self.k0) 

108 

109 @Property_RO 

110 def lam0(self): # True origin longitude C{radians} 

111 return Lam_(self.lon0) 

112 

113 @Property_RO 

114 def lat0(self): # True origin latitude, 49°N 

115 return Lat(49.0) 

116 

117 @Property_RO 

118 def lon0(self): # True origin longitude, 2°W 

119 return Lon(_N_2_0) 

120 

121 @Property_RO 

122 def Mabcd(self): # meridional coefficients (a, b, c, d) 

123 n, n2, n3 = fpowers(self.ellipsoid.n, 3) 

124 M = (Fsum(4, 4 * n, 5 * n2, 5 * n3) / 4, 

125 Fsum( 24 * n, 24 * n2, 21 * n3) / 8, 

126 Fsum( 15 * n2, 15 * n3) / 8, 

127 (35 * n3 / 24)) 

128 return M 

129 

130 def Mabcd0(self, a): # meridional arc, scaled 

131 c = a + self.phi0 

132 s = a - self.phi0 

133 R = Fdot(self.Mabcd, s, -sin(s) * cos(c), 

134 sin(s * 2) * cos(c * 2), 

135 -sin(s * 3) * cos(c * 3)) 

136 return float(R * self.b0) 

137 

138 @Property_RO 

139 def nor0(self): # False origin northing (C{meter}) 

140 return Northing(-_100km) 

141 

142 @Property_RO 

143 def nor0ffset(self): # like I{Karney}'s computenorthoffset 

144 # <https://GeographicLib.SourceForge.io/C++/doc/OSGB_8cpp_source.html> 

145 return self.nor0 - self.kTM.forward(self.lat0, 0).northing 

146 

147 @Property_RO 

148 def norX(self): # northing [0..extent] (C{meter}) 

149 return Northing(13 * _100km) 

150 

151 def nu_rho_eta3(self, sa): # 3-tuple (nu, nu / rho, eta2) 

152 E = self.ellipsoid # rho, nu = E.roc2_(sa) # .k0? 

153 s = E.e2s2(sa) # == 1 - E.e2 * sa**2 

154 v = self.a0 / sqrt(s) # == nu, transverse roc 

155 # rho = .a0 * E.e21 / s**1.5 == v * E.e21 / s 

156 # r = v * E.e21 / s # == rho, meridional roc 

157 # nu / rho == v / (v * E.e21 / s) == s / E.e21 == ... 

158 s *= E._1_e21 # ... s * E._1_e21 == s * E.a2_b2 

159 return v, s, (s - _1_0) # η2 = nu / rho - 1 

160 

161 @Property_RO 

162 def phi0(self): # True origin latitude C{radians} 

163 return Phi_(self.lat0) 

164 

165 def reverse(self, osgr): # convert C{osgr} to (ellipsoidal} LatLon, as I{Karney}'s 

166 # U{Reverse<https://GeographicLib.SourceForge.io/C++/doc/OSGB_8hpp_source.html>} 

167 r = osgr._latlonTM 

168 if r is None: 

169 x = osgr.easting - self.eas0 

170 y = osgr.northing - self.nor0ffset 

171 t = self.kTM.reverse(x, y, lon0=self.lon0) 

172 r = _LLEB(t.lat, t.lon, datum=self.datum, name=osgr.name) 

173 osgr._latlonTM = r 

174 return r 

175 

176_NG = _NG() # PYCHOK singleton 

177 

178 

179class OSGRError(_ValueError): 

180 '''Error raised for a L{parseOSGR}, L{Osgr} or other OSGR issue. 

181 ''' 

182 pass 

183 

184 

185class Osgr(_NamedBase): 

186 '''Ordnance Survey Grid References (OSGR) coordinates on 

187 the U{National Grid<https://www.OrdnanceSurvey.co.UK/ 

188 documents/resources/guide-to-nationalgrid.pdf>}. 

189 ''' 

190 _datum = _NG.datum # default datum (L{Datums.OSGB36}) 

191 _easting = 0 # Easting (C{meter}) 

192 _latlon = None # cached B{C{_toLatlon}} 

193 _latlonTM = None # cached B{C{_toLatlon kTM}} 

194 _northing = 0 # Nothing (C{meter}) 

195 _resolution = 0 # from L{parseOSGR} (C{meter}) 

196 

197 def __init__(self, easting, northing, datum=None, name=NN, 

198 resolution=0): 

199 '''New L{Osgr} coordinate. 

200 

201 @arg easting: Easting from the OS C{National Grid} 

202 origin (C{meter}). 

203 @arg northing: Northing from the OS C{National Grid} 

204 origin (C{meter}). 

205 @kwarg datum: Override default datum (C{Datums.OSGB36}). 

206 @kwarg name: Optional name (C{str}). 

207 @kwarg resolution: Optional resolution (C{meter}), 

208 C{0} for default. 

209 

210 @raise OSGRError: Invalid or negative B{C{easting}} or 

211 B{C{northing}} or B{C{datum}} not an 

212 C{Datums.OSGB36} equivalent. 

213 ''' 

214 if datum: # PYCHOK no cover 

215 try: 

216 self._datum = _ellipsoidal_datum(datum) 

217 if self.datum != _NG.datum: 

218 raise ValueError(_not_(_NG.datum.name, _equivalent_)) 

219 except (TypeError, ValueError) as x: 

220 raise OSGRError(datum=datum, cause=x) 

221 

222 self._easting = Easting( easting, Error=OSGRError, high=_NG.easX) 

223 self._northing = Northing(northing, Error=OSGRError, high=_NG.norX) 

224 

225 if name: 

226 self.name = name 

227 if resolution: 

228 self._resolution = _resolution10(resolution, Error=OSGRError) 

229 

230 def __str__(self): 

231 return self.toStr(GD=True, sep=_SPACE_) 

232 

233 @Property_RO 

234 def datum(self): 

235 '''Get the datum (L{Datum}). 

236 ''' 

237 return self._datum 

238 

239 @Property_RO 

240 def easting(self): 

241 '''Get the easting (C{meter}). 

242 ''' 

243 return self._easting 

244 

245 @Property_RO 

246 def falsing0(self): 

247 '''Get the C{OS National Grid} falsing (L{EasNor2Tuple}). 

248 ''' 

249 return EasNor2Tuple(_NG.eas0, _NG.nor0, name=_OSGR_) 

250 

251 @property_RO 

252 def iteration(self): 

253 '''Get the most recent C{Osgr.toLatLon} iteration number 

254 (C{int}) or C{None} if not available/applicable. 

255 ''' 

256 return self._iteration 

257 

258 @Property_RO 

259 def latlon0(self): 

260 '''Get the C{OS National Grid} origin (L{LatLon2Tuple}). 

261 ''' 

262 return LatLon2Tuple(_NG.lat, _NG.lon0, name=_OSGR_) 

263 

264 @Property_RO 

265 def northing(self): 

266 '''Get the northing (C{meter}). 

267 ''' 

268 return self._northing 

269 

270 def parse(self, strOSGR, name=NN): 

271 '''Parse an OSGR reference to a similar L{Osgr} instance. 

272 

273 @arg strOSGR: The OSGR reference (C{str}), see function L{parseOSGR}. 

274 @kwarg name: Optional instance name (C{str}), overriding this name. 

275 

276 @return: The similar instance (L{Osgr}) 

277 

278 @raise OSGRError: Invalid B{C{strOSGR}}. 

279 ''' 

280 return parseOSGR(strOSGR, Osgr=self.classof, name=name or self.name) 

281 

282 @property_RO 

283 def resolution(self): 

284 '''Get the OSGR resolution (C{meter}, power of 10) or C{0} if undefined. 

285 ''' 

286 return self._resolution 

287 

288 def toLatLon(self, LatLon=None, datum=_WGS84, kTM=False, eps=_10um, **LatLon_kwds): 

289 '''Convert this L{Osgr} coordinate to an (ellipsoidal) geodetic 

290 point. 

291 

292 @kwarg LatLon: Optional ellipsoidal class to return the 

293 geodetic point (C{LatLon}) or C{None}. 

294 @kwarg datum: Optional datum to convert to (L{Datum}, 

295 L{Ellipsoid}, L{Ellipsoid2}, L{Ellipsoid2} 

296 or L{a_f2Tuple}). 

297 @kwarg kTM: If C{True} use I{Karney}'s Krüger method from 

298 module L{ktm}, otherwise use the Ordnance 

299 Survey formulation (C{bool}). 

300 @kwarg eps: Tolerance for OS convergence (C{meter}). 

301 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword 

302 arguments, ignored if C{B{LatLon} is None}. 

303 

304 @return: A B{C{LatLon}} instance or if B{C{LatLon}} is C{None} 

305 a L{LatLonDatum3Tuple}C{(lat, lon, datum)}. 

306 

307 @note: While OS grid references are based on the OSGB36 datum, 

308 the Ordnance Survey have deprecated the use of OSGB36 for 

309 lat-/longitude coordinates (in favour of WGS84). Hence, 

310 this method returns WGS84 by default with OSGB36 as an 

311 option, U{see<https://www.OrdnanceSurvey.co.UK/blog/2014/12/2>}. 

312 

313 @note: The formulation implemented here due to Thomas, Redfearn, 

314 etc. is as published by the Ordnance Survey, but is 

315 inferior to Krüger as used by e.g. Karney 2011. 

316 

317 @raise OSGRError: No convergence. 

318 

319 @raise TypeError: If B{C{LatLon}} is not ellipsoidal or B{C{datum}} 

320 is invalid or conversion to B{C{datum}} failed. 

321 

322 @example: 

323 

324 >>> from pygeodesy import Datums, ellipsoidalVincenty as eV, Osgr 

325 >>> g = Osgr(651409.903, 313177.270) 

326 >>> p = g.toLatLon(eV.LatLon) # 52°39′28.723″N, 001°42′57.787″E 

327 >>> p = g.toLatLon(eV.LatLon, kTM=True) # 52°39′28.723″N, 001°42′57.787″E 

328 >>> # to obtain (historical) OSGB36 lat-/longitude point 

329 >>> p = g.toLatLon(eV.LatLon, datum=Datums.OSGB36) # 52°39′27.253″N, 001°43′04.518″E 

330 ''' 

331 NG = _NG 

332 if kTM: 

333 r = NG.reverse(self) 

334 

335 elif self._latlon is None: 

336 e0 = self.easting - NG.eas0 

337 n0 = m = self.northing - NG.nor0 

338 

339 _M = NG.Mabcd0 

340 a0 = NG.a0 

341 a = NG.phi0 

342 _A = Fsum(a).fsum_ 

343 for self._iteration in range(1, _TRIPS): 

344 a = _A(m / a0) 

345 m = n0 - _M(a) # meridional arc 

346 if fabs(m) < eps: 

347 break 

348 else: # PYCHOK no cover 

349 t = str(self) 

350 t = Fmt.PAREN(self.classname, repr(t)) 

351 t = _DOT_(t, self.toLatLon.__name__) 

352 t = unstr(t, eps=eps, kTM=kTM) 

353 raise OSGRError(Fmt.no_convergence(m), txt=t) 

354 

355 sa, ca, ta = sincostan3(a) 

356 v, v_r, n2 = NG.nu_rho_eta3(sa) 

357 

358 ta2 = ta**2 

359 ta4 = ta2**2 

360 

361 ta *= v_r / 2 

362 d = e0 / v 

363 d2 = d**2 

364 

365 a = (d2 * ta * (-1 + # Horner-like 

366 d2 / 12 * (Fsum( 5, 3 * ta2, -9 * ta2 * n2, n2) - 

367 d2 / 30 * Fsum(61, 90 * ta2, 45 * ta4)))).fsum_(a) 

368 

369 b = (d / ca * (1 - # Horner-like 

370 d2 / 6 * (Fsum(v_r, 2 * ta2) - 

371 d2 / 20 * (Fsum( 5, 28 * ta2, 24 * ta4) + 

372 d2 / 42 * Fsum(61, 662 * ta2, 1320 * ta4, 

373 720 * ta2 * ta4))))).fsum_(NG.lam0) 

374 

375 r = _LLEB(degrees90(a), degrees180(b), datum=self.datum, name=self.name) 

376 r._iteration = self._iteration # only ellipsoidal LatLon 

377 self._latlon = r 

378 else: 

379 r = self._latlon 

380 

381 return _ll2LatLon3(r, LatLon, datum, LatLon_kwds) 

382 

383 @Property_RO 

384 def scale0(self): 

385 '''Get the C{OS National Grid} central scale (C{scalar}). 

386 ''' 

387 return _NG.k0 

388 

389 def toRepr(self, GD=None, fmt=Fmt.SQUARE, sep=_COMMASPACE_, **prec): # PYCHOK expected 

390 '''Return a string representation of this L{Osgr} coordinate. 

391 

392 @kwarg GD: If C{bool}, in- or exclude the 2-letter grid designation and get 

393 the new B{C{prec}} behavior, otherwise if C{None}, default to the 

394 DEPRECATED definition C{B{prec}=5} I{for backward compatibility}. 

395 @kwarg fmt: Enclosing backets format (C{str}). 

396 @kwarg sep: Separator to join (C{str}) or C{None} to return an unjoined 2- or 

397 3-C{tuple} of C{str}s. 

398 @kwarg prec: Precison C{B{prec}=0}, the number of I{decimal} digits (C{int}) or 

399 if negative, the number of I{units to drop}, like MGRS U{PRECISION 

400 <https://GeographicLib.SourceForge.io/C++/doc/GeoConvert.1.html#PRECISION>}. 

401 

402 @return: This OSGR as (C{str}), C{"[G:GD, E:meter, N:meter]"} or if C{B{GD}=False} 

403 C{"[OSGR:easting,northing]"} or C{B{GD}=False} and C{B{prec} > 0} if 

404 C{"[OSGR:easting.d,northing.d]"}. 

405 

406 @note: OSGR easting and northing values are truncated, not rounded. 

407 

408 @raise OSGRError: If C{B{GD} not in (None, True, False)} or if C{B{prec} < -4} 

409 and C{B{GD}=False}. 

410 

411 @raise ValueError: Invalid B{C{prec}}. 

412 ''' 

413 GD, prec = _GD_prec2(GD, fmt=fmt, sep=sep, **prec) 

414 

415 if GD: 

416 t = self.toStr(GD=True, prec=prec, sep=None) 

417 t = _xzipairs('GEN', t, sep=sep, fmt=fmt) 

418 else: 

419 t = _COLON_(_OSGR_, self.toStr(GD=False, prec=prec)) 

420 if fmt: 

421 t = fmt % (t,) 

422 return t 

423 

424 def toStr(self, GD=None, sep=NN, **prec): # PYCHOK expected 

425 '''Return this L{Osgr} coordinate as a string. 

426 

427 @kwarg GD: If C{bool}, in- or exclude the 2-letter grid designation and get 

428 the new B{C{prec}} behavior, otherwise if C{None}, default to the 

429 DEPRECATED definition C{B{prec}=5} I{for backward compatibility}. 

430 @kwarg sep: Separator to join (C{str}) or C{None} to return an unjoined 2- or 

431 3-C{tuple} of C{str}s. 

432 @kwarg prec: Precison C{B{prec}=0}, the number of I{decimal} digits (C{int}) or 

433 if negative, the number of I{units to drop}, like MGRS U{PRECISION 

434 <https://GeographicLib.SourceForge.io/C++/doc/GeoConvert.1.html#PRECISION>}. 

435 

436 @return: This OSGR as (C{str}), C{"GD meter meter"} or if C{B{GD}=False} 

437 C{"easting,northing"} or if C{B{GD}=False} and C{B{prec} > 0} 

438 C{"easting.d,northing.d"} 

439 

440 @note: OSGR easting and northing values are truncated, not rounded. 

441 

442 @raise OSGRError: If C{B{GD} not in (None, True, False)} or if C{B{prec} 

443 < -4} and C{B{GD}=False}. 

444 

445 @raise ValueError: Invalid B{C{prec}}. 

446 

447 @example: 

448 

449 >>> r = Osgr(651409, 313177) 

450 >>> str(r) # 'TG 5140 1317' 

451 >>> r.toStr() # 'TG5140913177' 

452 >>> r.toStr(GD=False) # '651409,313177' 

453 ''' 

454 def _i2c(i): 

455 if i > 7: 

456 i += 1 

457 return chr(_ord_A + i) 

458 

459 GD, prec = _GD_prec2(GD, sep=sep, **prec) 

460 

461 if GD: 

462 E, e = divmod(self.easting, _100km) 

463 N, n = divmod(self.northing, _100km) 

464 E, N = int(E), int(N) 

465 if 0 > E or E > 6 or \ 

466 0 > N or N > 12: 

467 raise OSGRError(E=E, e=e, N=N, n=n, prec=prec, sep=sep) 

468 N = 19 - N 

469 EN = _i2c( N - (N % 5) + (E + 10) // 5) + \ 

470 _i2c((N * 5) % 25 + (E % 5)) 

471 t = enstr2(e, n, prec, EN) 

472 s = sep 

473 

474 elif prec <= -_EN_WIDE: 

475 raise OSGRError(GD=GD, prec=prec, sep=sep) 

476 else: 

477 t = enstr2(self.easting, self.northing, prec, dot=True, 

478 wide=_EN_WIDE + 1) 

479 s = sep if sep is None else (sep or _COMMA_) 

480 

481 return t if s is None else s.join(t) 

482 

483 

484def _GD_prec2(GD, **prec_et_al): 

485 '''(INTERNAL) Handle C{prec} backward compatibility. 

486 ''' 

487 if GD is None: # old C{prec} 5+ or neg 

488 prec = _xkwds_get(prec_et_al, prec=_EN_WIDE) 

489 GD = prec > 0 

490 prec = (prec - _EN_WIDE) if GD else -prec 

491 elif isbool(GD): 

492 prec = _xkwds_get(prec_et_al, prec=0) 

493 else: 

494 raise OSGRError(GD=GD, **prec_et_al) 

495 return GD, prec 

496 

497 

498def _ll2datum(ll, datum, name): 

499 '''(INTERNAL) Convert datum if needed. 

500 ''' 

501 if datum: 

502 try: 

503 if ll.datum != datum: 

504 ll = ll.toDatum(datum) 

505 except (AttributeError, TypeError, ValueError) as x: 

506 raise _TypeError(cause=x, datum=datum.name, **{name: ll}) 

507 return ll 

508 

509 

510def _ll2LatLon3(ll, LatLon, datum, LatLon_kwds): 

511 '''(INTERNAL) Convert C{ll} to C{LatLon} 

512 ''' 

513 if LatLon is None: 

514 r = _ll2datum(ll, datum, LatLonDatum3Tuple.__name__) 

515 r = LatLonDatum3Tuple(r.lat, r.lon, r.datum) 

516 else: # must be ellipsoidal 

517 _xsubclassof(_LLEB, LatLon=LatLon) 

518 r = _ll2datum(ll, datum, LatLon.__name__) 

519 r = LatLon(r.lat, r.lon, datum=r.datum, **LatLon_kwds) 

520 if r._iteration != ll._iteration: 

521 r._iteration = ll._iteration 

522 return _xnamed(r, nameof(ll)) 

523 

524 

525def parseOSGR(strOSGR, Osgr=Osgr, name=NN, **Osgr_kwds): 

526 '''Parse a string representing an OS Grid Reference, consisting 

527 of C{"[GD] easting northing"}. 

528 

529 Accepts standard OS Grid References like "SU 387 148", with 

530 or without whitespace separators, from 2- up to 22-digit 

531 references, or all-numeric, comma-separated references in 

532 meters, for example "438700,114800". 

533 

534 @arg strOSGR: An OSGR coordinate (C{str}). 

535 @kwarg Osgr: Optional class to return the OSGR coordinate 

536 (L{Osgr}) or C{None}. 

537 @kwarg name: Optional B{C{Osgr}} name (C{str}). 

538 @kwarg Osgr_kwds: Optional, additional B{C{Osgr}} keyword 

539 arguments, ignored if C{B{Osgr} is None}. 

540 

541 @return: An (B{C{Osgr}}) instance or if B{C{Osgr}} is 

542 C{None} an L{EasNor2Tuple}C{(easting, northing)}. 

543 

544 @raise OSGRError: Invalid B{C{strOSGR}}. 

545 

546 @example: 

547 

548 >>> r = parseOSGR('TG5140913177') 

549 >>> str(r) # 'TG 51409 13177' 

550 >>> r = parseOSGR('TG 51409 13177') 

551 >>> r.toStr() # 'TG5140913177' 

552 >>> r = parseOSGR('651409,313177') 

553 >>> r.toStr(sep=' ') # 'TG 51409 13177' 

554 >>> r.toStr(GD=False) # '651409,313177' 

555 ''' 

556 def _c2i(G): 

557 g = ord(G.upper()) - _ord_A 

558 if g > 7: 

559 g -= 1 

560 if g < 0 or g > 25: 

561 raise ValueError 

562 return g 

563 

564 def _OSGR(strOSGR, Osgr, name): 

565 s = _splituple(strOSGR.strip()) 

566 p = len(s) 

567 if not p: 

568 raise ValueError 

569 g = s[0] 

570 if p == 2 and isfloat(g): # "easting,northing" 

571 e, n, m = _enstr2m3(*s, wide=_EN_WIDE + 1) 

572 

573 else: 

574 if p == 1: # "GReastingnorthing" 

575 s = halfs2(g[2:]) 

576 g = g[:2] 

577 elif p == 2: # "GReasting northing" 

578 s = g[2:], s[1] # for backward ... 

579 g = g[:2] # ... compatibility 

580 elif p != 3: 

581 raise ValueError 

582 else: # "GR easting northing" 

583 s = s[1:] 

584 

585 e, n = map(_c2i, g) 

586 n, m = divmod(n, 5) 

587 E = ((e - 2) % 5) * 5 + m 

588 N = 19 - (e // 5) * 5 - n 

589 if 0 > E or E > 6 or \ 

590 0 > N or N > 12: 

591 raise ValueError 

592 

593 e, n, m = _enstr2m3(*s, wide=_EN_WIDE) 

594 e += E * _100km 

595 n += N * _100km 

596 

597 if Osgr is None: 

598 _ = _MODS.osgr.Osgr(e, n, resolution=m) # validate 

599 r = EasNor2Tuple(e, n, name=name) 

600 else: 

601 r = Osgr(e, n, name=name, 

602 **_xkwds(Osgr_kwds, resolution=m)) 

603 return r 

604 

605 return _parseX(_OSGR, strOSGR, Osgr, name, 

606 strOSGR=strOSGR, Error=OSGRError) 

607 

608 

609def toOsgr(latlon, lon=None, kTM=False, datum=_WGS84, Osgr=Osgr, name=NN, # MCCABE 14 

610 **prec_Osgr_kwds): 

611 '''Convert a lat-/longitude point to an OSGR coordinate. 

612 

613 @arg latlon: Latitude (C{degrees}) or an (ellipsoidal) geodetic 

614 C{LatLon} point. 

615 @kwarg lon: Optional longitude in degrees (scalar or C{None}). 

616 @kwarg kTM: If C{True} use I{Karney}'s Krüger method from 

617 module L{ktm}, otherwise use the Ordnance 

618 Survey formulation (C{bool}). 

619 @kwarg datum: Optional datum to convert B{C{lat, lon}} from 

620 (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or 

621 L{a_f2Tuple}). 

622 @kwarg Osgr: Optional class to return the OSGR coordinate 

623 (L{Osgr}) or C{None}. 

624 @kwarg name: Optional B{C{Osgr}} name (C{str}). 

625 @kwarg prec_Osgr_kwds: Optional L{truncate} precision 

626 C{B{prec}=ndigits} and/or additional 

627 B{C{Osgr}} keyword arguments, ignored 

628 if C{B{Osgr} is None}. 

629 

630 @return: An (B{C{Osgr}}) instance or if B{C{Osgr}} is C{None} 

631 an L{EasNor2Tuple}C{(easting, northing)}. 

632 

633 @note: If L{isint}C{(B{prec})} both easting and northing are 

634 L{truncate}d to the given number of digits. 

635 

636 @raise OSGRError: Invalid B{C{latlon}} or B{C{lon}}. 

637 

638 @raise TypeError: Non-ellipsoidal B{C{latlon}} or invalid 

639 B{C{datum}}, B{C{Osgr}}, B{C{Osgr_kwds}} 

640 or conversion to C{Datums.OSGB36} failed. 

641 

642 @example: 

643 

644 >>> p = LatLon(52.65798, 1.71605) 

645 >>> r = toOsgr(p) # [G:TG, E:51409, N:13177] 

646 >>> # for conversion of (historical) OSGB36 lat-/longitude: 

647 >>> r = toOsgr(52.65757, 1.71791, datum=Datums.OSGB36) 

648 >>> # alternatively and using Krüger 

649 >>> r = p.toOsgr(kTM=True) # [G:TG, E:51409, N:13177] 

650 ''' 

651 def _prec_kwds2(prec=MISSING, **kwds): 

652 return prec, kwds 

653 

654 if lon is not None: 

655 try: 

656 lat, lon = _MODS.dms.parseDMS2(latlon, lon) 

657 latlon = _LLEB(lat, lon, datum=datum) 

658 except Exception as x: 

659 raise OSGRError(latlon=latlon, lon=lon, datum=datum, cause=x) 

660 elif not isinstance(latlon, _LLEB): 

661 raise _TypeError(latlon=latlon, txt=_not_(_ellipsoidal_)) 

662 elif not name: # use latlon.name 

663 name = nameof(latlon) 

664 

665 NG = _NG 

666 # convert latlon to OSGB36 first 

667 ll = _ll2datum(latlon, NG.datum, _latlon_) 

668 

669 if kTM: 

670 e, n = NG.forward2(ll) 

671 

672 else: 

673 try: 

674 a, b = ll.philam 

675 except AttributeError: 

676 a, b = map1(radians, ll.lat, ll.lon) 

677 

678 sa, ca, ta = sincostan3(a) 

679 v, v_r, n2 = NG.nu_rho_eta3(sa) 

680 

681 m0 = NG.Mabcd0(a) 

682 b -= NG.lam0 

683 t = b * sa * v / 2 

684 d = b * ca 

685 d2 = d**2 

686 

687 ta2 = -(ta**2) 

688 ta4 = ta2**2 

689 

690 e = (d * v * (1 + # Horner-like 

691 d2 / 6 * (Fsum(v_r, ta2) + 

692 d2 / 20 * Fsum(5, 18 * ta2, ta4, 14 * n2, 

693 58 * n2 * ta2)))).fsum_(NG.eas0) 

694 

695 n = (d * t * (1 + # Horner-like 

696 d2 / 12 * (Fsum( 5, ta2, 9 * n2) + 

697 d2 / 30 * Fsum(61, ta4, 58 * ta2)))).fsum_(m0, NG.nor0) 

698 

699 p, kwds = _prec_kwds2(**prec_Osgr_kwds) 

700 if p is not MISSING: 

701 e = truncate(e, p) 

702 n = truncate(n, p) 

703 

704 if Osgr is None: 

705 _ = _MODS.osgr.Osgr(e, n) # validate 

706 r = EasNor2Tuple(e, n) 

707 else: 

708 r = Osgr(e, n, **kwds) # datum=NG.datum 

709 if lon is None and isinstance(latlon, _LLEB): 

710 if kTM: 

711 r._latlonTM = latlon # XXX weakref(latlon)? 

712 else: 

713 r._latlon = latlon # XXX weakref(latlon)? 

714 return _xnamed(r, name or nameof(latlon)) 

715 

716 

717if __name__ == '__main__': 

718 

719 from pygeodesy.lazily import printf 

720 from random import random, seed 

721 from time import localtime 

722 

723 seed(localtime().tm_yday) 

724 

725 def _rnd(X, n): 

726 X -= 2 

727 d = set() 

728 while len(d) < n: 

729 r = 1 + int(random() * X) 

730 if r not in d: 

731 d.add(r) 

732 yield r 

733 

734 D = _NG.datum 

735 i = t = 0 

736 t1 = t2 = 0, 0, 0, 0 

737 for e in _rnd(_NG.easX, 256): 

738 for n in _rnd(_NG.norX, 512): 

739 p = False 

740 t += 1 

741 

742 g = Osgr(e, n) 

743 v = g.toLatLon(kTM=False, datum=D) 

744 k = g.toLatLon(kTM=True, datum=D) 

745 d = max(fabs(v.lat - k.lat), fabs(v.lon - k.lon)) 

746 if d > t1[2]: 

747 t1 = e, n, d, t 

748 p = True 

749 

750 ll = _LLEB((v.lat + k.lat) / 2, 

751 (v.lon + k.lon) / 2, datum=D) 

752 v = ll.toOsgr(kTM=False) 

753 k = ll.toOsgr(kTM=True) 

754 d = max(fabs(v.easting - k.easting), 

755 fabs(v.northing - k.northing)) 

756 if d > t2[2]: 

757 t2 = ll.lat, ll.lon, d, t 

758 p = True 

759 

760 if p: 

761 i += 1 

762 printf('%5d: %s %s', i, 

763 'll(%.2f, %.2f) %.3e %d' % t2, 

764 'en(%d, %d) %.3e %d' % t1) 

765 printf('%d total %s', t, D.name) 

766 

767# **) MIT License 

768# 

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

770# 

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

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

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

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

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

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

777# 

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

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

780# 

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

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

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

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

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

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

787# OTHER DEALINGS IN THE SOFTWARE. 

788 

789# % python3 -m pygeodesy.osgr 

790# 1: ll(53.42, -0.59) 4.672e-07 1 en(493496, 392519) 2.796e-11 1 

791# 2: ll(60.86, -0.28) 2.760e-05 2 en(493496, 1220986) 2.509e-10 2 

792# 3: ll(61.41, -0.25) 3.045e-05 13 en(493496, 1281644) 2.774e-10 13 

793# 4: ll(61.41, -0.25) 3.045e-05 13 en(493496, 1192797) 3.038e-10 20 

794# 5: ll(61.41, -0.25) 3.045e-05 13 en(493496, 1192249) 3.073e-10 120 

795# 6: ll(61.55, -0.24) 3.120e-05 160 en(493496, 1192249) 3.073e-10 120 

796# 7: ll(61.55, -0.24) 3.122e-05 435 en(493496, 1192249) 3.073e-10 120 

797# 8: ll(61.57, -0.24) 3.130e-05 473 en(493496, 1192249) 3.073e-10 120 

798# 9: ll(58.66, -8.56) 8.084e-04 513 en(19711, 993800) 3.020e-06 513 

799# 10: ll(52.83, -7.65) 8.156e-04 518 en(19711, 993800) 3.020e-06 513 

800# 11: ll(51.55, -7.49) 8.755e-04 519 en(19711, 993800) 3.020e-06 513 

801# 12: ll(60.20, -8.87) 9.439e-04 521 en(19711, 1165686) 4.318e-06 521 

802# 13: ll(60.45, -8.92) 9.668e-04 532 en(19711, 1194002) 4.588e-06 532 

803# 14: ll(61.17, -9.08) 1.371e-03 535 en(19711, 1274463) 5.465e-06 535 

804# 15: ll(61.31, -9.11) 1.463e-03 642 en(19711, 1290590) 5.663e-06 642 

805# 16: ll(61.35, -9.12) 1.488e-03 807 en(19711, 1294976) 5.718e-06 807 

806# 17: ll(61.38, -9.13) 1.510e-03 929 en(19711, 1298667) 5.765e-06 929 

807# 18: ll(61.11, -9.24) 1.584e-03 11270 en(10307, 1268759) 6.404e-06 11270 

808# 19: ll(61.20, -9.26) 1.650e-03 11319 en(10307, 1278686) 6.545e-06 11319 

809# 20: ll(61.23, -9.27) 1.676e-03 11383 en(10307, 1282514) 6.600e-06 11383 

810# 21: ll(61.36, -9.30) 1.776e-03 11437 en(10307, 1297037) 6.816e-06 11437 

811# 22: ll(61.38, -9.30) 1.789e-03 11472 en(10307, 1298889) 6.844e-06 11472 

812# 23: ll(61.25, -9.39) 1.885e-03 91137 en(4367, 1285831) 7.392e-06 91137 

813# 24: ll(61.32, -9.40) 1.944e-03 91207 en(4367, 1293568) 7.519e-06 91207 

814# 25: ll(61.34, -9.41) 1.963e-03 91376 en(4367, 1296061) 7.561e-06 91376 

815# 26: ll(61.37, -9.41) 1.986e-03 91595 en(4367, 1298908) 7.608e-06 91595 

816# 131072 total OSGB36