source: trunk/npemap.org.uk/scripts/freethepostcode.org-importer/geo_helper.py @ 33

Last change on this file since 33 was 33, checked in by Nick Burch, 14 years ago

Make a start on the freethepostcode import script

File size: 13.2 KB
Line 
1# Geographical helper functions for nmea_info.py and friends
2#
3# Helps with geographic functions, including:
4#  Lat+Long+Height -> XYZ
5#  XYZ -> Lat+Long+Height
6#  Lat+Long -> other Lat+Long (Helmert Transform)
7#  Lat+Long -> easting/northing (OS Only)
8#  easting/northing -> Lat+Long (OS Only)
9#  OS easting/northing -> OS 6 figure ref
10#
11# See http://gagravarr.org/code/ for updates and information
12#
13# GPL
14#
15# Nick Burch - v0.03 (25/06/2006)
16
17import math
18
19# The earth's radius, in meters, at the equator (should be close enough)
20earths_radius = 6335.437 * 1000
21
22# For each co-ordinate system we do, what are the A, B and E2 values?
23# List is A, B, E^2 (E^2 calculated after)
24abe_values = {
25        'wgs84': [ 6378137.0, 6356752.3141, -1 ],
26        'osgb' : [ 6377563.396, 6356256.91, -1 ],
27        'osie' : [ 6377340.189, 6356034.447, -1 ]
28}
29
30# Calculate the E2 values
31for system in abe_values.keys():
32        a = abe_values[system][0]
33        b = abe_values[system][1]
34        e2 = (a*a - b*b) / (a*a)
35        abe_values[system][2] = e2
36
37# For each co-ordinate system we can translate between, what are
38#  the tx, ty, tz, s, rx, ry and rz values?
39# List is tx, ty, tz, s, rx, ry, rz
40transform_values = {
41        'wgs84_to_osgb' : [ -446.448, 125.157, -542.060, 
42                                                        20.4894 / 1000.0 / 1000.0, # given as ppm
43                                                        -0.1502 / 206265.0, # given as seconds of arc
44                                                        -0.2470 / 206265.0, # given as seconds of arc
45                                                        -0.8421 / 206265.0  # given as seconds of arc
46                                        ],
47        'wgs84_to_osie' : [ -482.530, 130.596, -564.557, 
48                                                        -8.1500 / 1000.0 / 1000.0, # given as ppm
49                                                        -1.0420 / 206265.0, # given as seconds of arc
50                                                        -0.2140 / 206265.0, # given as seconds of arc
51                                                        -0.6310 / 206265.0  # given as seconds of arc
52                                        ],
53        'itrs2000_to_etrs89' : [ 0.054, 0.051, -0.048, 0,
54                                                         0.000081 / 206265.0, # given as seconds of arc
55                                                         0.00049 / 206265.0, # given as seconds of arc
56                                                         0.000792 / 206265.0  # given as seconds of arc
57                                        ]
58}
59
60# Calculate reverse transforms
61for systems in [('wgs84','osgb'), ('wgs84','osie'), ('itrs2000','etrs89')]:
62        fs = systems[0] + "_to_" + systems[1]
63        rs = systems[1] + "_to_" + systems[0]
64        ra = []
65        for val in transform_values[fs]:
66                ra.append(-1.0 * val)
67        transform_values[rs] = ra
68
69# Easting and Northin system values, for the systems we work with.
70# List is n0, e0, F0, theta0 and landa0
71en_values = {
72        'osgb' : [ -100000.0, 400000.0, 0.9996012717,
73                                        49.0 /360.0 *2.0*math.pi,
74                                        -2.0 /360.0 *2.0*math.pi
75                         ],
76        'osie' : [ 250000.0, 200000.0, 1.000035,
77                                        53.5 /360.0 *2.0*math.pi,
78                                        -8.0 /360.0 *2.0*math.pi
79                         ]
80}
81
82##############################################################
83#           OS Specific Helpers for Generic Methods          #
84##############################################################
85
86def turn_wgs84_into_osgb36(lat_dec,long_dec,height):
87        """See http://www.gps.gov.uk/guide6.asp#6.2 and http://www.gps.gov.uk/guide6.asp#6.6 for the calculations, and http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34h.html for some background."""
88
89        wgs84_xyz = turn_llh_into_xyz(lat_dec,long_dec,height,'wgs84')
90
91        osgb_xyz = turn_xyz_into_other_xyz(
92                                        wgs84_xyz[0],wgs84_xyz[1],wgs84_xyz[2],'wgs84','osgb')
93
94        osgb_latlong = turn_xyz_into_llh(
95                                        osgb_xyz[0],osgb_xyz[1],osgb_xyz[2],'osgb')
96
97        return osgb_latlong
98def turn_osgb36_into_wgs84(lat_dec,long_dec,height):
99        """See http://www.gps.gov.uk/guide6.asp#6.2 and http://www.gps.gov.uk/guide6.asp#6.6 for the calculations, and http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34h.html for some background."""
100
101        osgb_xyz = turn_llh_into_xyz(lat_dec,long_dec,height,'osgb')
102
103        wgs84_xyz = turn_xyz_into_other_xyz(
104                                        osgb_xyz[0],osgb_xyz[1],osgb_xyz[2],'osgb','wgs84')
105
106        wgs84_latlong = turn_xyz_into_llh(
107                                        wgs84_xyz[0],wgs84_xyz[1],wgs84_xyz[2],'wgs84')
108
109        return wgs84_latlong
110
111def turn_osgb36_into_eastingnorthing(lat_dec,long_dec):
112        """Turn OSGB36 (decimal) lat/long values into OS easting and northing values."""
113        return turn_latlong_into_eastingnorthing(lat_dec,long_dec,'osgb')
114
115def turn_eastingnorthing_into_osgb36(easting,northing):
116        """Turn OSGB36 easting and northing values into (decimal) lat/long values inOSGB36."""
117        return turn_eastingnorthing_into_latlong(easting,northing,'osgb')
118
119##############################################################
120#             Generic Transform Functions                    #
121##############################################################
122
123def turn_llh_into_xyz(lat_dec,long_dec,height,system):
124        """Convert Lat, Long and Height into 3D Cartesian x,y,z
125       See http://www.ordnancesurvey.co.uk/gps/docs/convertingcoordinates3D.pdf"""
126
127        a = abe_values[system][0]
128        b = abe_values[system][1]
129        e2 = abe_values[system][2]
130
131        theta = float(lat_dec)  / 360.0 * 2.0 * math.pi
132        landa = float(long_dec) / 360.0 * 2.0 * math.pi
133        height = float(height)
134
135        v = a / math.sqrt( 1.0 - e2 * (math.sin(theta) * math.sin(theta)) )
136        x = (v + height) * math.cos(theta) * math.cos(landa)
137        y = (v + height) * math.cos(theta) * math.sin(landa)
138        z = ( (1.0 - e2) * v + height ) * math.sin(theta)
139
140        return [x,y,z]
141
142def turn_xyz_into_llh(x,y,z,system):
143        """Convert 3D Cartesian x,y,z into Lat, Long and Height
144       See http://www.ordnancesurvey.co.uk/gps/docs/convertingcoordinates3D.pdf"""
145
146        a = abe_values[system][0]
147        b = abe_values[system][1]
148        e2 = abe_values[system][2]
149
150        p = math.sqrt(x*x + y*y)
151
152        long = math.atan(y/x)
153        lat_init = math.atan( z / (p * (1.0 - e2)) )
154        v = a / math.sqrt( 1.0 - e2 * (math.sin(lat_init) * math.sin(lat_init)) )
155        lat = math.atan( (z + e2*v*math.sin(lat_init)) / p )
156
157        height = (p / math.cos(lat)) - v # Ignore if a bit out
158
159        # Turn from radians back into degrees
160        long = long / 2 / math.pi * 360
161        lat = lat / 2 / math.pi * 360
162
163        return [lat,long,height]
164
165def turn_xyz_into_other_xyz(old_x,old_y,old_z,from_scheme,to_scheme):
166        """Helmert Transformation between one lat+long system and another
167See http://www.ordnancesurvey.co.uk/oswebsite/gps/information/coordinatesystemsinfo/guidecontents/guide6.html for the calculations, and http://www.movable-type.co.uk/scripts/LatLongConvertCoords.html for a friendlier version with examples"""
168
169        transform = from_scheme + "_to_" + to_scheme
170        tx = transform_values[transform][0]
171        ty = transform_values[transform][1]
172        tz = transform_values[transform][2]
173        s =  transform_values[transform][3]
174        rx = transform_values[transform][4]
175        ry = transform_values[transform][5]
176        rz = transform_values[transform][6]
177
178        # Do the transform
179        new_x = tx + ((1.0+s) * old_x) + (-rz * old_y) + (ry * old_z)
180        new_y = ty + (rz * old_x) + ((1.0+s) * old_y) + (-rx * old_z)
181        new_z = tz + (-ry * old_x) + (rx * old_y) + ((1.0+s) * old_z)
182
183        return [new_x,new_y,new_z]
184
185def calculate_distance_and_bearing(from_lat_dec,from_long_dec,to_lat_dec,to_long_dec):
186        """Uses the spherical law of cosines to calculate the distance and bearing between two positions"""
187
188        # Turn them all into radians
189        from_theta = float(from_lat_dec)  / 360.0 * 2.0 * math.pi
190        from_landa = float(from_long_dec) / 360.0 * 2.0 * math.pi
191        to_theta = float(to_lat_dec)  / 360.0 * 2.0 * math.pi
192        to_landa = float(to_long_dec) / 360.0 * 2.0 * math.pi
193
194        d = math.acos(
195                        math.sin(from_theta) * math.sin(to_theta) +
196                        math.cos(from_theta) * math.cos(to_theta) * math.cos(to_landa-from_landa)
197                ) * earths_radius
198
199        bearing = math.atan2(
200                        math.sin(to_landa-from_landa) * math.cos(to_theta),
201                        math.cos(from_theta) * math.sin(to_theta) -
202                        math.sin(from_theta) * math.cos(to_theta) * math.cos(to_landa-from_landa)
203                )
204        bearing = bearing / 2.0 / math.pi * 360.0
205
206        return [d,bearing]
207
208##############################################################
209#            Easting/Northing Transform Methods              #
210##############################################################
211
212def turn_latlong_into_eastingnorthing(lat_dec,long_dec,scheme):
213        """Turn OSGB36 or OSIE36 (decimal) lat/long values into OS easting and northing values. See http://www.ordnancesurvey.co.uk/oswebsite/gps/information/coordinatesystemsinfo/guidecontents/guide7.html for the calculations, and http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34h.html for some background."""
214
215        n0 = en_values[scheme][0]
216        e0 = en_values[scheme][1]
217        f0 = en_values[scheme][2]
218
219        theta0 = en_values[scheme][3]
220        landa0 = en_values[scheme][4]
221
222        a = abe_values[scheme][0]
223        b = abe_values[scheme][1]
224        e2 = abe_values[scheme][2]
225
226        theta = float(lat_dec)  /360.0 *2.0*math.pi
227        landa = float(long_dec) /360.0 *2.0*math.pi
228
229        n = (a-b) / (a+b)
230        v = a * f0 * math.pow( (1 - e2 * math.sin(theta)*math.sin(theta)), -0.5 )
231        ro = a * f0 * (1 - e2) * math.pow( (1 - e2 * math.sin(theta)*math.sin(theta)), -1.5 )
232        nu2 = v/ro - 1
233
234        M = b * f0 * ( \
235                (1.0 + n + 5.0/4.0 *n*n + 5.0/4.0 *n*n*n) * (theta-theta0) - \
236                (3.0*n + 3.0*n*n + 21.0/8.0 *n*n*n) *math.sin(theta-theta0) *math.cos(theta+theta0) + \
237                (15.0/8.0*n*n + 15.0/8.0*n*n*n) *math.sin(2.0*(theta-theta0)) *math.cos(2.0*(theta+theta0)) - \
238                35.0/24.0*n*n*n *math.sin(3.0*(theta-theta0)) *math.cos(3.0*(theta+theta0)) \
239        )
240
241        I = M + n0
242        II = v/2.0 * math.sin(theta) * math.cos(theta)
243        III = v/24.0 * math.sin(theta) * math.pow( math.cos(theta),3 ) * \
244                (5.0 - math.pow(math.tan(theta),2) + 9.0*nu2)
245        IIIa = v/720.0 * math.sin(theta) * math.pow( math.cos(theta),5 ) * \
246                ( 61.0 - 58.0 *math.pow(math.tan(theta),2) + math.pow(math.tan(theta),4) )
247        IV = v * math.cos(theta)
248        V = v/6.0 * math.pow( math.cos(theta),3 ) * \
249                ( v/ro - math.pow(math.tan(theta),2) )
250        VI = v/120.0 * math.pow(math.cos(theta),5) * \
251                ( 5.0 - 18.0 *math.pow(math.tan(theta),2) + \
252                math.pow(math.tan(theta),4) + 14.0*nu2 - \
253                58.0 * math.pow(math.tan(theta),2)*nu2 )
254
255        northing = I + II*math.pow(landa-landa0,2) + \
256                III*math.pow(landa-landa0,4) + \
257                IIIa*math.pow(landa-landa0,6)
258        easting = e0 + IV*(landa-landa0) + V*math.pow(landa-landa0,3) + \
259                VI*math.pow(landa-landa0,5)
260       
261        return (easting,northing)
262
263def turn_eastingnorthing_into_latlong(easting,northing,scheme):
264        """Turn OSGB36 or OSIE36 easting and northing values into (decimal) lat/long values in OSGB36 / OSIE36. See http://www.ordnancesurvey.co.uk/oswebsite/gps/information/coordinatesystemsinfo/guidecontents/guide7.html for the calculations, and http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34h.html for some background."""
265
266        n0 = en_values[scheme][0]
267        e0 = en_values[scheme][1]
268        f0 = en_values[scheme][2]
269
270        theta0 = en_values[scheme][3]
271        landa0 = en_values[scheme][4]
272
273        a = abe_values[scheme][0]
274        b = abe_values[scheme][1]
275        e2 = abe_values[scheme][2]
276
277        n = (a-b) / (a+b)
278
279        # Prepare to iterate
280        M = 0
281        theta = theta0
282        # Iterate, 4 times should be enough
283        for i in range(4):
284                theta = ((northing - n0 - M) / (a * f0)) + theta
285                M = b * f0 * ( \
286                        (1.0 + n + 5.0/4.0 *n*n + 5.0/4.0 *n*n*n) * (theta-theta0) - \
287                        (3.0*n + 3.0*n*n + 21.0/8.0 *n*n*n) *math.sin(theta-theta0) *math.cos(theta+theta0) + \
288                        (15.0/8.0*n*n + 15.0/8.0*n*n*n) *math.sin(2.0*(theta-theta0)) *math.cos(2.0*(theta+theta0)) - \
289                        35.0/24.0*n*n*n *math.sin(3.0*(theta-theta0)) *math.cos(3.0*(theta+theta0)) \
290                )
291
292        # Compute intermediate values
293        v = a * f0 * math.pow( (1 - e2 * math.sin(theta)*math.sin(theta)), -0.5 )
294        ro = a * f0 * (1 - e2) * math.pow( (1 - e2 * math.sin(theta)*math.sin(theta)), -1.5 )
295        nu2 = v/ro - 1
296        tantheta2 = math.pow(math.tan(theta),2)
297
298        VII = math.tan(theta) / (2 * ro * v)
299        VIII = math.tan(theta) / (24 * ro * math.pow(v,3)) \
300                        * (5 + 3 * tantheta2 + nu2 - \
301                           9 * tantheta2 *  nu2 )
302        IX = math.tan(theta) / (720 * ro * math.pow(v,5)) \
303                        * (61 + 90 * tantheta2 + 45 * tantheta2 * tantheta2)
304        X = 1 / (math.cos(theta) * v)
305        XI = 1 / (math.cos(theta) * 6 * math.pow(v,3)) \
306                        * (v/ro + 2*tantheta2)
307        XII = 1 / (math.cos(theta) * 120 * math.pow(v,5)) \
308                        * (5 + 28 * tantheta2 + 24 * tantheta2 * tantheta2)
309        XIIa = 1 / (math.cos(theta) * 5040 * math.pow(v,7)) \
310                        * (61 + 662 * tantheta2 + 1320 * tantheta2 * tantheta2 \
311                                + 720 * tantheta2 * tantheta2 * tantheta2)
312
313        lat_rad = theta - VII * math.pow((easting-e0),2) \
314                                + VIII * math.pow((easting-e0),4) \
315                                - IX * math.pow((easting-e0),6)
316        long_rad = landa0 + X * (easting-e0) \
317                                - XI * math.pow((easting-e0),3) \
318                                + XII * math.pow((easting-e0),5) \
319                                - XIIa * math.pow((easting-e0),7)
320
321        lat = lat_rad / 2.0 / math.pi * 360.0
322        long = long_rad / 2.0 / math.pi * 360.0
323
324        return (lat,long)
325
326##############################################################
327#             OS Specific Methods Follow                     #
328##############################################################
329
330def turn_easting_northing_into_six_fig(easting,northing):
331        """Turn OS easting and northing values into the six figure OS grid refecence. See http://www.jstott.me.uk/jscoord/"""
332        first_letter = ""
333        second_letter = ""
334
335        easting = int(easting)
336        northing = int(northing)
337
338        # Get the 100 km part
339        hkm_east = int( math.floor(easting / 100000.0) )
340        hkm_north = int( math.floor(northing / 100000.0) )
341        if hkm_north < 5:
342                if hkm_east < 5:
343                        first_letter = "S"
344                else:
345                        first_letter = "T"
346        elif hkm_north < 10:
347                if hkm_east < 5:
348                        first_letter = "N"
349                else:
350                        first_letter = "O"
351        else:
352                first_letter = "H"
353
354        # Get the 10km part
355        index = 65 + ((4 - (hkm_north % 5)) * 5) + (hkm_east % 5)
356        ti = index
357        if index >= 73:
358                index += 1
359        second_letter = chr(index)
360
361        # Get digits 2-4 on easting and northing
362        e = math.floor( (easting  - (100000.0 * hkm_east))  / 100.0)
363        n = math.floor( (northing - (100000.0 * hkm_north)) / 100.0)
364        e = "%03d" % e
365        n = "%03d" % n
366
367        return first_letter + second_letter + e + n
Note: See TracBrowser for help on using the repository browser.