1 # -*- coding: utf-8 -*-
3 # Adapted from http://positio.rista.net/en/pys60gps/src/KKJWGS84.py
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), \
18 ###########################################################################
19 # Function: KKJ_Zone_I
20 ###########################################################################
22 ZoneNumber = math.floor((KKJI/1000000.0))
23 if ZoneNumber < 0 or ZoneNumber > 5:
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
36 while ZoneNumber >= 0:
37 if math.fabs(KKJlo - KKJ_ZONE_INFO[ZoneNumber][0]) <= 1.5:
39 ZoneNumber = ZoneNumber - 1
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
52 ee = (a * a - bb) / bb
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)
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
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
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.
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)
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
98 MaxLa = MinLa + 0.55 * DeltaLa
99 if (KKJt['I'] < KKJ['I']):
100 MinLo = MinLo + 0.45 * DeltaLo
102 MaxLo = MinLo + 0.55 * DeltaLo
107 ###########################################################################
108 # Function: KKJlalo_to_WGS84lalo
109 ###########################################################################
110 def KKJlalo_to_WGS84lalo(KKJ):
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
126 WGS['La'] = math.degrees(math.radians(KKJ['La']) + dLa)
127 WGS['Lo'] = math.degrees(math.radians(KKJ['Lo']) + dLo)
130 ###########################################################################
131 # Function: WGS84lalo_to_KKJlalo
132 ###########################################################################
133 def WGS84lalo_to_KKJlalo(WGS):
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
144 KKJ['La'] = math.degrees(math.radians(WGS['La']) + dLa)
145 KKJ['Lo'] = math.degrees(math.radians(WGS['Lo']) + dLo)
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)
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)
179 class testCoordinate:
180 def __init__(self, x, y, lon, lat):
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))
211 def testKKJxytoWGS84lalo(x, y):
212 test = { 'P': y, 'I': x }
213 result = KKJxy_to_WGS84lalo(test)
214 return [result['Lo'], result['La']]
218 # Test transforming from KKJxy to WGS84latlon
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:
224 print "Got: (",lon,lat,"), expected: (",t.lon,t.lat,")"
228 print "All tests in testKKJxytoWGS84lalo passed"
231 def testWGS84lalotoKKJxy(lon, lat):
232 test = { 'La': lat, 'Lo': lon }
233 result = WGS84lalo_to_KKJxy(test)
234 return [result['I'], result['P']]
238 # Test transforming from WGS84latlon to KKJxy
240 [x, y] = testWGS84lalotoKKJxy(t.lon, t.lat)
241 if abs(t.x - x) < 2 and abs(t.y - y) < 2:
244 print "Got: (",x,y,"), expected: (",t.x,t.y,")"
248 print "All tests in testWGS84lalotoKKJxy passed"