vcs.maemo.org Git - ptas/blob - tests/coord_trans_proto.py

   1 # -*- coding: utf-8 -*-
   2
   3 # Adapted from http://positio.rista.net/en/pys60gps/src/KKJWGS84.py
   4
   5 import math
   6
   7
   8 # Constants
   9 # Longitude0 and Center meridian of KKJ bands
  10 KKJ_ZONE_INFO = { 0: (18.0,  500000.0), \
  11                   1: (21.0, 1500000.0), \
  12                   2: (24.0, 2500000.0), \
  13                   3: (27.0, 3500000.0), \
  14                   4: (30.0, 4500000.0), \
  15                   5: (33.0, 5500000.0), \
  16                 }
  17
  18 ###########################################################################
  19 # Function:  KKJ_Zone_I
  20 ###########################################################################
  21 def KKJ_Zone_I(KKJI):
  22     ZoneNumber = math.floor((KKJI/1000000.0))
  23     if ZoneNumber < 0 or ZoneNumber > 5:
  24         ZoneNumber = -1
  25     return ZoneNumber
  26
  27 ###########################################################################
  28 # Function:  KKJ_Zone_Lo
  29 ###########################################################################
  30 def KKJ_Zone_Lo(KKJlo):
  31     # determine the zonenumber from KKJ easting
  32     # takes KKJ zone which has center meridian
  33     # longitude nearest (in math value) to
  34     # the given KKJ longitude
  35     ZoneNumber = 5
  36     while ZoneNumber >= 0:
  37         if math.fabs(KKJlo - KKJ_ZONE_INFO[ZoneNumber][0]) <= 1.5:
  38             break
  39         ZoneNumber = ZoneNumber - 1
  40     return ZoneNumber
  41
  42 ###########################################################################
  43 # Function:  KKJlalo_to_KKJxy
  44 ###########################################################################
  45 def KKJlalo_to_KKJxy(INP, ZoneNumber):
  46     Lo = math.radians(INP['Lo']) - math.radians(KKJ_ZONE_INFO[ZoneNumber][0])
  47     a  = 6378388.0            # Hayford ellipsoid
  48     f  = 1/297.0
  49     b  = (1.0 - f) * a
  50     bb = b * b
  51     c  = (a / b) * a
  52     ee = (a * a - bb) / bb
  53     n = (a - b)/(a + b)
  54     nn = n * n
  55     cosLa = math.cos(math.radians(INP['La']))
  56     NN = ee * cosLa * cosLa
  57     LaF = math.atan(math.tan(math.radians(INP['La'])) / math.cos(Lo * math.sqrt(1 + NN)))
  58     cosLaF = math.cos(LaF)
  59     t   = (math.tan(Lo) * cosLaF) / math.sqrt(1 + ee * cosLaF * cosLaF)
  60     A   = a / ( 1 + n )
  61     A1  = A * (1 + nn / 4 + nn * nn / 64)
  62     A2  = A * 1.5 * n * (1 - nn / 8)
  63     A3  = A * 0.9375 * nn * (1 - nn / 4)
  64     A4  = A * 35/48.0 * nn * n
  65     OUT = {}
  66     OUT['P'] = A1 * LaF - \
  67         A2 * math.sin(2 * LaF) + \
  68             A3 * math.sin(4 * LaF) - \
  69                 A4 * math.sin(6 * LaF)
  70     OUT['I'] = c * math.log(t + math.sqrt(1+t*t)) + \
  71         500000.0 + ZoneNumber * 1000000.0
  72     return OUT
  73
  74 ###########################################################################
  75 # Function:  KKJxy_to_KKJlalo
  76 ###########################################################################
  77 def KKJxy_to_KKJlalo(KKJ):
  78     # Scan iteratively the target area, until find matching
  79     # KKJ coordinate value.  Area is defined with Hayford Ellipsoid.
  80
  81     LALO = {}
  82     ZoneNumber = KKJ_Zone_I(KKJ['I'])
  83     MinLa = math.radians(59.0)
  84     MaxLa = math.radians(70.5)
  85     MinLo = math.radians(18.5)
  86     MaxLo = math.radians(32.0)
  87
  88     i = 1
  89     while (i < 35):
  90         DeltaLa = MaxLa - MinLa
  91         DeltaLo = MaxLo - MinLo
  92         LALO['La'] = math.degrees(MinLa + 0.5 * DeltaLa)
  93         LALO['Lo'] = math.degrees(MinLo + 0.5 * DeltaLo)
  94         KKJt = KKJlalo_to_KKJxy(LALO, ZoneNumber)
  95         if (KKJt['P'] < KKJ['P']):
  96             MinLa = MinLa + 0.45 * DeltaLa
  97         else:
  98             MaxLa = MinLa + 0.55 * DeltaLa
  99         if (KKJt['I'] < KKJ['I']):
 100             MinLo = MinLo + 0.45 * DeltaLo
 101         else:
 102             MaxLo = MinLo + 0.55 * DeltaLo
 103         i = i + 1
 104
 105     return LALO
 106
 107 ###########################################################################
 108 # Function:  KKJlalo_to_WGS84lalo
 109 ###########################################################################
 110 def KKJlalo_to_WGS84lalo(KKJ):
 111     La = KKJ['La']
 112     Lo = KKJ['Lo']
 113     dLa = math.radians( 0.124867E+01      + \
 114                         -0.269982E+00 * La + \
 115                         0.191330E+00 * Lo + \
 116                         0.356119E-02 * La * La + \
 117                         -0.122312E-02 * La * Lo + \
 118                         -0.335514E-03 * Lo * Lo ) / 3600.0
 119     dLo = math.radians(-0.286111E+02      + \
 120                         0.114183E+01 * La + \
 121                         -0.581428E+00 * Lo + \
 122                         -0.152421E-01 * La * La + \
 123                         0.118177E-01 * La * Lo + \
 124                         0.826646E-03 * Lo * Lo ) / 3600.0
 125     WGS = {}
 126     WGS['La'] = math.degrees(math.radians(KKJ['La']) + dLa)
 127     WGS['Lo'] = math.degrees(math.radians(KKJ['Lo']) + dLo)
 128     return WGS
 129
 130 ###########################################################################
 131 # Function:  WGS84lalo_to_KKJlalo
 132 ###########################################################################
 133 def WGS84lalo_to_KKJlalo(WGS):
 134     La = WGS['La']
 135     Lo = WGS['Lo']
 136     dLa = math.radians(-0.124766E+01      +           0.269941E+00 * La +                          -0.191342E+00 * Lo +  -0.356086E-02 * La * La +                          0.122353E-02 * La * Lo +                           0.335456E-03 * Lo * Lo ) / 3600.0
 137     dLo = math.radians( 0.286008E+02      + \
 138                         -0.114139E+01 * La + \
 139                         0.581329E+00 * Lo + \
 140                         0.152376E-01 * La * La + \
 141                         -0.118166E-01 * La * Lo + \
 142                         -0.826201E-03 * Lo * Lo ) / 3600.0
 143     KKJ = {}
 144     KKJ['La'] = math.degrees(math.radians(WGS['La']) + dLa)
 145     KKJ['Lo'] = math.degrees(math.radians(WGS['Lo']) + dLo)
 146     return KKJ
 147
 148 ###########################################################################
 149 # Function:  KKJxy_to_WGS84lalo
 150 ###########################################################################
 151 # Input:     dictionary with ['P'] is KKJ Northing
 152 #                            ['I'] in KKJ Eeasting
 153 # Output:    dictionary with ['La'] is latitude in degrees (WGS84)
 154 #                            ['Lo'] is longitude in degrees (WGS84)
 155 ###########################################################################
 156 def KKJxy_to_WGS84lalo(KKJin):
 157     KKJz = KKJxy_to_KKJlalo(KKJin)
 158     WGS = KKJlalo_to_WGS84lalo(KKJz)
 159     return WGS
 160
 161 ###########################################################################
 162 # Function:  WGS84lalo_to_KKJxy
 163 ###########################################################################
 164 # Input:     dictionary with ['La'] is latitude in degrees (WGS84)
 165 #                            ['Lo'] is longitude in degrees (WGS84)
 166 # Output:    dictionary with ['P'] is KKJ Northing
 167 #                            ['I'] in KKJ Eeasting
 168 ###########################################################################
 169 def WGS84lalo_to_KKJxy(WGSin):
 170     KKJlalo = WGS84lalo_to_KKJlalo(WGSin)
 171     ZoneNumber = KKJ_Zone_Lo(KKJlalo['Lo'])
 172     KKJxy = KKJlalo_to_KKJxy(KKJlalo, ZoneNumber)
 173     return KKJxy
 174
 175 ###########
 176 # Test code
 177 ###########
 178
 179 class testCoordinate:
 180     def __init__(self, x, y, lon, lat):
 181         self.x = x
 182         self.y = y
 183         self.lon = lon
 184         self.lat = lat
 185
 186 testData = []
 187 # Test data extracted from example on page
 188 # http://developer.reittiopas.fi/pages/fi/http-get-interface.php
 189 testData.append(testCoordinate(2556686, 6682815, 25.02051, 60.2528))
 190 testData.append(testCoordinate(2546340, 6675352, 24.832, 60.18713))
 191 testData.append(testCoordinate(2557985, 6685213, 25.04465, 60.27414))
 192 testData.append(testCoordinate(2556532, 6682578, 25.01767, 60.2507))
 193 testData.append(testCoordinate(2524959, 6686629, 24.44804, 60.2902))
 194 testData.append(testCoordinate(2559094, 6693721, 25.06718, 60.35033))
 195 testData.append(testCoordinate(2556861, 6683030, 25.02373, 60.25471))
 196 testData.append(testCoordinate(2556888, 6682971, 25.0242, 60.25417))
 197 testData.append(testCoordinate(2560257, 6698983, 25.08981, 60.39737))
 198 testData.append(testCoordinate(2562518, 6686969, 25.12709, 60.28923))
 199 testData.append(testCoordinate(2536615, 6673635, 24.65643, 60.1727))
 200 testData.append(testCoordinate(2559118, 6693833, 25.06764, 60.35133))
 201 testData.append(testCoordinate(2559182, 6693629, 25.06874, 60.34949))
 202 testData.append(testCoordinate(2556947, 6682640, 25.02518, 60.25119))
 203 testData.append(testCoordinate(2556822, 6682723, 25.02294, 60.25196))
 204 testData.append(testCoordinate(2559089, 6693605, 25.06705, 60.34929))
 205 testData.append(testCoordinate(2546445, 6675512, 24.83393, 60.18855))
 206 testData.append(testCoordinate(2556964, 6682609, 25.02547, 60.25091))
 207 testData.append(testCoordinate(2556740, 6682861, 25.0215, 60.25321))
 208 testData.append(testCoordinate(2559002, 6694007, 25.06559, 60.35291))
 209
 210
 211 def testKKJxytoWGS84lalo(x, y):
 212     test = { 'P': y, 'I': x }
 213     result = KKJxy_to_WGS84lalo(test)
 214     return [result['Lo'], result['La']]
 215
 216 testsPass = True
 217
 218 # Test transforming from KKJxy to WGS84latlon
 219 for t in testData:
 220     [lon, lat] = testKKJxytoWGS84lalo(t.x, t.y)
 221     if math.fabs(t.lon - lon) < 0.001 and math.fabs(t.lat - lat) < 0.001:
 222         pass
 223     else:
 224         print "Got: (",lon,lat,"), expected: (",t.lon,t.lat,")"
 225         testsPass = False
 226
 227 if testsPass:
 228     print "All tests in testKKJxytoWGS84lalo passed"
 229
 230
 231 def testWGS84lalotoKKJxy(lon, lat):
 232     test = { 'La': lat, 'Lo': lon }
 233     result = WGS84lalo_to_KKJxy(test)
 234     return [result['I'], result['P']]
 235
 236 testsPass = True
 237
 238 # Test transforming from WGS84latlon to KKJxy
 239 for t in testData:
 240     [x, y] = testWGS84lalotoKKJxy(t.lon, t.lat)
 241     if abs(t.x - x) < 2 and abs(t.y - y) < 2:
 242         pass
 243     else:
 244         print "Got: (",x,y,"), expected: (",t.x,t.y,")"
 245         testsPass = False
 246
 247 if testsPass:
 248     print "All tests in testWGS84lalotoKKJxy passed"