1

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

3

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

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

8

9 A 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

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

13

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

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

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

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

18

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

20 Karney'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>},

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

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

24 and 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

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

28

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

30 _splituple, _xsubclassof

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

32 from pygeodesy.datums import Datums, _ellipsoidal_datum, _WGS84

33# from pygeodesy.dms import parseDMS2 # _MODS

34 from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase as _LLEB

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

36 _xkwds, _xkwds_get

37 from pygeodesy.fmath import Fdot, fpowers

38 from pygeodesy.fsums import _Fsumf_

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

40 _COMMASPACE_, _DOT_, _ellipsoidal_, \

41 _latlon_, _not_, _SPACE_

42 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS

43 from pygeodesy.named import _NamedBase, nameof, _xnamed

44 from pygeodesy.namedTuples import EasNor2Tuple, LatLon2Tuple, \

45 LatLonDatum3Tuple

46 from pygeodesy.props import Property_RO, property_RO

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

48 _resolution10, unstr, _xzipairs

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

50 Phi_, Scalar, _10um, _100km

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

52

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

54

55__all__ = _ALL_LAZY.osgr

56__version__ = '24.05.13'

57

58_equivalent_ = 'equivalent'

59_OSGR_ = 'OSGR'

60_ord_A = ord(_A_)

61_TRIPS = 33 # .toLatLon convergence

62

63

64 class _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

179 class OSGRError(_ValueError):

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

181 '''

182 pass

183

184

185 class 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 NG = _NG

323 if kTM:

324 r = NG.reverse(self)

325

326 elif self._latlon is None:

327 _F = _Fsumf_

328 e0 = self.easting - NG.eas0

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

330

331 a0 = NG.a0

332 _M = NG.Mabcd0

333 a = NG.phi0

334 _a = fabs

335 _A = _F(a).fsum_

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

337 a = _A(m / a0)

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

339 if _a(m) < eps:

340 break

341 else: # PYCHOK no cover

342 t = str(self)

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

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

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

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

347

348 sa, ca, ta = sincostan3(a)

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

350

351 ta2 = ta**2

352 ta4 = ta2**2

353

354 ta *= v_r / 2

355 d = e0 / v

356 d2 = d**2

357

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

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

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

361

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

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

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

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

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

367

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

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

370 self._latlon = r

371 else:

372 r = self._latlon

373

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

375

376 @Property_RO

377 def scale0(self):

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

379 '''

380 return _NG.k0

381

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

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

384

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

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

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

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

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

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

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

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

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

394

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

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

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

398

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

400

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

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

403

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

405 '''

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

407

408 if GD:

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

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

411 else:

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

413 if fmt:

414 t = fmt % (t,)

415 return t

416

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

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

419

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

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

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

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

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

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

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

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

428

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

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

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

432

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

434

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

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

437

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

439 '''

440 def _i2c(i):

441 if i > 7:

442 i += 1

443 return chr(_ord_A + i)

444

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

446

447 if GD:

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

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

450 E, N = int(E), int(N)

451 if 0 > E or E > 6 or \

452 0 > N or N > 12:

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

454 N = 19 - N

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

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

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

458 s = sep

459

460 elif prec <= -_EN_WIDE:

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

462 else:

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

464 wide=_EN_WIDE + 1)

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

466

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

468

469

470 def _GD_prec2(GD, **prec_et_al):

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

472 '''

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

474 prec = _xkwds_get(prec_et_al, prec=_EN_WIDE)

475 GD = prec > 0

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

477 elif isbool(GD):

478 prec = _xkwds_get(prec_et_al, prec=0)

479 else:

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

481 return GD, prec

482

483

484 def _ll2datum(ll, datum, name):

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

486 '''

487 if datum:

488 try:

489 if ll.datum != datum:

490 ll = ll.toDatum(datum)

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

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

493 return ll

494

495

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

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

498 '''

499 if LatLon is None:

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

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

502 else: # must be ellipsoidal

503 _xsubclassof(_LLEB, LatLon=LatLon)

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

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

506 if r._iteration != ll._iteration:

507 r._iteration = ll._iteration

508 return _xnamed(r, nameof(ll))

509

510

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

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

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

514

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

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

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

518 meters, for example "438700,114800".

519

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

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

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

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

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

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

526

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

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

529

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

531 '''

532 def _c2i(G):

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

534 if g > 7:

535 g -= 1

536 if g < 0 or g > 25:

537 raise ValueError

538 return g

539

540 def _OSGR(strOSGR, Osgr, name):

541 s = _splituple(strOSGR.strip())

542 p = len(s)

543 if not p:

544 raise ValueError

545 g = s[0]

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

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

548

549 else:

550 if p == 1: # "GReastingnorthing"

551 s = halfs2(g[2:])

552 g = g[:2]

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

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

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

556 elif p != 3:

557 raise ValueError

558 else: # "GR easting northing"

559 s = s[1:]

560

561 e, n = map(_c2i, g)

562 n, m = divmod(n, 5)

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

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

565 if 0 > E or E > 6 or \

566 0 > N or N > 12:

567 raise ValueError

568

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

570 e += E * _100km

571 n += N * _100km

572

573 if Osgr is None:

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

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

576 else:

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

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

579 return r

580

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

582 strOSGR=strOSGR, Error=OSGRError)

583

584

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

586 **prec_Osgr_kwds):

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

588

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

590 C{LatLon} point.

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

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

593 module L{ktm}, otherwise use the Ordnance

594 Survey formulation (C{bool}).

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

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

597 L{a_f2Tuple}).

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

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

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

601 @kwarg prec_Osgr_kwds: Optional L{truncate} precision

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

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

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

605

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

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

608

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

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

611

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

613

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

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

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

617 '''

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

619 return prec, kwds

620

621 if lon is not None:

622 try:

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

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

625 except Exception as x:

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

627 elif not isinstance(latlon, _LLEB):

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

629 elif not name: # use latlon.name

630 name = nameof(latlon)

631

632 NG = _NG

633 # convert latlon to OSGB36 first

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

635

636 if kTM:

637 e, n = NG.forward2(ll)

638

639 else:

640 try:

641 a, b = ll.philam

642 except AttributeError:

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

644

645 sa, ca, ta = sincostan3(a)

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

647

648 m0 = NG.Mabcd0(a)

649 b -= NG.lam0

650 t = b * sa * v / 2

651 d = b * ca

652 d2 = d**2

653

654 ta2 = -(ta**2)

655 ta4 = ta2**2

656

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

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

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

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

661

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

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

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

665

666 p, kwds = _prec_kwds2(**prec_Osgr_kwds)

667 if p is not MISSING:

668 e = truncate(e, p)

669 n = truncate(n, p)

670

671 if Osgr is None:

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

673 r = EasNor2Tuple(e, n)

674 else:

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

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

677 if kTM:

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

679 else:

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

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

682

683

684 if __name__ == '__main__':

685

686 from pygeodesy import printf

687 from random import random, seed

688 from time import localtime

689

690 seed(localtime().tm_yday)

691

692 def _rnd(X, n):

693 X -= 2

694 d = set()

695 while len(d) < n:

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

697 if r not in d:

698 d.add(r)

699 yield r

700

701 D = _NG.datum

702 i = t = 0

703 t1 = t2 = 0, 0, 0, 0

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

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

706 p = False

707 t += 1

708

709 g = Osgr(e, n)

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

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

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

713 if d > t1[2]:

714 t1 = e, n, d, t

715 p = True

716

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

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

719 v = ll.toOsgr(kTM=False)

720 k = ll.toOsgr(kTM=True)

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

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

723 if d > t2[2]:

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

725 p = True

726

727 if p:

728 i += 1

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

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

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

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

733

734# **) MIT License

735#

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

737#

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

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

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

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

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

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

744#

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

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

747#

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

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

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

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

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

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

754# OTHER DEALINGS IN THE SOFTWARE.

755

756# % python3 -m pygeodesy.osgr

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

783# 131072 total OSGB36