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

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

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, Lamd, Lat, Lon, Northing, \

50 Phid, 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.11.06'

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 Lamd(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 Phid(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, 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 Survey

294 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 C{B{LatLon} is 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

465 def _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

479 def _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

491 def _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

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

508 '''Parse a string representing an OS Grid Reference, consisting of C{"[GD]

509 easting northing"}.

510

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

512 whitespace separators, from 2- up to 22-digit references, or all-numeric,

513 comma-separated references in meters, for example "438700,114800".

514

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

516 @kwarg Osgr: Optional class to return the OSGR coordinate (L{Osgr}) or C{None}.

517 @kwarg name_Osgr_kwds: Optional C{B{name}=NN} (C{str}) and optionally, additional

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

519

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

521 L{EasNor2Tuple}C{(easting, northing)}.

522

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

524 '''

525 def _c2i(G):

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

527 if g > 7:

528 g -= 1

529 if g < 0 or g > 25:

530 raise ValueError

531 return g

532

533 def _OSGR(strOSGR, Osgr, kwds):

534 s = _splituple(strOSGR.strip())

535 p = len(s)

536 if not p:

537 raise ValueError

538 g = s[0]

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

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

541

542 else:

543 if p == 1: # "GReastingnorthing"

544 s = halfs2(g[2:])

545 g = g[:2]

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

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

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

549 elif p != 3:

550 raise ValueError

551 else: # "GR easting northing"

552 s = s[1:]

553

554 e, n = map(_c2i, g)

555 n, m = divmod(n, 5)

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

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

558 if 0 > E or E > 6 or \

559 0 > N or N > 12:

560 raise ValueError

561

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

563 e += E * _100km

564 n += N * _100km

565

566 name, kwds = _name2__(**kwds)

567 if Osgr is None:

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

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

570 else:

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

572 return r

573

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

575 strOSGR=strOSGR, Error=OSGRError)

576

577

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

579 **prec_name_Osgr_kwds):

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

581

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

583 C{LatLon} point.

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

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

586 module L{ktm}, otherwise use the Ordnance Survey

587 formulation (C{bool}).

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

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

590 L{a_f2Tuple}).

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

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

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

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

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

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

597

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

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

600

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

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

603

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

605

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

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

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

609 '''

610 if lon is not None:

611 try:

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

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

614 except Exception as x:

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

616 elif not isinstance(latlon, _LLEB):

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

618

619 NG = _NG

620 # convert latlon to OSGB36 first

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

622

623 if kTM:

624 e, n = NG.forward2(ll)

625

626 else:

627 try:

628 a, b = ll.philam

629 except AttributeError:

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

631

632 sa, ca, ta = sincostan3(a)

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

634

635 m0 = NG.Mabcd0(a)

636 b -= NG.lam0

637 t = b * sa * v / 2

638 d = b * ca

639 d2 = d**2

640

641 ta2 = -(ta**2)

642 ta4 = ta2**2

643

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

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

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

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

648

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

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

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

652

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

654 if kwds:

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

656 if p is not MISSING:

657 e = truncate(e, p)

658 n = truncate(n, p)

659

660 if Osgr is None:

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

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

663 else:

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

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

666 if kTM:

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

668 else:

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

670 return r

671

672

673 if __name__ == '__main__':

674

675 from pygeodesy import printf

676 from random import random, seed

677 from time import localtime

678

679 seed(localtime().tm_yday)

680

681 def _rnd(X, n):

682 X -= 2

683 d = set()

684 while len(d) < n:

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

686 if r not in d:

687 d.add(r)

688 yield r

689

690 D = _NG.datum

691 i = t = 0

692 t1 = t2 = 0, 0, 0, 0

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

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

695 p = False

696 t += 1

697

698 g = Osgr(e, n)

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

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

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

702 if d > t1[2]:

703 t1 = e, n, d, t

704 p = True

705

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

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

708 v = ll.toOsgr(kTM=False)

709 k = ll.toOsgr(kTM=True)

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

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

712 if d > t2[2]:

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

714 p = True

715

716 if p:

717 i += 1

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

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

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

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

722

723# **) MIT License

724#

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

726#

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

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

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

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

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

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

733#

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

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

736#

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

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

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

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

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

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

743# OTHER DEALINGS IN THE SOFTWARE.

744

745# % python3 -m pygeodesy.osgr

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

772# 131072 total OSGB36