Peter Landgren: 02659 support for RT90 in PlaceUtils.py

svn: r11778

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:
@@ -337,11 +340,15 @@ 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"
+        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)
     
     
@@ -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()