Coverage for pygeodesy/osgr.py : 97%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
# -*- coding: utf-8 -*-
<https://www.OrdnanceSurvey.co.UK/documents/resources/guide-to-nationalgrid.pdf>}.
Classes L{Osgr} and L{OSGRError} and functions L{parseOSGR} and L{toOsgr}.
A pure Python implementation, transcoded from I{Chris Veness}' JavaScript originals U{OS National Grid <https://www.Movable-Type.co.UK/scripts/latlong-os-gridref.html>} and U{Module osgridref <https://www.Movable-Type.co.UK/scripts/geodesy/docs/module-osgridref.html>} and I{Charles Karney}'s C++ class U{OSGB<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1OSGB.html>}.
OSGR provides geocoordinate references for UK mapping purposes, converted in 2015 to work with the C{WGS84} or the original C{OSGB36} datum. In addition, this implementation includes both the OS recommended and the Krüger-based method to convert between OSGR and geodetic coordinates (with keyword argument C{kTM} of function L{toOsgr}, method L{Osgr.toLatLon} and method C{toOsgr} of any ellipsoidal C{LatLon} class).
See U{Transverse Mercator: Redfearn series<https://WikiPedia.org/wiki/Transverse_Mercator:_Redfearn_series>}, Karney's U{"Transverse Mercator with an accuracy of a few nanometers", 2011<https://Arxiv.org/pdf/1002.1417v3.pdf>} (building on U{"Konforme Abbildung des Erdellipsoids in der Ebene", 1912<https://bib.GFZ-Potsdam.DE/pub/digi/krueger2.pdf>}, U{"Die Mathematik der Gauß-Krueger-Abbildung", 2006<https://DE.WikiPedia.org/wiki/Gauß-Krüger-Koordinatensystem>}, U{A Guide<https://www.OrdnanceSurvey.co.UK/documents/resources/guide-coordinate-systems-great-britain.pdf>} and U{Ordnance Survey National Grid<https://WikiPedia.org/wiki/Ordnance_Survey_National_Grid>}. ''' # make sure int/int division yields float quotient, see .basics
# from pygeodesy.dms import parseDMS2 # _MODS _xkwds, _xkwds_get # from pygeodesy.fsums import Fsum # from .fmath _COMMASPACE_, _DOT_, _ellipsoidal_, \ _latlon_, _not_, _SPACE_, _splituple LatLonDatum3Tuple _resolution10, unstr, _xzipairs Phi_, Scalar, _10um, _100km
'''Ordnance Survey National Grid parameters. '''
# U{Forward<https://GeographicLib.SourceForge.io/C++/doc/OSGB_8hpp_source.html>}
# <https://GeographicLib.SourceForge.io/C++/doc/OSGB_8hpp_source.html>
# <https://GeographicLib.SourceForge.io/C++/doc/OSGB_8cpp_source.html>
Fsum( 24 * n, 24 * n2, 21 * n3) / 8, Fsum( 15 * n2, 15 * n3) / 8, (35 * n3 / 24))
sin(s * 2) * cos(c * 2), -sin(s * 3) * cos(c * 3))
# <https://GeographicLib.SourceForge.io/C++/doc/OSGB_8cpp_source.html>
# rho = .a0 * E.e21 / s**1.5 == v * E.e21 / s # r = v * E.e21 / s # == rho, meridional roc # nu / rho == v / (v * E.e21 / s) == s / E.e21 == ...
# U{Reverse<https://GeographicLib.SourceForge.io/C++/doc/OSGB_8hpp_source.html>}
'''Error raised for a L{parseOSGR}, L{Osgr} or other OSGR issue. '''
'''Ordnance Survey Grid References (OSGR) coordinates on the U{National Grid<https://www.OrdnanceSurvey.co.UK/ documents/resources/guide-to-nationalgrid.pdf>}. '''
resolution=0): '''New L{Osgr} coordinate.
@arg easting: Easting from the OS C{National Grid} origin (C{meter}). @arg northing: Northing from the OS C{National Grid} origin (C{meter}). @kwarg datum: Override default datum (C{Datums.OSGB36}). @kwarg name: Optional name (C{str}). @kwarg resolution: Optional resolution (C{meter}), C{0} for default.
@raise OSGRError: Invalid or negative B{C{easting}} or B{C{northing}} or B{C{datum}} not an C{Datums.OSGB36} equivalent. ''' if datum: # PYCHOK no cover try: self._datum = _ellipsoidal_datum(datum) if self.datum != _NG.datum: raise ValueError(_not_(_NG.datum.name, _equivalent_)) except (TypeError, ValueError) as x: raise OSGRError(datum=datum, txt=str(x))
self.name = name
'''Get the datum (L{Datum}). '''
'''Get the easting (C{meter}). '''
'''Get the C{OS National Grid} falsing (L{EasNor2Tuple}). ''' return EasNor2Tuple(_NG.eas0, _NG.nor0, name=_OSGR_)
'''Get the most recent C{Osgr.toLatLon} iteration number (C{int}) or C{None} if not available/applicable. '''
'''Get the C{OS National Grid} origin (L{LatLon2Tuple}). ''' return LatLon2Tuple(_NG.lat, _NG.lon0, name=_OSGR_)
'''Get the northing (C{meter}). '''
'''Parse an OSGR reference to a similar L{Osgr} instance.
@arg strOSGR: The OSGR reference (C{str}), see function L{parseOSGR}. @kwarg name: Optional instance name (C{str}), overriding this name.
@return: The similar instance (L{Osgr})
@raise OSGRError: Invalid B{C{strOSGR}}. '''
'''Get the OSGR resolution (C{meter}, power of 10) or C{0} if undefined. '''
'''Convert this L{Osgr} coordinate to an (ellipsoidal) geodetic point.
@kwarg LatLon: Optional ellipsoidal class to return the geodetic point (C{LatLon}) or C{None}. @kwarg datum: Optional datum to convert to (L{Datum}, L{Ellipsoid}, L{Ellipsoid2}, L{Ellipsoid2} or L{a_f2Tuple}). @kwarg kTM: If C{True} use I{Karney}'s Krüger method from module L{ktm}, otherwise use the Ordnance Survey formulation (C{bool}). @kwarg eps: Tolerance for OS convergence (C{meter}). @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword arguments, ignored if C{B{LatLon} is None}.
@return: A B{C{LatLon}} instance or if B{C{LatLon}} is C{None} a L{LatLonDatum3Tuple}C{(lat, lon, datum)}.
@note: While OS grid references are based on the OSGB36 datum, the Ordnance Survey have deprecated the use of OSGB36 for lat-/longitude coordinates (in favour of WGS84). Hence, this method returns WGS84 by default with OSGB36 as an option, U{see<https://www.OrdnanceSurvey.co.UK/blog/2014/12/2>}.
@note: The formulation implemented here due to Thomas, Redfearn, etc. is as published by the Ordnance Survey, but is inferior to Krüger as used by e.g. Karney 2011.
@raise OSGRError: No convergence.
@raise TypeError: If B{C{LatLon}} is not ellipsoidal or B{C{datum}} is invalid or conversion to B{C{datum}} failed.
@example:
>>> from pygeodesy import Datums, ellipsoidalVincenty as eV, Osgr >>> g = Osgr(651409.903, 313177.270) >>> p = g.toLatLon(eV.LatLon) # 52°39′28.723″N, 001°42′57.787″E >>> p = g.toLatLon(eV.LatLon, kTM=True) # 52°39′28.723″N, 001°42′57.787″E >>> # to obtain (historical) OSGB36 lat-/longitude point >>> p = g.toLatLon(eV.LatLon, datum=Datums.OSGB36) # 52°39′27.253″N, 001°43′04.518″E '''
else: # PYCHOK no cover t = str(self) t = Fmt.PAREN(self.classname, repr(t)) t = _DOT_(t, self.toLatLon.__name__) t = unstr(t, eps=eps, kTM=kTM) raise OSGRError(Fmt.no_convergence(m), txt=t)
d2 / 12 * (Fsum( 5, 3 * ta2, -9 * ta2 * n2, n2) - d2 / 30 * Fsum(61, 90 * ta2, 45 * ta4)))).fsum_(a)
d2 / 6 * (Fsum(v_r, 2 * ta2) - d2 / 20 * (Fsum( 5, 28 * ta2, 24 * ta4) + d2 / 42 * Fsum(61, 662 * ta2, 1320 * ta4, 720 * ta2 * ta4))))).fsum_(NG.lam0)
else:
'''Get the C{OS National Grid} central scale (C{scalar}). ''' return _NG.k0
'''Return a string representation of this L{Osgr} coordinate.
@kwarg GD: If C{bool}, in- or exclude the 2-letter grid designation and get the new B{C{prec}} behavior, otherwise if C{None}, default to the DEPRECATED definition C{B{prec}=5} I{for backward compatibility}. @kwarg fmt: Enclosing backets format (C{str}). @kwarg sep: Separator to join (C{str}) or C{None} to return an unjoined 2- or 3-C{tuple} of C{str}s. @kwarg prec: Precison C{B{prec}=0}, the number of I{decimal} digits (C{int}) or if negative, the number of I{units to drop}, like MGRS U{PRECISION <https://GeographicLib.SourceForge.io/C++/doc/GeoConvert.1.html#PRECISION>}.
@return: This OSGR as (C{str}), C{"[G:GD, E:meter, N:meter]"} or if C{B{GD}=False} C{"[OSGR:easting,northing]"} or C{B{GD}=False} and C{B{prec} > 0} if C{"[OSGR:easting.d,northing.d]"}.
@note: OSGR easting and northing values are truncated, not rounded.
@raise OSGRError: If C{B{GD} not in (None, True, False)} or if C{B{prec} < -4} and C{B{GD}=False}.
@raise ValueError: Invalid B{C{prec}}. '''
else:
'''Return this L{Osgr} coordinate as a string.
@kwarg GD: If C{bool}, in- or exclude the 2-letter grid designation and get the new B{C{prec}} behavior, otherwise if C{None}, default to the DEPRECATED definition C{B{prec}=5} I{for backward compatibility}. @kwarg sep: Separator to join (C{str}) or C{None} to return an unjoined 2- or 3-C{tuple} of C{str}s. @kwarg prec: Precison C{B{prec}=0}, the number of I{decimal} digits (C{int}) or if negative, the number of I{units to drop}, like MGRS U{PRECISION <https://GeographicLib.SourceForge.io/C++/doc/GeoConvert.1.html#PRECISION>}.
@return: This OSGR as (C{str}), C{"GD meter meter"} or if C{B{GD}=False} C{"easting,northing"} or if C{B{GD}=False} and C{B{prec} > 0} C{"easting.d,northing.d"}
@note: OSGR easting and northing values are truncated, not rounded.
@raise OSGRError: If C{B{GD} not in (None, True, False)} or if C{B{prec} < -4} and C{B{GD}=False}.
@raise ValueError: Invalid B{C{prec}}.
@example:
>>> r = Osgr(651409, 313177) >>> str(r) # 'TG 5140 1317' >>> r.toStr() # 'TG5140913177' >>> r.toStr(GD=False) # '651409,313177' '''
0 > N or N > 12: raise OSGRError(E=E, e=e, N=N, n=n, prec=prec, sep=sep) _i2c((N * 5) % 25 + (E % 5))
raise OSGRError(GD=GD, prec=prec, sep=sep) else: wide=_EN_WIDE + 1)
'''(INTERNAL) Handle C{prec} backward compatibility. ''' else: raise OSGRError(GD=GD, **prec_et_al)
'''(INTERNAL) Convert datum if needed. ''' except (AttributeError, TypeError, ValueError) as x: raise _TypeError(txt=str(x), datum=datum.name, **{name: ll})
'''(INTERNAL) Convert C{ll} to C{LatLon} ''' else: # must be ellipsoidal
'''Parse a string representing an OS Grid Reference, consisting of C{"[GD] easting northing"}.
Accepts standard OS Grid References like "SU 387 148", with or without whitespace separators, from 2- up to 22-digit references, or all-numeric, comma-separated references in meters, for example "438700,114800".
@arg strOSGR: An OSGR coordinate (C{str}). @kwarg Osgr: Optional class to return the OSGR coordinate (L{Osgr}) or C{None}. @kwarg name: Optional B{C{Osgr}} name (C{str}). @kwarg Osgr_kwds: Optional, additional B{C{Osgr}} keyword arguments, ignored if C{B{Osgr} is None}.
@return: An (B{C{Osgr}}) instance or if B{C{Osgr}} is C{None} an L{EasNor2Tuple}C{(easting, northing)}.
@raise OSGRError: Invalid B{C{strOSGR}}.
@example:
>>> r = parseOSGR('TG5140913177') >>> str(r) # 'TG 51409 13177' >>> r = parseOSGR('TG 51409 13177') >>> r.toStr() # 'TG5140913177' >>> r = parseOSGR('651409,313177') >>> r.toStr(sep=' ') # 'TG 51409 13177' >>> r.toStr(GD=False) # '651409,313177' ''' raise ValueError
raise ValueError
else: raise ValueError else: # "GR easting northing"
0 > N or N > 12: raise ValueError
else: **_xkwds(Osgr_kwds, resolution=m))
strOSGR=strOSGR, Error=OSGRError)
**prec_Osgr_kwds): '''Convert a lat-/longitude point to an OSGR coordinate.
@arg latlon: Latitude (C{degrees}) or an (ellipsoidal) geodetic C{LatLon} point. @kwarg lon: Optional longitude in degrees (scalar or C{None}). @kwarg kTM: If C{True} use I{Karney}'s Krüger method from module L{ktm}, otherwise use the Ordnance Survey formulation (C{bool}). @kwarg datum: Optional datum to convert B{C{lat, lon}} from (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}). @kwarg Osgr: Optional class to return the OSGR coordinate (L{Osgr}) or C{None}. @kwarg name: Optional B{C{Osgr}} name (C{str}). @kwarg prec_Osgr_kwds: Optional L{truncate} precision C{B{prec}=ndigits} and/or additional B{C{Osgr}} keyword arguments, ignored if C{B{Osgr} is None}.
@return: An (B{C{Osgr}}) instance or if B{C{Osgr}} is C{None} an L{EasNor2Tuple}C{(easting, northing)}.
@note: If L{isint}C{(B{prec})} both easting and northing are L{truncate}d to the given number of digits.
@raise OSGRError: Invalid B{C{latlon}} or B{C{lon}}.
@raise TypeError: Non-ellipsoidal B{C{latlon}} or invalid B{C{datum}}, B{C{Osgr}}, B{C{Osgr_kwds}} or conversion to C{Datums.OSGB36} failed.
@example:
>>> p = LatLon(52.65798, 1.71605) >>> r = toOsgr(p) # [G:TG, E:51409, N:13177] >>> # for conversion of (historical) OSGB36 lat-/longitude: >>> r = toOsgr(52.65757, 1.71791, datum=Datums.OSGB36) >>> # alternatively and using Krüger >>> r = p.toOsgr(kTM=True) # [G:TG, E:51409, N:13177] '''
try: lat, lon = _MODS.dms.parseDMS2(latlon, lon) latlon = _LLEB(lat, lon, datum=datum) except Exception as x: raise OSGRError(latlon=latlon, lon=lon, datum=datum, txt=str(x)) raise _TypeError(latlon=latlon, txt=_not_(_ellipsoidal_))
# convert latlon to OSGB36 first
else: except AttributeError: a, b = map1(radians, ll.lat, ll.lon)
d2 / 6 * (Fsum(v_r, ta2) + d2 / 20 * Fsum(5, 18 * ta2, ta4, 14 * n2, 58 * n2 * ta2)))).fsum_(NG.eas0)
d2 / 12 * (Fsum( 5, ta2, 9 * n2) + d2 / 30 * Fsum(61, ta4, 58 * ta2)))).fsum_(m0, NG.nor0)
_ = _MODS.osgr.Osgr(e, n) # validate r = EasNor2Tuple(e, n) else: else:
if __name__ == '__main__':
from pygeodesy.lazily import printf from random import random, seed from time import localtime
seed(localtime().tm_yday)
def _rnd(X, n): X -= 2 d = set() while len(d) < n: r = 1 + int(random() * X) if r not in d: d.add(r) yield r
D = _NG.datum i = t = 0 t1 = t2 = 0, 0, 0, 0 for e in _rnd(_NG.easX, 256): for n in _rnd(_NG.norX, 512): p = False t += 1
g = Osgr(e, n) v = g.toLatLon(kTM=False, datum=D) k = g.toLatLon(kTM=True, datum=D) d = max(abs(v.lat - k.lat), abs(v.lon - k.lon)) if d > t1[2]: t1 = e, n, d, t p = True
ll = _LLEB((v.lat + k.lat) / 2, (v.lon + k.lon) / 2, datum=D) v = ll.toOsgr(kTM=False) k = ll.toOsgr(kTM=True) d = max(abs(v.easting - k.easting), abs(v.northing - k.northing)) if d > t2[2]: t2 = ll.lat, ll.lon, d, t p = True
if p: i += 1 printf('%5d: %s %s', i, 'll(%.2f, %.2f) %.3e %d' % t2, 'en(%d, %d) %.3e %d' % t1) printf('%d total %s', t, D.name)
# **) MIT License # # Copyright (C) 2016-2022 -- mrJean1 at Gmail -- All Rights Reserved. # # Permission is hereby granted, free of charge, to any person obtaining a # copy of this software and associated documentation files (the "Software"), # to deal in the Software without restriction, including without limitation # the rights to use, copy, modify, merge, publish, distribute, sublicense, # and/or sell copies of the Software, and to permit persons to whom the # Software is furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included # in all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS # OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR # OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, # ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR # OTHER DEALINGS IN THE SOFTWARE.
# % python3 -m pygeodesy.osgr # 1: ll(53.42, -0.59) 4.672e-07 1 en(493496, 392519) 2.796e-11 1 # 2: ll(60.86, -0.28) 2.760e-05 2 en(493496, 1220986) 2.509e-10 2 # 3: ll(61.41, -0.25) 3.045e-05 13 en(493496, 1281644) 2.774e-10 13 # 4: ll(61.41, -0.25) 3.045e-05 13 en(493496, 1192797) 3.038e-10 20 # 5: ll(61.41, -0.25) 3.045e-05 13 en(493496, 1192249) 3.073e-10 120 # 6: ll(61.55, -0.24) 3.120e-05 160 en(493496, 1192249) 3.073e-10 120 # 7: ll(61.55, -0.24) 3.122e-05 435 en(493496, 1192249) 3.073e-10 120 # 8: ll(61.57, -0.24) 3.130e-05 473 en(493496, 1192249) 3.073e-10 120 # 9: ll(58.66, -8.56) 8.084e-04 513 en(19711, 993800) 3.020e-06 513 # 10: ll(52.83, -7.65) 8.156e-04 518 en(19711, 993800) 3.020e-06 513 # 11: ll(51.55, -7.49) 8.755e-04 519 en(19711, 993800) 3.020e-06 513 # 12: ll(60.20, -8.87) 9.439e-04 521 en(19711, 1165686) 4.318e-06 521 # 13: ll(60.45, -8.92) 9.668e-04 532 en(19711, 1194002) 4.588e-06 532 # 14: ll(61.17, -9.08) 1.371e-03 535 en(19711, 1274463) 5.465e-06 535 # 15: ll(61.31, -9.11) 1.463e-03 642 en(19711, 1290590) 5.663e-06 642 # 16: ll(61.35, -9.12) 1.488e-03 807 en(19711, 1294976) 5.718e-06 807 # 17: ll(61.38, -9.13) 1.510e-03 929 en(19711, 1298667) 5.765e-06 929 # 18: ll(61.11, -9.24) 1.584e-03 11270 en(10307, 1268759) 6.404e-06 11270 # 19: ll(61.20, -9.26) 1.650e-03 11319 en(10307, 1278686) 6.545e-06 11319 # 20: ll(61.23, -9.27) 1.676e-03 11383 en(10307, 1282514) 6.600e-06 11383 # 21: ll(61.36, -9.30) 1.776e-03 11437 en(10307, 1297037) 6.816e-06 11437 # 22: ll(61.38, -9.30) 1.789e-03 11472 en(10307, 1298889) 6.844e-06 11472 # 23: ll(61.25, -9.39) 1.885e-03 91137 en(4367, 1285831) 7.392e-06 91137 # 24: ll(61.32, -9.40) 1.944e-03 91207 en(4367, 1293568) 7.519e-06 91207 # 25: ll(61.34, -9.41) 1.963e-03 91376 en(4367, 1296061) 7.561e-06 91376 # 26: ll(61.37, -9.41) 1.986e-03 91595 en(4367, 1298908) 7.608e-06 91595 # 131072 total OSGB36 |