Coverage for pygeodesy/osgr.py: 97%

306 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-06-01 11:43 -0400

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, _xkwds_pop2 

37from pygeodesy.fmath import Fdot, fpowers 

38from pygeodesy.fsums import _Fsumf_ 

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 _name2__, _NamedBase, nameof 

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__ = '24.05.31' 

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 = (_Fsumf_(4, 4 * n, 5 * n2, 5 * n3) / 4, 

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

126 _Fsumf_( 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, resolution=0, **name): 

198 '''New L{Osgr} coordinate. 

199 

200 @arg easting: Easting from the OS C{National Grid} origin (C{meter}). 

201 @arg northing: Northing from the OS C{National Grid} origin (C{meter}). 

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

203 @kwarg resolution: Optional resolution (C{meter}), C{0} for default. 

204 @kwarg name: Optional C{B{name}=NN} (C{str}). 

205 

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

207 or B{C{datum}} not an C{Datums.OSGB36} equivalent. 

208 ''' 

209 if datum: # PYCHOK no cover 

210 try: 

211 self._datum = _ellipsoidal_datum(datum) 

212 if self.datum != _NG.datum: 

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

214 except (TypeError, ValueError) as x: 

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

216 

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

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

219 

220 if name: 

221 self.name = name 

222 if resolution: 

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

224 

225 def __str__(self): 

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

227 

228 @Property_RO 

229 def datum(self): 

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

231 ''' 

232 return self._datum 

233 

234 @Property_RO 

235 def easting(self): 

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

237 ''' 

238 return self._easting 

239 

240 @Property_RO 

241 def falsing0(self): 

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

243 ''' 

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

245 

246 @property_RO 

247 def iteration(self): 

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

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

250 ''' 

251 return self._iteration 

252 

253 @Property_RO 

254 def latlon0(self): 

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

256 ''' 

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

258 

259 @Property_RO 

260 def northing(self): 

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

262 ''' 

263 return self._northing 

264 

265 def parse(self, strOSGR, **name): 

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

267 

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

269 @kwarg name: Optional C{B{name}=NN} (C{str}), overriding this name. 

270 

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

272 

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

274 ''' 

275 return parseOSGR(strOSGR, Osgr=self.classof, name=self._name__(name)) 

276 

277 @property_RO 

278 def resolution(self): 

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

280 ''' 

281 return self._resolution 

282 

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

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

285 point. 

286 

287 @kwarg LatLon: Optional ellipsoidal class to return the 

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

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

290 L{Ellipsoid}, L{Ellipsoid2}, L{Ellipsoid2} 

291 or L{a_f2Tuple}). 

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

293 module L{ktm}, otherwise use the Ordnance 

294 Survey formulation (C{bool}). 

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

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

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

298 

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

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

301 

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

303 the Ordnance Survey have deprecated the use of OSGB36 for 

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

305 this method returns WGS84 by default with OSGB36 as an 

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

307 

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

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

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

311 

312 @raise OSGRError: No convergence. 

313 

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

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

316 ''' 

317 NG = _NG 

318 if kTM: 

319 r = NG.reverse(self) 

320 

321 elif self._latlon is None: 

322 _F = _Fsumf_ 

323 e0 = self.easting - NG.eas0 

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

325 

326 a0 = NG.a0 

327 _M = NG.Mabcd0 

328 a = NG.phi0 

329 _a = fabs 

330 _A = _F(a).fsum_ 

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

332 a = _A(m / a0) 

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

334 if _a(m) < eps: 

335 break 

336 else: # PYCHOK no cover 

337 t = str(self) 

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

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

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

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

342 

343 sa, ca, ta = sincostan3(a) 

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

345 

346 ta2 = ta**2 

347 ta4 = ta2**2 

348 

349 ta *= v_r / 2 

350 d = e0 / v 

351 d2 = d**2 

352 

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

354 d2 / 12 * (_F( 5, 3 * ta2, -9 * ta2 * n2, n2) - 

355 d2 / 30 * _F(61, 90 * ta2, 45 * ta4)))).fsum_(a) 

356 

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

358 d2 / 6 * (_F(v_r, 2 * ta2) - 

359 d2 / 20 * (_F( 5, 28 * ta2, 24 * ta4) + 

360 d2 / 42 * _F(61, 662 * ta2, 1320 * ta4, 

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

362 

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

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

365 self._latlon = r 

366 else: 

367 r = self._latlon 

368 

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

370 

371 @Property_RO 

372 def scale0(self): 

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

374 ''' 

375 return _NG.k0 

376 

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

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

379 

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

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

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

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

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

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

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

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

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

389 

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

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

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

393 

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

395 

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

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

398 

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

400 ''' 

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

402 

403 if GD: 

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

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

406 else: 

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

408 if fmt: 

409 t = fmt % (t,) 

410 return t 

411 

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

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

414 

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

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

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

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

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

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

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

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

423 

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

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

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

427 

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

429 

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

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

432 

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

434 ''' 

435 def _i2c(i): 

436 if i > 7: 

437 i += 1 

438 return chr(_ord_A + i) 

439 

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

441 

442 if GD: 

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

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

445 E, N = int(E), int(N) 

446 if 0 > E or E > 6 or \ 

447 0 > N or N > 12: 

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

449 N = 19 - N 

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

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

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

453 s = sep 

454 

455 elif prec <= -_EN_WIDE: 

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

457 else: 

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

459 wide=_EN_WIDE + 1) 

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

461 

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

463 

464 

465def _GD_prec2(GD, **prec_et_al): 

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

467 ''' 

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

469 prec = _xkwds_get(prec_et_al, prec=_EN_WIDE) 

470 GD = prec > 0 

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

472 elif isbool(GD): 

473 prec = _xkwds_get(prec_et_al, prec=0) 

474 else: 

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

476 return GD, prec 

477 

478 

479def _ll2datum(ll, datum, name): 

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

481 ''' 

482 if datum: 

483 try: 

484 if ll.datum != datum: 

485 ll = ll.toDatum(datum) 

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

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

488 return ll 

489 

490 

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

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

493 ''' 

494 n = nameof(ll) 

495 if LatLon is None: 

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

497 r = LatLonDatum3Tuple(r.lat, r.lon, r.datum, name=n) 

498 else: # must be ellipsoidal 

499 _xsubclassof(_LLEB, LatLon=LatLon) 

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

501 r = LatLon(r.lat, r.lon, datum=r.datum, **_xkwds(LatLon_kwds, name=n)) 

502 if r._iteration != ll._iteration: 

503 r._iteration = ll._iteration 

504 return r 

505 

506 

507def parseOSGR(strOSGR, Osgr=Osgr, **name_Osgr_kwds): 

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

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

510 

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

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

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

514 meters, for example "438700,114800". 

515 

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

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

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

519 @kwarg name_Osgr_kwds: Optional C{B{name}=NN} (C{str}) and 

520 optional, additional B{C{Osgr}} keyword arguments, 

521 ignored if C{B{Osgr} is None}. 

522 

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

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

525 

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

527 ''' 

528 def _c2i(G): 

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

530 if g > 7: 

531 g -= 1 

532 if g < 0 or g > 25: 

533 raise ValueError 

534 return g 

535 

536 def _OSGR(strOSGR, Osgr, kwds): 

537 s = _splituple(strOSGR.strip()) 

538 p = len(s) 

539 if not p: 

540 raise ValueError 

541 g = s[0] 

542 if p == 2 and isfloat(g, both=True): # "easting,northing" 

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

544 

545 else: 

546 if p == 1: # "GReastingnorthing" 

547 s = halfs2(g[2:]) 

548 g = g[:2] 

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

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

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

552 elif p != 3: 

553 raise ValueError 

554 else: # "GR easting northing" 

555 s = s[1:] 

556 

557 e, n = map(_c2i, g) 

558 n, m = divmod(n, 5) 

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

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

561 if 0 > E or E > 6 or \ 

562 0 > N or N > 12: 

563 raise ValueError 

564 

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

566 e += E * _100km 

567 n += N * _100km 

568 

569 name, kwds = _name2__(**kwds) 

570 if Osgr is None: 

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

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

573 else: 

574 r = Osgr(e, n, name=name, **_xkwds(kwds, resolution=m)) 

575 return r 

576 

577 return _parseX(_OSGR, strOSGR, Osgr, name_Osgr_kwds, 

578 strOSGR=strOSGR, Error=OSGRError) 

579 

580 

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

582 **prec_name_Osgr_kwds): 

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

584 

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

586 C{LatLon} point. 

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

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

589 module L{ktm}, otherwise use the Ordnance 

590 Survey formulation (C{bool}). 

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

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

593 L{a_f2Tuple}). 

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

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

596 @kwarg prec_name_Osgr_kwds: Optional C{B{name}=NN} (C{str}), 

597 optional L{truncate} precision C{B{prec}=ndigits} 

598 and additional B{C{Osgr}} keyword arguments, 

599 ignored if C{B{Osgr} is None}. 

600 

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

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

603 

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

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

606 

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

608 

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

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

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

612 ''' 

613 if lon is not None: 

614 try: 

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

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

617 except Exception as x: 

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

619 elif not isinstance(latlon, _LLEB): 

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

621 

622 NG = _NG 

623 # convert latlon to OSGB36 first 

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

625 

626 if kTM: 

627 e, n = NG.forward2(ll) 

628 

629 else: 

630 try: 

631 a, b = ll.philam 

632 except AttributeError: 

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

634 

635 sa, ca, ta = sincostan3(a) 

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

637 

638 m0 = NG.Mabcd0(a) 

639 b -= NG.lam0 

640 t = b * sa * v / 2 

641 d = b * ca 

642 d2 = d**2 

643 

644 ta2 = -(ta**2) 

645 ta4 = ta2**2 

646 

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

648 d2 / 6 * (_Fsumf_(v_r, ta2) + 

649 d2 / 20 * _Fsumf_(5, 18 * ta2, ta4, 14 * n2, 

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

651 

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

653 d2 / 12 * (_Fsumf_( 5, ta2, 9 * n2) + 

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

655 

656 t, kwds = _name2__(prec_name_Osgr_kwds, _or_nameof=latlon) 

657 if kwds: 

658 p, kwds = _xkwds_pop2(kwds, prec=MISSING) 

659 if p is not MISSING: 

660 e = truncate(e, p) 

661 n = truncate(n, p) 

662 

663 if Osgr is None: 

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

665 r = EasNor2Tuple(e, n, name=t) 

666 else: 

667 r = Osgr(e, n, name=t, **kwds) # datum=NG.datum 

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

669 if kTM: 

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

671 else: 

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

673 return r 

674 

675 

676if __name__ == '__main__': 

677 

678 from pygeodesy import printf 

679 from random import random, seed 

680 from time import localtime 

681 

682 seed(localtime().tm_yday) 

683 

684 def _rnd(X, n): 

685 X -= 2 

686 d = set() 

687 while len(d) < n: 

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

689 if r not in d: 

690 d.add(r) 

691 yield r 

692 

693 D = _NG.datum 

694 i = t = 0 

695 t1 = t2 = 0, 0, 0, 0 

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

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

698 p = False 

699 t += 1 

700 

701 g = Osgr(e, n) 

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

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

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

705 if d > t1[2]: 

706 t1 = e, n, d, t 

707 p = True 

708 

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

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

711 v = ll.toOsgr(kTM=False) 

712 k = ll.toOsgr(kTM=True) 

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

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

715 if d > t2[2]: 

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

717 p = True 

718 

719 if p: 

720 i += 1 

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

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

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

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

725 

726# **) MIT License 

727# 

728# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved. 

729# 

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

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

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

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

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

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

736# 

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

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

739# 

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

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

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

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

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

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

746# OTHER DEALINGS IN THE SOFTWARE. 

747 

748# % python3 -m pygeodesy.osgr 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

775# 131072 total OSGB36