77
88# Author: Josh Batchelor <josh.batchelor@ga.gov.au>
99
10-
11- from decimal import *
12- from math import sqrt , log , degrees , sin , cos , sinh , cosh , atan
13- import os
1410import argparse
1511import csv
12+ import os
13+ from decimal import *
14+ from math import atan , cos , cosh , degrees , log , sin , sinh , sqrt
1615
1716# Universal Transverse Mercator Projection Parameters
18- proj = [6378137 , Decimal ('298.25722210088' ), 500000 ,
19- 10000000 , Decimal ('0.9996' ), 6 , - 177 ]
17+ proj = [
18+ 6378137 ,
19+ Decimal ("298.25722210088" ),
20+ 500000 ,
21+ 10000000 ,
22+ Decimal ("0.9996" ),
23+ 6 ,
24+ - 177 ,
25+ ]
2026
2127# Ellipsoidal Constants
2228f = 1 / proj [1 ]
2329semi_maj = proj [0 ]
2430semi_min = float (semi_maj * (1 - f ))
2531ecc1sq = float (f * (2 - f ))
26- ecc2sq = float (ecc1sq / (1 - ecc1sq ))
32+ ecc2sq = float (ecc1sq / (1 - ecc1sq ))
2733ecc1 = sqrt (ecc1sq )
2834n = f / (2 - f )
2935n = float (n )
30- n2 = n ** 2
36+ n2 = n ** 2
3137
3238# Rectifying Radius (Horner Form)
33- A = proj [0 ] / (1 + n ) * ((n2 *
34- (n2 *
35- (n2 *
36- (25 * n2 + 64 )
37- + 256 )
38- + 4096 )
39- + 16384 )
40- / 16384. )
39+ A1 = n2 * (25 * n2 + 64 ) + 256
40+ A2 = n2 * A1 + 4096
41+ A3 = n2 * A2 + 16384
42+ A = proj [0 ] / (1 + n ) * (A3 / 16384.0 )
4143
4244# Beta Coefficients (Horner Form)
43- b2 = ((n *
44- (n *
45- (n *
46- (n *
47- (n *
48- (n *
49- ((37845269 - 31777436 * n ) - 43097152 )
50- + 42865200 )
51- + 752640 )
52- - 104428800 )
53- + 180633600 )
54- - 135475200 ))
55- / 270950400. )
56-
57- b4 = ((n ** 2 *
58- (n *
59- (n *
60- (n *
61- (n *
62- ((- 24749483 * n - 14930208 ) * n + 100683990 )
63- - 152616960 )
64- + 105719040 )
65- - 23224320 )
66- - 7257600 ))
67- / 348364800. )
68-
69- b6 = ((n ** 3 *
70- (n *
71- (n *
72- (n *
73- (n *
74- (232468668 * n - 101880889 )
75- - 39205760 )
76- + 29795040 )
77- + 28131840 )
78- - 22619520 ))
79- / 638668800. )
80-
81- b8 = ((n ** 4 *
82- (n *
83- (n *
84- ((- 324154477 * n - 1433121792 ) * n + 876745056 )
85- + 167270400 )
86- - 208945440 ))
87- / 7664025600. )
88-
89- b10 = ((n ** 5 *
90- (n *
91- (n *
92- (312227409 - 457888660 * n )
93- + 67920528 )
94- - 70779852 ))
95- / 2490808320. )
96-
97- b12 = ((n ** 6 *
98- (n *
99- (19841813847 * n + 3665348512 )
100- - 3758062126 ))
101- / 116237721600. )
102-
103- b14 = ((n ** 7 *
104- (1989295244 * n - 1979471673 ))
105- / 49816166400. )
106-
107- b16 = ((- 191773887257 * n ** 8 ) / 3719607091200. )
45+ b2a = (37845269 - 31777436 * n ) - 43097152
46+ b2b = n * b2a + 42865200
47+ b2c = n * b2b + 752640
48+ b2d = n * b2c - 104428800
49+ b2e = n * b2d + 180633600
50+ b2f = n * b2e - 135475200
51+ b2 = n * b2f / 270950400.0
52+
53+ b4a = (- 24749483 * n - 14930208 ) * n + 100683990
54+ b4b = n * b4a - 152616960
55+ b4c = n * b4b + 105719040
56+ b4d = n * b4c - 23224320
57+ b4e = n * b4d - 7257600
58+ b4 = n ** 2 * b4e / 348364800.0
59+
60+ b6a = 232468668 * n - 101880889
61+ b6b = n * b6a - 39205760
62+ b6c = n * b6b + 29795040
63+ b6d = n * b6c + 28131840
64+ b6e = n * b6d - 22619520
65+ b6 = n ** 3 * b6e / 638668800.0
66+
67+ b8a = (- 324154477 * n - 1433121792 ) * n + 876745056
68+ b8b = n * b8a + 167270400
69+ b8c = n * b8b - 208945440
70+ b8 = n ** 4 * b8c / 7664025600.0
71+
72+ b10a = 312227409 - 457888660 * n
73+ b10b = n * b10a + 67920528
74+ b10c = n * b10b - 70779852
75+ b10 = n ** 5 * b10c / 2490808320.0
76+
77+ b12a = 19841813847 * n + 3665348512
78+ b12b = n * b12a - 3758062126
79+ b12 = n ** 6 * b12b / 116237721600.0
80+
81+ b14a = 1989295244 * n - 1979471673
82+ b14 = n ** 7 * b14a / 49816166400.0
83+
84+ b16 = (- 191773887257 * n ** 8 ) / 3719607091200.0
10885
10986
11087def dd2dms (dd ):
@@ -154,16 +131,23 @@ def grid2geo(zone, easting, northing):
154131 # Finding t using Newtons Method
155132 def sigma (t ):
156133 sigma = sinh (
157- ecc1 * 0.5 * log ((1 + ((ecc1 * t ) / (sqrt (1 + t ** 2 )))) / (1 - ((ecc1 * t ) / (sqrt (1 + t ** 2 ))))))
134+ ecc1
135+ * 0.5
136+ * log (
137+ (1 + ((ecc1 * t ) / (sqrt (1 + t ** 2 ))))
138+ / (1 - ((ecc1 * t ) / (sqrt (1 + t ** 2 ))))
139+ )
140+ )
158141 return sigma
159142
160143 def ftn (t ):
161- ftn = t * sqrt (1 + (sigma (t )) ** 2 ) - sigma (t ) * sqrt (1 + t ** 2 ) - t1
144+ ftn = t * sqrt (1 + (sigma (t )) ** 2 ) - sigma (t ) * sqrt (1 + t ** 2 ) - t1
162145 return ftn
163146
164147 def f1tn (t ):
165- f1tn = (sqrt (1 + (sigma (t )) ** 2 ) * sqrt (1 + t ** 2 ) - sigma (t ) * t ) * (
166- ((1 - float (ecc1sq )) * sqrt (1 + t ** 2 )) / (1 + (1 - float (ecc1sq )) * t ** 2 ))
148+ f1tn = (sqrt (1 + (sigma (t )) ** 2 ) * sqrt (1 + t ** 2 ) - sigma (t ) * t ) * (
149+ ((1 - float (ecc1sq )) * sqrt (1 + t ** 2 )) / (1 + (1 - float (ecc1sq )) * t ** 2 )
150+ )
167151 return f1tn
168152
169153 t2 = t1 - (ftn (t1 )) / (f1tn (t1 ))
@@ -196,9 +180,9 @@ def grid2geoio(fn):
196180 csvfile = open (fn )
197181 csvreader = csv .reader (csvfile )
198182 # Create Output File
199- fn_part = ( os .path .splitext (fn ) )
200- fn_out = fn_part [0 ] + ' _out' + fn_part [1 ]
201- outfile = open (fn_out , 'w' , newline = '' )
183+ fn_part = os .path .splitext (fn )
184+ fn_out = fn_part [0 ] + " _out" + fn_part [1 ]
185+ outfile = open (fn_out , "w" , newline = "" )
202186 # Write Output
203187 outfilewriter = csv .writer (outfile )
204188 # outfilewriter.writerow(['Pt', 'Latitude', 'Longitude'])
@@ -216,15 +200,17 @@ def grid2geoio(fn):
216200 # Close Files
217201 outfile .close ()
218202 csvfile .close ()
219- return ' Output saved as ' + str (fn_out )
203+ return " Output saved as " + str (fn_out )
220204
221205
222206if __name__ == "__main__" :
223- parser = argparse .ArgumentParser (description = 'Batch Converter of Map Grid of Australia grid co-ordinates to '
224- 'Geodetic Datum of Australia geographic co-ordinates. Files must '
225- 'be .csv and of the format Pt ID, Zone, Easting, Northing with no '
226- 'header line.' )
227- parser .add_argument ('f' , metavar = '--file' , type = str , help = 'Input Filename' )
207+ parser = argparse .ArgumentParser (
208+ description = "Batch Converter of Map Grid of Australia grid co-ordinates to "
209+ "Geodetic Datum of Australia geographic co-ordinates. Files must "
210+ "be .csv and of the format Pt ID, Zone, Easting, Northing with no "
211+ "header line."
212+ )
213+ parser .add_argument ("f" , metavar = "--file" , type = str , help = "Input Filename" )
228214 args = parser .parse_args ()
229215 fn = args .f
230216
0 commit comments