Coverage for pygeodesy/osgr.py: 97%

303 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-04-11 14:35 -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, _xsubclassof 

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

31from pygeodesy.datums import Datums, _ellipsoidal_datum, _WGS84 

32# from pygeodesy.dms import parseDMS2 # _MODS 

33from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase as _LLEB 

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

35 _xkwds, _xkwds_get 

36from pygeodesy.fmath import Fdot, fpowers, Fsum 

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

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

39 _COMMASPACE_, _DOT_, _ellipsoidal_, \ 

40 _latlon_, _not_, _SPACE_, _splituple 

41from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS 

42from pygeodesy.named import _NamedBase, nameof, _xnamed 

43from pygeodesy.namedTuples import EasNor2Tuple, LatLon2Tuple, \ 

44 LatLonDatum3Tuple 

45from pygeodesy.props import Property_RO, property_RO 

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

47 _resolution10, unstr, _xzipairs 

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

49 Phi_, Scalar, _10um, _100km 

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

51 

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

53 

54__all__ = _ALL_LAZY.osgr 

55__version__ = '23.03.19' 

56 

57_equivalent_ = 'equivalent' 

58_OSGR_ = 'OSGR' 

59_ord_A = ord(_A_) 

60_TRIPS = 33 # .toLatLon convergence 

61 

62 

63class _NG(object): 

64 '''Ordnance Survey National Grid parameters. 

65 ''' 

66 @Property_RO 

67 def a0(self): # equatoradius, scaled 

68 return self.ellipsoid.a * self.k0 

69 

70 @Property_RO 

71 def b0(self): # polaradius, scaled 

72 return self.ellipsoid.b * self.k0 

73 

74 @Property_RO 

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

76 return Datums.OSGB36 

77 

78 @Property_RO 

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

80 return Easting(4 * _100km) 

81 

82 @Property_RO 

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

84 return Easting(7 * _100km) 

85 

86 @Property_RO 

87 def ellipsoid(self): # ellipsoid, Airy130 

88 return self.datum.ellipsoid 

89 

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

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

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

93 e = t.easting + self.eas0 

94 n = t.northing + self.nor0ffset 

95 return e, n 

96 

97 @Property_RO 

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

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

100 _0_9998268 = (9998268 - 10000000) / 10000000 

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

102 

103 @Property_RO 

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

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

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

107 

108 @Property_RO 

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

110 return Lam_(self.lon0) 

111 

112 @Property_RO 

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

114 return Lat(49.0) 

115 

116 @Property_RO 

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

118 return Lon(_N_2_0) 

119 

120 @Property_RO 

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

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

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

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

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

126 (35 * n3 / 24)) 

127 return M 

128 

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

130 c = a + self.phi0 

131 s = a - self.phi0 

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

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

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

135 return float(R * self.b0) 

136 

137 @Property_RO 

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

139 return Northing(-_100km) 

140 

141 @Property_RO 

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

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

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

145 

146 @Property_RO 

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

148 return Northing(13 * _100km) 

149 

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

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

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

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

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

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

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

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

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

159 

160 @Property_RO 

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

162 return Phi_(self.lat0) 

163 

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

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

166 r = osgr._latlonTM 

167 if r is None: 

168 x = osgr.easting - self.eas0 

169 y = osgr.northing - self.nor0ffset 

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

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

172 osgr._latlonTM = r 

173 return r 

174 

175_NG = _NG() # PYCHOK singleton 

176 

177 

178class OSGRError(_ValueError): 

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

180 ''' 

181 pass 

182 

183 

184class Osgr(_NamedBase): 

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

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

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

188 ''' 

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

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

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

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

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

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

195 

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

197 resolution=0): 

198 '''New L{Osgr} coordinate. 

199 

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

201 origin (C{meter}). 

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

203 origin (C{meter}). 

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

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

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

207 C{0} for default. 

208 

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

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

211 C{Datums.OSGB36} equivalent. 

212 ''' 

213 if datum: # PYCHOK no cover 

214 try: 

215 self._datum = _ellipsoidal_datum(datum) 

216 if self.datum != _NG.datum: 

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

218 except (TypeError, ValueError) as x: 

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

220 

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

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

223 

224 if name: 

225 self.name = name 

226 if resolution: 

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

228 

229 def __str__(self): 

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

231 

232 @Property_RO 

233 def datum(self): 

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

235 ''' 

236 return self._datum 

237 

238 @Property_RO 

239 def easting(self): 

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

241 ''' 

242 return self._easting 

243 

244 @Property_RO 

245 def falsing0(self): 

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

247 ''' 

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

249 

250 @property_RO 

251 def iteration(self): 

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

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

254 ''' 

255 return self._iteration 

256 

257 @Property_RO 

258 def latlon0(self): 

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

260 ''' 

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

262 

263 @Property_RO 

264 def northing(self): 

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

266 ''' 

267 return self._northing 

268 

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

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

271 

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

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

274 

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

276 

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

278 ''' 

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

280 

281 @property_RO 

282 def resolution(self): 

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

284 ''' 

285 return self._resolution 

286 

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

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

289 point. 

290 

291 @kwarg LatLon: Optional ellipsoidal class to return the 

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

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

294 L{Ellipsoid}, L{Ellipsoid2}, L{Ellipsoid2} 

295 or L{a_f2Tuple}). 

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

297 module L{ktm}, otherwise use the Ordnance 

298 Survey formulation (C{bool}). 

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

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

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

302 

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

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

305 

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

307 the Ordnance Survey have deprecated the use of OSGB36 for 

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

309 this method returns WGS84 by default with OSGB36 as an 

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

311 

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

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

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

315 

316 @raise OSGRError: No convergence. 

317 

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

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

320 

321 @example: 

322 

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

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

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

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

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

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

329 ''' 

330 NG = _NG 

331 if kTM: 

332 r = NG.reverse(self) 

333 

334 elif self._latlon is None: 

335 e0 = self.easting - NG.eas0 

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

337 

338 _M = NG.Mabcd0 

339 a0 = NG.a0 

340 a = NG.phi0 

341 _A = Fsum(a).fsum_ 

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

343 a = _A(m / a0) 

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

345 if fabs(m) < eps: 

346 break 

347 else: # PYCHOK no cover 

348 t = str(self) 

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

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

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

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

353 

354 sa, ca, ta = sincostan3(a) 

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

356 

357 ta2 = ta**2 

358 ta4 = ta2**2 

359 

360 ta *= v_r / 2 

361 d = e0 / v 

362 d2 = d**2 

363 

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

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

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

367 

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

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

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

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

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

373 

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

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

376 self._latlon = r 

377 else: 

378 r = self._latlon 

379 

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

381 

382 @Property_RO 

383 def scale0(self): 

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

385 ''' 

386 return _NG.k0 

387 

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

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

390 

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

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

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

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

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

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

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

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

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

400 

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

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

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

404 

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

406 

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

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

409 

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

411 ''' 

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

413 

414 if GD: 

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

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

417 else: 

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

419 if fmt: 

420 t = fmt % (t,) 

421 return t 

422 

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

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

425 

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

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

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

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

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

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

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

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

434 

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

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

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

438 

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

440 

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

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

443 

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

445 

446 @example: 

447 

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

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

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

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

452 ''' 

453 def _i2c(i): 

454 if i > 7: 

455 i += 1 

456 return chr(_ord_A + i) 

457 

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

459 

460 if GD: 

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

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

463 E, N = int(E), int(N) 

464 if 0 > E or E > 6 or \ 

465 0 > N or N > 12: 

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

467 N = 19 - N 

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

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

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

471 s = sep 

472 

473 elif prec <= -_EN_WIDE: 

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

475 else: 

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

477 wide=_EN_WIDE + 1) 

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

479 

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

481 

482 

483def _GD_prec2(GD, **prec_et_al): 

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

485 ''' 

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

487 prec = _xkwds_get(prec_et_al, prec=_EN_WIDE) 

488 GD = prec > 0 

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

490 elif isbool(GD): 

491 prec = _xkwds_get(prec_et_al, prec=0) 

492 else: 

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

494 return GD, prec 

495 

496 

497def _ll2datum(ll, datum, name): 

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

499 ''' 

500 if datum: 

501 try: 

502 if ll.datum != datum: 

503 ll = ll.toDatum(datum) 

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

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

506 return ll 

507 

508 

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

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

511 ''' 

512 if LatLon is None: 

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

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

515 else: # must be ellipsoidal 

516 _xsubclassof(_LLEB, LatLon=LatLon) 

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

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

519 if r._iteration != ll._iteration: 

520 r._iteration = ll._iteration 

521 return _xnamed(r, nameof(ll)) 

522 

523 

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

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

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

527 

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

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

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

531 meters, for example "438700,114800". 

532 

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

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

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

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

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

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

539 

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

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

542 

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

544 

545 @example: 

546 

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

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

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

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

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

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

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

554 ''' 

555 def _c2i(G): 

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

557 if g > 7: 

558 g -= 1 

559 if g < 0 or g > 25: 

560 raise ValueError 

561 return g 

562 

563 def _OSGR(strOSGR, Osgr, name): 

564 s = _splituple(strOSGR.strip()) 

565 p = len(s) 

566 if not p: 

567 raise ValueError 

568 g = s[0] 

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

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

571 

572 else: 

573 if p == 1: # "GReastingnorthing" 

574 s = halfs2(g[2:]) 

575 g = g[:2] 

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

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

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

579 elif p != 3: 

580 raise ValueError 

581 else: # "GR easting northing" 

582 s = s[1:] 

583 

584 e, n = map(_c2i, g) 

585 n, m = divmod(n, 5) 

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

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

588 if 0 > E or E > 6 or \ 

589 0 > N or N > 12: 

590 raise ValueError 

591 

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

593 e += E * _100km 

594 n += N * _100km 

595 

596 if Osgr is None: 

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

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

599 else: 

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

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

602 return r 

603 

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

605 strOSGR=strOSGR, Error=OSGRError) 

606 

607 

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

609 **prec_Osgr_kwds): 

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

611 

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

613 C{LatLon} point. 

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

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

616 module L{ktm}, otherwise use the Ordnance 

617 Survey formulation (C{bool}). 

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

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

620 L{a_f2Tuple}). 

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

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

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

624 @kwarg prec_Osgr_kwds: Optional L{truncate} precision 

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

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

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

628 

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

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

631 

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

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

634 

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

636 

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

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

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

640 

641 @example: 

642 

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

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

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

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

647 >>> # alternatively and using Krüger 

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

649 ''' 

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

651 return prec, kwds 

652 

653 if lon is not None: 

654 try: 

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

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

657 except Exception as x: 

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

659 elif not isinstance(latlon, _LLEB): 

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

661 elif not name: # use latlon.name 

662 name = nameof(latlon) 

663 

664 NG = _NG 

665 # convert latlon to OSGB36 first 

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

667 

668 if kTM: 

669 e, n = NG.forward2(ll) 

670 

671 else: 

672 try: 

673 a, b = ll.philam 

674 except AttributeError: 

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

676 

677 sa, ca, ta = sincostan3(a) 

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

679 

680 m0 = NG.Mabcd0(a) 

681 b -= NG.lam0 

682 t = b * sa * v / 2 

683 d = b * ca 

684 d2 = d**2 

685 

686 ta2 = -(ta**2) 

687 ta4 = ta2**2 

688 

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

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

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

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

693 

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

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

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

697 

698 p, kwds = _prec_kwds2(**prec_Osgr_kwds) 

699 if p is not MISSING: 

700 e = truncate(e, p) 

701 n = truncate(n, p) 

702 

703 if Osgr is None: 

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

705 r = EasNor2Tuple(e, n) 

706 else: 

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

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

709 if kTM: 

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

711 else: 

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

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

714 

715 

716if __name__ == '__main__': 

717 

718 from pygeodesy.lazily import printf 

719 from random import random, seed 

720 from time import localtime 

721 

722 seed(localtime().tm_yday) 

723 

724 def _rnd(X, n): 

725 X -= 2 

726 d = set() 

727 while len(d) < n: 

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

729 if r not in d: 

730 d.add(r) 

731 yield r 

732 

733 D = _NG.datum 

734 i = t = 0 

735 t1 = t2 = 0, 0, 0, 0 

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

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

738 p = False 

739 t += 1 

740 

741 g = Osgr(e, n) 

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

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

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

745 if d > t1[2]: 

746 t1 = e, n, d, t 

747 p = True 

748 

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

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

751 v = ll.toOsgr(kTM=False) 

752 k = ll.toOsgr(kTM=True) 

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

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

755 if d > t2[2]: 

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

757 p = True 

758 

759 if p: 

760 i += 1 

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

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

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

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

765 

766# **) MIT License 

767# 

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

769# 

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

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

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

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

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

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

776# 

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

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

779# 

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

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

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

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

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

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

786# OTHER DEALINGS IN THE SOFTWARE. 

787 

788# % python3 -m pygeodesy.osgr 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

815# 131072 total OSGB36