Coverage for pygeodesy/osgr.py: 97%
303 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-03-31 10:52 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-03-31 10:52 -0400
2# -*- coding: utf-8 -*-
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>}.
7Classes L{Osgr} and L{OSGRError} and functions L{parseOSGR} and L{toOsgr}.
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>}.
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).
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
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
52from math import cos, fabs, radians, sin, sqrt
54__all__ = _ALL_LAZY.osgr
55__version__ = '23.03.19'
57_equivalent_ = 'equivalent'
58_OSGR_ = 'OSGR'
59_ord_A = ord(_A_)
60_TRIPS = 33 # .toLatLon convergence
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
70 @Property_RO
71 def b0(self): # polaradius, scaled
72 return self.ellipsoid.b * self.k0
74 @Property_RO
75 def datum(self): # datum, Airy130 ellipsoid
76 return Datums.OSGB36
78 @Property_RO
79 def eas0(self): # False origin easting (C{meter})
80 return Easting(4 * _100km)
82 @Property_RO
83 def easX(self): # easting [0..extent] (C{meter})
84 return Easting(7 * _100km)
86 @Property_RO
87 def ellipsoid(self): # ellipsoid, Airy130
88 return self.datum.ellipsoid
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
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...
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)
108 @Property_RO
109 def lam0(self): # True origin longitude C{radians}
110 return Lam_(self.lon0)
112 @Property_RO
113 def lat0(self): # True origin latitude, 49°N
114 return Lat(49.0)
116 @Property_RO
117 def lon0(self): # True origin longitude, 2°W
118 return Lon(_N_2_0)
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
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)
137 @Property_RO
138 def nor0(self): # False origin northing (C{meter})
139 return Northing(-_100km)
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
146 @Property_RO
147 def norX(self): # northing [0..extent] (C{meter})
148 return Northing(13 * _100km)
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
160 @Property_RO
161 def phi0(self): # True origin latitude C{radians}
162 return Phi_(self.lat0)
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
175_NG = _NG() # PYCHOK singleton
178class OSGRError(_ValueError):
179 '''Error raised for a L{parseOSGR}, L{Osgr} or other OSGR issue.
180 '''
181 pass
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})
196 def __init__(self, easting, northing, datum=None, name=NN,
197 resolution=0):
198 '''New L{Osgr} coordinate.
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.
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)
221 self._easting = Easting( easting, Error=OSGRError, high=_NG.easX)
222 self._northing = Northing(northing, Error=OSGRError, high=_NG.norX)
224 if name:
225 self.name = name
226 if resolution:
227 self._resolution = _resolution10(resolution, Error=OSGRError)
229 def __str__(self):
230 return self.toStr(GD=True, sep=_SPACE_)
232 @Property_RO
233 def datum(self):
234 '''Get the datum (L{Datum}).
235 '''
236 return self._datum
238 @Property_RO
239 def easting(self):
240 '''Get the easting (C{meter}).
241 '''
242 return self._easting
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_)
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
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_)
263 @Property_RO
264 def northing(self):
265 '''Get the northing (C{meter}).
266 '''
267 return self._northing
269 def parse(self, strOSGR, name=NN):
270 '''Parse an OSGR reference to a similar L{Osgr} instance.
272 @arg strOSGR: The OSGR reference (C{str}), see function L{parseOSGR}.
273 @kwarg name: Optional instance name (C{str}), overriding this name.
275 @return: The similar instance (L{Osgr})
277 @raise OSGRError: Invalid B{C{strOSGR}}.
278 '''
279 return parseOSGR(strOSGR, Osgr=self.classof, name=name or self.name)
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
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.
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}.
303 @return: A B{C{LatLon}} instance or if B{C{LatLon}} is C{None}
304 a L{LatLonDatum3Tuple}C{(lat, lon, datum)}.
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>}.
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.
316 @raise OSGRError: No convergence.
318 @raise TypeError: If B{C{LatLon}} is not ellipsoidal or B{C{datum}}
319 is invalid or conversion to B{C{datum}} failed.
321 @example:
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)
334 elif self._latlon is None:
335 e0 = self.easting - NG.eas0
336 n0 = m = self.northing - NG.nor0
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)
354 sa, ca, ta = sincostan3(a)
355 v, v_r, n2 = NG.nu_rho_eta3(sa)
357 ta2 = ta**2
358 ta4 = ta2**2
360 ta *= v_r / 2
361 d = e0 / v
362 d2 = d**2
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)
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)
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
380 return _ll2LatLon3(r, LatLon, datum, LatLon_kwds)
382 @Property_RO
383 def scale0(self):
384 '''Get the C{OS National Grid} central scale (C{scalar}).
385 '''
386 return _NG.k0
388 def toRepr(self, GD=None, fmt=Fmt.SQUARE, sep=_COMMASPACE_, **prec): # PYCHOK expected
389 '''Return a string representation of this L{Osgr} coordinate.
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>}.
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]"}.
405 @note: OSGR easting and northing values are truncated, not rounded.
407 @raise OSGRError: If C{B{GD} not in (None, True, False)} or if C{B{prec} < -4}
408 and C{B{GD}=False}.
410 @raise ValueError: Invalid B{C{prec}}.
411 '''
412 GD, prec = _GD_prec2(GD, fmt=fmt, sep=sep, **prec)
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
423 def toStr(self, GD=None, sep=NN, **prec): # PYCHOK expected
424 '''Return this L{Osgr} coordinate as a string.
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>}.
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"}
439 @note: OSGR easting and northing values are truncated, not rounded.
441 @raise OSGRError: If C{B{GD} not in (None, True, False)} or if C{B{prec}
442 < -4} and C{B{GD}=False}.
444 @raise ValueError: Invalid B{C{prec}}.
446 @example:
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)
458 GD, prec = _GD_prec2(GD, sep=sep, **prec)
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
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_)
480 return t if s is None else s.join(t)
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
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
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))
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"}.
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".
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}.
540 @return: An (B{C{Osgr}}) instance or if B{C{Osgr}} is
541 C{None} an L{EasNor2Tuple}C{(easting, northing)}.
543 @raise OSGRError: Invalid B{C{strOSGR}}.
545 @example:
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
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)
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:]
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
592 e, n, m = _enstr2m3(*s, wide=_EN_WIDE)
593 e += E * _100km
594 n += N * _100km
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
604 return _parseX(_OSGR, strOSGR, Osgr, name,
605 strOSGR=strOSGR, Error=OSGRError)
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.
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}.
629 @return: An (B{C{Osgr}}) instance or if B{C{Osgr}} is C{None}
630 an L{EasNor2Tuple}C{(easting, northing)}.
632 @note: If L{isint}C{(B{prec})} both easting and northing are
633 L{truncate}d to the given number of digits.
635 @raise OSGRError: Invalid B{C{latlon}} or B{C{lon}}.
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.
641 @example:
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
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)
664 NG = _NG
665 # convert latlon to OSGB36 first
666 ll = _ll2datum(latlon, NG.datum, _latlon_)
668 if kTM:
669 e, n = NG.forward2(ll)
671 else:
672 try:
673 a, b = ll.philam
674 except AttributeError:
675 a, b = map1(radians, ll.lat, ll.lon)
677 sa, ca, ta = sincostan3(a)
678 v, v_r, n2 = NG.nu_rho_eta3(sa)
680 m0 = NG.Mabcd0(a)
681 b -= NG.lam0
682 t = b * sa * v / 2
683 d = b * ca
684 d2 = d**2
686 ta2 = -(ta**2)
687 ta4 = ta2**2
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)
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)
698 p, kwds = _prec_kwds2(**prec_Osgr_kwds)
699 if p is not MISSING:
700 e = truncate(e, p)
701 n = truncate(n, p)
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))
716if __name__ == '__main__':
718 from pygeodesy.lazily import printf
719 from random import random, seed
720 from time import localtime
722 seed(localtime().tm_yday)
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
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
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
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
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)
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.
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