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

Last change on this file since 508 was 508, checked in by Nick Burch, 15 years ago

Support IE stuff

File size: 18.6 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 GB+IE Only)
8#  easting/northing -> Lat+Long (OS GB+IE 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.05 (25/03/2007)
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# Cassini Projection Origins
83# List is lat (rad), long (rad), false easting, false northing
84cassini_values = {
85        'osgb' : [ (53.0 + (13.0 / 60.0) + (17.274 / 3600.0)) /360.0 *2.0*math.pi,
86                   -(2.0 + (41.0 / 60.0) +  (3.562 / 3600.0)) /360.0 *2.0*math.pi,
87                   0, 0 ]                               
88}
89
90# How many feet to the meter
91feet_per_meter = 1.0 / 0.3048007491   # 3.28083
92
93##############################################################
94#         OS GB Specific Helpers for Generic Methods         #
95##############################################################
96
97def turn_wgs84_into_osgb36(lat_dec,long_dec,height):
98        """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."""
99
100        wgs84_xyz = turn_llh_into_xyz(lat_dec,long_dec,height,'wgs84')
101
102        osgb_xyz = turn_xyz_into_other_xyz(
103                                        wgs84_xyz[0],wgs84_xyz[1],wgs84_xyz[2],'wgs84','osgb')
104
105        osgb_latlong = turn_xyz_into_llh(
106                                        osgb_xyz[0],osgb_xyz[1],osgb_xyz[2],'osgb')
107        return osgb_latlong
108
109def turn_osgb36_into_wgs84(lat_dec,long_dec,height):
110        """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."""
111
112        osgb_xyz = turn_llh_into_xyz(lat_dec,long_dec,height,'osgb')
113
114        wgs84_xyz = turn_xyz_into_other_xyz(
115                                        osgb_xyz[0],osgb_xyz[1],osgb_xyz[2],'osgb','wgs84')
116
117        wgs84_latlong = turn_xyz_into_llh(
118                                        wgs84_xyz[0],wgs84_xyz[1],wgs84_xyz[2],'wgs84')
119
120        return wgs84_latlong
121
122def turn_osgb36_into_eastingnorthing(lat_dec,long_dec):
123        """Turn OSGB36 (decimal) lat/long values into OS easting and northing values."""
124        return turn_latlong_into_eastingnorthing(lat_dec,long_dec,'osgb')
125
126def turn_eastingnorthing_into_osgb36(easting,northing):
127        """Turn OSGB36 easting and northing values into (decimal) lat/long values inOSGB36."""
128        return turn_eastingnorthing_into_latlong(easting,northing,'osgb')
129
130##############################################################
131#         OS IE Specific Helpers for Generic Methods         #
132##############################################################
133
134def turn_wgs84_into_osie36(lat_dec,long_dec,height):
135        """As per turn_wgs84_into_osgb36, but for Irish grid"""
136
137        wgs84_xyz = turn_llh_into_xyz(lat_dec,long_dec,height,'wgs84')
138
139        osie_xyz = turn_xyz_into_other_xyz(
140                                        wgs84_xyz[0],wgs84_xyz[1],wgs84_xyz[2],'wgs84','osie')
141
142        osie_latlong = turn_xyz_into_llh(
143                                        osie_xyz[0],osie_xyz[1],osie_xyz[2],'osie')
144        return osie_latlong
145
146def turn_osie36_into_wgs84(lat_dec,long_dec,height):
147        """As per turn_osgb36_into_wgs84, but for Irish grid"""
148
149        osie_xyz = turn_llh_into_xyz(lat_dec,long_dec,height,'osie')
150
151        wgs84_xyz = turn_xyz_into_other_xyz(
152                                        osie_xyz[0],osie_xyz[1],osie_xyz[2],'osie','wgs84')
153
154        wgs84_latlong = turn_xyz_into_llh(
155                                        wgs84_xyz[0],wgs84_xyz[1],wgs84_xyz[2],'wgs84')
156
157        return wgs84_latlong
158
159def turn_osie36_into_eastingnorthing(lat_dec,long_dec):
160        """Turn OSIE36 (decimal) lat/long values into OS IE easting and northing values."""
161        return turn_latlong_into_eastingnorthing(lat_dec,long_dec,'osie')
162
163def turn_eastingnorthing_into_osie36(easting,northing):
164        """Turn OSIE36 easting and northing values into (decimal) lat/long values inOSIE36."""
165        return turn_eastingnorthing_into_latlong(easting,northing,'osie')
166
167##############################################################
168#             Generic Transform Functions                    #
169##############################################################
170
171def turn_llh_into_xyz(lat_dec,long_dec,height,system):
172        """Convert Lat, Long and Height into 3D Cartesian x,y,z
173       See http://www.ordnancesurvey.co.uk/gps/docs/convertingcoordinates3D.pdf"""
174
175        a = abe_values[system][0]
176        b = abe_values[system][1]
177        e2 = abe_values[system][2]
178
179        theta = float(lat_dec)  / 360.0 * 2.0 * math.pi
180        landa = float(long_dec) / 360.0 * 2.0 * math.pi
181        height = float(height)
182
183        v = a / math.sqrt( 1.0 - e2 * (math.sin(theta) * math.sin(theta)) )
184        x = (v + height) * math.cos(theta) * math.cos(landa)
185        y = (v + height) * math.cos(theta) * math.sin(landa)
186        z = ( (1.0 - e2) * v + height ) * math.sin(theta)
187
188        return [x,y,z]
189
190def turn_xyz_into_llh(x,y,z,system):
191        """Convert 3D Cartesian x,y,z into Lat, Long and Height
192       See http://www.ordnancesurvey.co.uk/gps/docs/convertingcoordinates3D.pdf"""
193
194        a = abe_values[system][0]
195        b = abe_values[system][1]
196        e2 = abe_values[system][2]
197
198        p = math.sqrt(x*x + y*y)
199
200        long = math.atan(y/x)
201        lat_init = math.atan( z / (p * (1.0 - e2)) )
202        v = a / math.sqrt( 1.0 - e2 * (math.sin(lat_init) * math.sin(lat_init)) )
203        lat = math.atan( (z + e2*v*math.sin(lat_init)) / p )
204
205        height = (p / math.cos(lat)) - v # Ignore if a bit out
206
207        # Turn from radians back into degrees
208        long = long / 2 / math.pi * 360
209        lat = lat / 2 / math.pi * 360
210
211        return [lat,long,height]
212
213def turn_xyz_into_other_xyz(old_x,old_y,old_z,from_scheme,to_scheme):
214        """Helmert Transformation between one lat+long system and another
215See 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"""
216
217        transform = from_scheme + "_to_" + to_scheme
218        tx = transform_values[transform][0]
219        ty = transform_values[transform][1]
220        tz = transform_values[transform][2]
221        s =  transform_values[transform][3]
222        rx = transform_values[transform][4]
223        ry = transform_values[transform][5]
224        rz = transform_values[transform][6]
225
226        # Do the transform
227        new_x = tx + ((1.0+s) * old_x) + (-rz * old_y) + (ry * old_z)
228        new_y = ty + (rz * old_x) + ((1.0+s) * old_y) + (-rx * old_z)
229        new_z = tz + (-ry * old_x) + (rx * old_y) + ((1.0+s) * old_z)
230
231        return [new_x,new_y,new_z]
232
233def calculate_distance_and_bearing(from_lat_dec,from_long_dec,to_lat_dec,to_long_dec):
234        """Uses the spherical law of cosines to calculate the distance and bearing between two positions"""
235
236        # Turn them all into radians
237        from_theta = float(from_lat_dec)  / 360.0 * 2.0 * math.pi
238        from_landa = float(from_long_dec) / 360.0 * 2.0 * math.pi
239        to_theta = float(to_lat_dec)  / 360.0 * 2.0 * math.pi
240        to_landa = float(to_long_dec) / 360.0 * 2.0 * math.pi
241
242        d = math.acos(
243                        math.sin(from_theta) * math.sin(to_theta) +
244                        math.cos(from_theta) * math.cos(to_theta) * math.cos(to_landa-from_landa)
245                ) * earths_radius
246
247        bearing = math.atan2(
248                        math.sin(to_landa-from_landa) * math.cos(to_theta),
249                        math.cos(from_theta) * math.sin(to_theta) -
250                        math.sin(from_theta) * math.cos(to_theta) * math.cos(to_landa-from_landa)
251                )
252        bearing = bearing / 2.0 / math.pi * 360.0
253
254        return [d,bearing]
255
256##############################################################
257#            Easting/Northing Transform Methods              #
258##############################################################
259
260def turn_latlong_into_eastingnorthing(lat_dec,long_dec,scheme):
261        """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."""
262
263        n0 = en_values[scheme][0]
264        e0 = en_values[scheme][1]
265        f0 = en_values[scheme][2]
266
267        theta0 = en_values[scheme][3]
268        landa0 = en_values[scheme][4]
269
270        a = abe_values[scheme][0]
271        b = abe_values[scheme][1]
272        e2 = abe_values[scheme][2]
273
274        theta = float(lat_dec)  /360.0 *2.0*math.pi
275        landa = float(long_dec) /360.0 *2.0*math.pi
276
277        n = (a-b) / (a+b)
278        v = a * f0 * math.pow( (1 - e2 * math.sin(theta)*math.sin(theta)), -0.5 )
279        ro = a * f0 * (1 - e2) * math.pow( (1 - e2 * math.sin(theta)*math.sin(theta)), -1.5 )
280        nu2 = v/ro - 1
281
282        M = b * f0 * ( \
283                (1.0 + n + 5.0/4.0 *n*n + 5.0/4.0 *n*n*n) * (theta-theta0) - \
284                (3.0*n + 3.0*n*n + 21.0/8.0 *n*n*n) *math.sin(theta-theta0) *math.cos(theta+theta0) + \
285                (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)) - \
286                35.0/24.0*n*n*n *math.sin(3.0*(theta-theta0)) *math.cos(3.0*(theta+theta0)) \
287        )
288
289        I = M + n0
290        II = v/2.0 * math.sin(theta) * math.cos(theta)
291        III = v/24.0 * math.sin(theta) * math.pow( math.cos(theta),3 ) * \
292                (5.0 - math.pow(math.tan(theta),2) + 9.0*nu2)
293        IIIa = v/720.0 * math.sin(theta) * math.pow( math.cos(theta),5 ) * \
294                ( 61.0 - 58.0 *math.pow(math.tan(theta),2) + math.pow(math.tan(theta),4) )
295        IV = v * math.cos(theta)
296        V = v/6.0 * math.pow( math.cos(theta),3 ) * \
297                ( v/ro - math.pow(math.tan(theta),2) )
298        VI = v/120.0 * math.pow(math.cos(theta),5) * \
299                ( 5.0 - 18.0 *math.pow(math.tan(theta),2) + \
300                math.pow(math.tan(theta),4) + 14.0*nu2 - \
301                58.0 * math.pow(math.tan(theta),2)*nu2 )
302
303        northing = I + II*math.pow(landa-landa0,2) + \
304                III*math.pow(landa-landa0,4) + \
305                IIIa*math.pow(landa-landa0,6)
306        easting = e0 + IV*(landa-landa0) + V*math.pow(landa-landa0,3) + \
307                VI*math.pow(landa-landa0,5)
308       
309        return (easting,northing)
310
311def turn_eastingnorthing_into_latlong(easting,northing,scheme):
312        """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."""
313
314        n0 = en_values[scheme][0]
315        e0 = en_values[scheme][1]
316        f0 = en_values[scheme][2]
317
318        theta0 = en_values[scheme][3]
319        landa0 = en_values[scheme][4]
320
321        a = abe_values[scheme][0]
322        b = abe_values[scheme][1]
323        e2 = abe_values[scheme][2]
324
325        n = (a-b) / (a+b)
326
327        # Prepare to iterate
328        M = 0
329        theta = theta0
330        # Iterate, 4 times should be enough
331        for i in range(4):
332                theta = ((northing - n0 - M) / (a * f0)) + theta
333                M = b * f0 * ( \
334                        (1.0 + n + 5.0/4.0 *n*n + 5.0/4.0 *n*n*n) * (theta-theta0) - \
335                        (3.0*n + 3.0*n*n + 21.0/8.0 *n*n*n) *math.sin(theta-theta0) *math.cos(theta+theta0) + \
336                        (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)) - \
337                        35.0/24.0*n*n*n *math.sin(3.0*(theta-theta0)) *math.cos(3.0*(theta+theta0)) \
338                )
339
340        # Compute intermediate values
341        v = a * f0 * math.pow( (1 - e2 * math.sin(theta)*math.sin(theta)), -0.5 )
342        ro = a * f0 * (1 - e2) * math.pow( (1 - e2 * math.sin(theta)*math.sin(theta)), -1.5 )
343        nu2 = v/ro - 1
344        tantheta2 = math.pow(math.tan(theta),2)
345
346        VII = math.tan(theta) / (2 * ro * v)
347        VIII = math.tan(theta) / (24 * ro * math.pow(v,3)) \
348                        * (5 + 3 * tantheta2 + nu2 - \
349                           9 * tantheta2 *  nu2 )
350        IX = math.tan(theta) / (720 * ro * math.pow(v,5)) \
351                        * (61 + 90 * tantheta2 + 45 * tantheta2 * tantheta2)
352        X = 1 / (math.cos(theta) * v)
353        XI = 1 / (math.cos(theta) * 6 * math.pow(v,3)) \
354                        * (v/ro + 2*tantheta2)
355        XII = 1 / (math.cos(theta) * 120 * math.pow(v,5)) \
356                        * (5 + 28 * tantheta2 + 24 * tantheta2 * tantheta2)
357        XIIa = 1 / (math.cos(theta) * 5040 * math.pow(v,7)) \
358                        * (61 + 662 * tantheta2 + 1320 * tantheta2 * tantheta2 \
359                                + 720 * tantheta2 * tantheta2 * tantheta2)
360
361        lat_rad = theta - VII * math.pow((easting-e0),2) \
362                                + VIII * math.pow((easting-e0),4) \
363                                - IX * math.pow((easting-e0),6)
364        long_rad = landa0 + X * (easting-e0) \
365                                - XI * math.pow((easting-e0),3) \
366                                + XII * math.pow((easting-e0),5) \
367                                - XIIa * math.pow((easting-e0),7)
368
369        lat = lat_rad / 2.0 / math.pi * 360.0
370        long = long_rad / 2.0 / math.pi * 360.0
371
372        return (lat,long)
373
374##############################################################
375#         Cassini Easting/Northing Transform Methods         #
376##############################################################
377
378def turn_latlong_into_cassini_en(lat_dec,long_dec,scheme):
379        """Latitude and Longitude, into Cassini-Soldner easting and northing co-ordinates, in the given scheme. See http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34g.html for details of the calculation used"""
380       
381        a = abe_values[scheme][0]
382        b = abe_values[scheme][1]
383        e2 = abe_values[scheme][2]
384
385        e4 = e2 * e2
386        e6 = e2 * e2 * e2
387
388        theta = float(lat_dec)  /360.0 *2.0*math.pi
389        landa = float(long_dec) /360.0 *2.0*math.pi
390
391        theta0 = cassini_values[scheme][0]
392        landa0 = cassini_values[scheme][1]
393        false_easting = cassini_values[scheme][2]
394        false_northing = cassini_values[scheme][3]
395
396        # Compute intermediate values
397        A = (landa - landa0) * math.cos(theta)
398        T = math.tan(theta) * math.tan(theta)
399        C = e2 / (1.0 - e2) * math.cos(theta) * math.cos(theta)
400        v = a / math.sqrt( 1 - (e2 * math.sin(theta) * math.sin(theta)) )
401
402        A2 = A ** 2
403        A3 = A ** 3
404        A4 = A ** 4
405        A5 = A ** 5
406
407        # And M, which is how far along the meridian our latitude is from the origin
408        def makeM(picked_theta):
409                return a * (
410                        (1.0 - e2/4.0 - 3.0*e4/64.0 - 5.0*e6/256.0) * picked_theta
411                  - (3.0*e2/8.0 + 3.0*e4/32.0 + 45.0*e6/1024.0) * math.sin(2.0*picked_theta)
412                  + (15.0*e4/256.0 + 45.0*e6/1024.0) * math.sin(4.0*picked_theta)
413                  - (35.0*e6/3072.0) * math.sin(6.0*picked_theta)
414            )
415        M = makeM(theta)
416        M0 = makeM(theta0)
417
418        # Now calculate
419        easting = false_easting + v * (
420                                A - T * A3 / 6.0 - (8.0 - T + 8.0*C) * T * A5 / 120.0 )
421        northing = false_northing + M - M0 + v * math.tan(theta) * (
422                                A2 / 2.0 + (5.0 - T + 6.0*C) * A4 / 24.0 )
423
424        return (easting,northing)
425
426def turn_cassini_en_into_latlong(easting,northing,scheme):
427        """Cassini-Soldner easting and northing, into Latitude and Longitude, in the given scheme. See http://www.posc.org/Epicentre.2_2/DataModel/ExamplesofUsage/eu_cs34g.html for details of the calculation used"""
428       
429        a = abe_values[scheme][0]
430        b = abe_values[scheme][1]
431        e2 = abe_values[scheme][2]
432
433        e4 = e2 * e2
434        e6 = e2 * e2 * e2
435
436        theta0 = cassini_values[scheme][0]
437        landa0 = cassini_values[scheme][1]
438        false_easting = cassini_values[scheme][2]
439        false_northing = cassini_values[scheme][3]
440
441        def makeM(picked_theta):
442                return a * (
443                        (1.0 - e2/4.0 - 3.0*e4/64.0 - 5.0*e6/256.0) * picked_theta
444                  - (3.0*e2/8.0 + 3.0*e4/32.0 + 45.0*e6/1024.0) * math.sin(2.0*picked_theta)
445                  + (15.0*e4/256.0 + 45.0*e6/1024.0) * math.sin(4.0*picked_theta)
446                  - (35.0*e6/3072.0) * math.sin(6.0*picked_theta)
447            )
448
449        # Compute first batch of intermediate values
450        M1 = makeM(theta0) + (northing - false_northing)
451        mu1 = M1 / (a * (1.0 - e2/4.0 - 3.0*e4/64.0 - 5.0*e6/256.0) )
452        e1 = (1 - ((1-e2) ** 0.5)) / (1 + ((1-e2) ** 0.5))
453
454        e1_2 = e1 ** 2
455        e1_3 = e1 ** 3
456        e1_4 = e1 ** 4
457
458        # Now compute theta1 at T1
459        theta1 = mu1 + (
460        + (3.0*e1 / 2.0 - 27.0*e1_3 / 32.0) * math.sin(2.0*mu1)
461        + (21.0*e1_2 / 16.0 - 55.0*e1_4 / 32.0) * math.sin(4.0*mu1)
462        + (151.0*e1_3 / 96.0) * math.sin(6.0*mu1)
463        + (1097.0*e1_4 / 512.0) * math.sin(8.0*mu1) 
464        )
465        T1 = (math.tan(theta1)) ** 2
466
467        # Now we can find v1, ro1 and D
468        v1 = a / math.sqrt( 1.0 - (e2 * math.sin(theta1) * math.sin(theta1)) )
469        ro1 = a * (1 - e2) / ((1 - e2 * math.sin(theta1) * math.sin(theta1)) ** 1.5)
470        D = (easting - false_easting) / v1
471
472        # And finally the lat and long
473        lat = theta1 - (v1 * math.tan(theta1)) / ro1 * (
474                        D*D/2.0 - (1.0 + 3.0 * T1) * ( (D**4) / 24.0 ) )
475        long = landa0 + (
476                                D - T1 * (D**3) / 3.0 + (1 + 3.0 * T1) * T1 * (D**5) / 15.0
477                        ) / math.cos(theta1)
478
479        # Now make decimal versions
480        lat_dec = lat * 360.0 / 2.0 / math.pi
481        long_dec = long * 360.0 / 2.0 / math.pi
482
483        return (lat_dec,long_dec)
484
485##############################################################
486#             OS Specific Methods Follow                     #
487##############################################################
488
489def turn_easting_northing_into_six_fig(easting,northing):
490        """Turn OS easting and northing values into the six figure OS grid refecence. See http://www.jstott.me.uk/jscoord/"""
491        first_letter = ""
492        second_letter = ""
493
494        easting = int(easting)
495        northing = int(northing)
496
497        # Get the 100 km part
498        hkm_east = int( math.floor(easting / 100000.0) )
499        hkm_north = int( math.floor(northing / 100000.0) )
500        if hkm_north < 5:
501                if hkm_east < 5:
502                        first_letter = "S"
503                else:
504                        first_letter = "T"
505        elif hkm_north < 10:
506                if hkm_east < 5:
507                        first_letter = "N"
508                else:
509                        first_letter = "O"
510        else:
511                first_letter = "H"
512
513        # Get the 10km part
514        index = 65 + ((4 - (hkm_north % 5)) * 5) + (hkm_east % 5)
515        ti = index
516        if index >= 73:
517                index += 1
518        second_letter = chr(index)
519
520        # Get digits 2-4 on easting and northing
521        e = math.floor( (easting  - (100000.0 * hkm_east))  / 100.0)
522        n = math.floor( (northing - (100000.0 * hkm_north)) / 100.0)
523        e = "%03d" % e
524        n = "%03d" % n
525
526        return first_letter + second_letter + e + n
Note: See TracBrowser for help on using the repository browser.