From 31ac581c3a09f4d042da4408582f9dd9ac9b8b69 Mon Sep 17 00:00:00 2001 From: Benny Malengier Date: Sat, 31 Jan 2009 23:09:27 +0000 Subject: [PATCH] Peter Landgren: 02659 support for RT90 in PlaceUtils.py svn: r11778 --- src/PlaceUtils.py | 136 ++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 132 insertions(+), 4 deletions(-) diff --git a/src/PlaceUtils.py b/src/PlaceUtils.py index 6eab09f88..54ffadede 100644 --- a/src/PlaceUtils.py +++ b/src/PlaceUtils.py @@ -30,6 +30,7 @@ # #------------------------------------------------------------------------- from gettext import gettext as _ +import math #------------------------------------------------------------------------- # @@ -285,6 +286,8 @@ def conv_lat_lon(latitude, longitude, format="D.D4"): i.e. ±DDMM.MMM±DDDMM.MMM 'ISO-DMS' : ISO 6709 degree, minutes, seconds notation i.e. ±DDMMSS.SS±DDDMMSS.SS + 'RT90' : Output format for the Swedish coordinate system RT90 + Return value: a tuple of 2 strings, or a string (for ISO formats) If conversion fails: returns: (None, None) or None (for ISO formats) Some generalities: @@ -292,7 +295,7 @@ def conv_lat_lon(latitude, longitude, format="D.D4"): -180 <= longitude < +180 with +000 prime meridian and -180 180th meridian """ - + # we start the function changing latitude/longitude in english if latitude.find('N') == -1 and latitude.find('S') == -1: # entry is not in english, convert to english @@ -337,12 +340,16 @@ def conv_lat_lon(latitude, longitude, format="D.D4"): str_lon ="-180.0000" return ("%.4f" % lat_float , str_lon) - if format == "D.D8": + if format == "D.D8" or format == "RT90": # correct possible roundoff error str_lon = "%.8f" % (lon_float) if str_lon == "180.00000000": str_lon ="-180.00000000" - return ("%.8f" % lat_float , str_lon) + if format == "RT90": + tx = __conv_WGS84_SWED_RT90(lat_float, lon_float) + return ("%i" % tx[0], "%i" % tx[1]) + else: + return ("%.8f" % lat_float , str_lon) deg_lat = int(lat_float) @@ -481,6 +488,112 @@ def conv_lat_lon(latitude, longitude, format="D.D4"): +def atanh(x): + """arctangent hyperbolicus""" + return 1.0/2.0*math.log((1.0 + x)/(1.0 -x)) + + +def __conv_WGS84_SWED_RT90(lat, lon): + """ + Input is lat and lon as two float numbers + Output is X and Y coordinates in RT90 + as a tuple of float numbers + + The code below converts to/from the Swedish RT90 koordinate + system. The converion functions use "Gauss Conformal Projection + (Transverse Marcator)" Krüger Formulas. + The constanst are for the Swedish RT90-system. + With other constants the conversion should be useful for + other geographical areas. + + """ + # Some constants used for conversion to/from Swedish RT90 + f = 1.0/298.257222101 + e2 = f*(2.0-f) + n = f/(2.0-f) + L0 = math.radians(15.8062845294) # 15 deg 48 min 22.624306 sec + k0 = 1.00000561024 + a = 6378137.0 # meter + at = a/(1.0+n)*(1.0+ 1.0/4.0* pow(n,2)+1.0/64.0*pow(n,4)) + FN = -667.711 # m + FE = 1500064.274 # m + + #the conversion + lat_rad = math.radians(lat) + lon_rad = math.radians(lon) + A = e2 + B = 1.0/6.0*(5.0*pow(e2,2) - pow(e2,3)) + C = 1.0/120.0*(104.0*pow(e2,3) - 45.0*pow(e2,4)) + D = 1.0/1260.0*(1237.0*pow(e2,4)) + DL = lon_rad - L0 + E = A + B*pow(math.sin(lat_rad),2) + \ + C*pow(math.sin(lat_rad),4) + \ + D*pow(math.sin(lat_rad),6) + psi = lat_rad - math.sin(lat_rad)*math.cos(lat_rad)*E + xi = math.atan2(math.tan(psi),math.cos(DL)) + eta = atanh(math.cos(psi)*math.sin(DL)) + B1 = 1.0/2.0*n - 2.0/3.0*pow(n,2) + 5.0/16.0*pow(n,3) + 41.0/180.0*pow(n,4) + B2 = 13.0/48.0*pow(n,2) - 3.0/5.0*pow(n,3) + 557.0/1440.0*pow(n,4) + B3 = 61.0/240.0*pow(n,3) - 103.0/140.0*pow(n,4) + B4 = 49561.0/161280.0*pow(n,4) + X = xi + B1*math.sin(2.0*xi)*math.cosh(2.0*eta) + \ + B2*math.sin(4.0*xi)*math.cosh(4.0*eta) + \ + B3*math.sin(6.0*xi)*math.cosh(6.0*eta) + \ + B4*math.sin(8.0*xi)*math.cosh(8.0*eta) + Y = eta + B1*math.cos(2.0*xi)*math.sinh(2.0*eta) + \ + B2*math.cos(4.0*xi)*math.sinh(4.0*eta) + \ + B3*math.cos(6.0*xi)*math.sinh(6.0*eta) + \ + B4*math.cos(8.0*xi)*math.sinh(8.0*eta) + X = X*k0*at + FN + Y = Y*k0*at + FE + return (X, Y) + +def __conv_SWED_RT90_WGS84(X, Y): + """ + Input is X and Y coordinates in RT90 as float + Output is lat and long in degrees, float as tuple + """ + # Some constants used for conversion to/from Swedish RT90 + f = 1.0/298.257222101 + e2 = f*(2.0-f) + n = f/(2.0-f) + L0 = math.radians(15.8062845294) # 15 deg 48 min 22.624306 sec + k0 = 1.00000561024 + a = 6378137.0 # meter + at = a/(1.0+n)*(1.0+ 1.0/4.0* pow(n,2)+1.0/64.0*pow(n,4)) + FN = -667.711 # m + FE = 1500064.274 # m + + xi = (X - FN)/(k0*at) + eta = (Y - FE)/(k0*at) + D1 = 1.0/2.0*n - 2.0/3.0*pow(n,2) + 37.0/96.0*pow(n,3) - 1.0/360.0*pow(n,4) + D2 = 1.0/48.0*pow(n,2) + 1.0/15.0*pow(n,3) - 437.0/1440.0*pow(n,4) + D3 = 17.0/480.0*pow(n,3) - 37.0/840.0*pow(n,4) + D4 = 4397.0/161280.0*pow(n,4) + xip = xi - D1*math.sin(2.0*xi)*math.cosh(2.0*eta) - \ + D2*math.sin(4.0*xi)*math.cosh(4.0*eta) - \ + D3*math.sin(6.0*xi)*math.cosh(6.0*eta) - \ + D4*math.sin(8.0*xi)*math.cosh(8.0*eta) + etap = eta - D1*math.cos(2.0*xi)*math.sinh(2.0*eta) - \ + D2*math.cos(4.0*xi)*math.sinh(4.0*eta) - \ + D3*math.cos(6.0*xi)*math.sinh(6.0*eta) - \ + D4*math.cos(8.0*xi)*math.sinh(8.0*eta) + psi = math.asin(math.sin(xip)/math.cosh(etap)) + DL = math.atan2(math.sinh(etap),math.cos(xip)) + LON = L0 + DL + A = e2 + pow(e2,2) + pow(e2,3) + pow(e2,4) + B = -1.0/6.0*(7.0*pow(e2,2) + 17*pow(e2,3) + 30*pow(e2,4)) + C = 1.0/120.0*(224.0*pow(e2,3) + 889.0*pow(e2,4)) + D = 1.0/1260.0*(4279.0*pow(e2,4)) + E = A + B*pow(math.sin(psi),2) + \ + C*pow(math.sin(psi),4) + \ + D*pow(math.sin(psi),6) + LAT = psi + math.sin(psi)*math.cos(psi)*E + LAT = math.degrees(LAT) + LON = math.degrees(LON) + return LAT, LON + + #------------------------------------------------------------------------- # # For Testing the convert function in this module, apply it as a script: @@ -497,6 +610,7 @@ if __name__ == '__main__': format4 = "ISO-D" format5 = "ISO-DM" format6 = "ISO-DMS" + format7 = "RT90" print "Testing conv_lat_lon function, "+text+':' res1, res2 = conv_lat_lon(lat1,lon1,format0) print lat1,lon1,"in format",format0, "is ",res1,res2 @@ -511,13 +625,26 @@ if __name__ == '__main__': res = conv_lat_lon(lat1,lon1,format5) print lat1,lon1,"in format",format5, "is",res res = conv_lat_lon(lat1,lon1,format6) - print lat1,lon1,"in format",format6, "is",res,"\n" + print lat1,lon1,"in format",format6, "is",res + res1, res2 = conv_lat_lon(lat1,lon1,format7) + print lat1,lon1,"in format",format7, "is",res1,res2,"\n" def test_formats_fail(lat1,lon1,text=''): print "This test should make conv_lat_lon function fail, "+text+":" res1, res2 = conv_lat_lon(lat1,lon1) print lat1,lon1," fails to convert, result=", res1,res2,"\n" + def test_RT90_conversion(): + """ + a given lat/lon is converted to RT90 and back as a test: + """ + la = 59.0 + 40.0/60. + 9.09/3600.0 + lo = 12.0 + 58.0/60.0 + 57.74/3600.0 + x, y = __conv_WGS84_SWED_RT90(la, lo) + lanew, lonew = __conv_SWED_RT90_WGS84(x,y) + assert math.fabs(lanew - la) < 1e-6, math.fabs(lanew - la) + assert math.fabs(lonew - lo) < 1e-6, math.fabs(lonew - lo) + lat, lon = '50.849888888888', '2.885897222222' test_formats_success(lat,lon) lat, lon = u' 50°50\'59.60"N', u' 2°53\'9.23"E' @@ -636,3 +763,4 @@ if __name__ == '__main__': lat, lon = u'S 50º52\'21.92"', u'W 124º52\'21.92"' test_formats_success(lat,lon, 'New format with S/W first and another º - character') + test_RT90_conversion()