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

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

31 from pygeodesy.datums import Datums, _ellipsoidal_datum, _WGS84

32# from pygeodesy.dms import parseDMS2 # _MODS

33 from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase as _LLEB

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

35 _xkwds, _xkwds_get

36 from pygeodesy.fmath import Fdot, fpowers, Fsum

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

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

39 _COMMASPACE_, _DOT_, _ellipsoidal_, \

40 _latlon_, _not_, _SPACE_, _splituple

41 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS

42 from pygeodesy.named import _NamedBase, nameof, _xnamed

43 from pygeodesy.namedTuples import EasNor2Tuple, LatLon2Tuple, \

44 LatLonDatum3Tuple

45 from pygeodesy.props import Property_RO, property_RO

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

47 _resolution10, unstr, _xzipairs

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

49 Phi_, Scalar, _10um, _100km

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

51

52 from 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

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

178 class OSGRError(_ValueError):

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

180 '''

181 pass

182

183

184 class 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

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

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

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

524 def 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

608 def 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

716 if __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