-
Notifications
You must be signed in to change notification settings - Fork 0
/
geodetic.py
64 lines (53 loc) · 2.06 KB
/
geodetic.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
import math, time
import xml.etree.ElementTree as ET
a = 6378137
f = 1/298.257222101
e = math.sqrt((2*f)-(f**2))
class geodetic:
@staticmethod
def geo(lat, long):
global a, f, e
tree = ET.parse('ui.xml')
root = tree.getroot()
zone = list(tree.iter('zone'))[0].text
if zone == '3601':
lambdacm = 120.5
LatO = 43.0 + 40.0 / 60.0
LatS = 44.0 + 20.0 / 60.0
LatN = 46.0
Eo = 2500000 #meters
Nb = 0
elif zone == '3602':
lambdacm = 120.5
LatO = 41.0 + 40.0 / 60.0
LatS = 42.0 + 20.0 / 60.0
LatN = 44.0
Eo = 1500000 #meters
Nb = 0
else:
print "ERROR: Zone number not recognized"
lat = math.radians(lat)
long = math.radians(long)
LatO = math.radians(LatO)
LatS = math.radians(LatS)
LatN = math.radians(LatN)
W = math.sqrt(1-(e**2)*((math.sin(lat))**2))
M = math.cos(lat)/W
T = math.sqrt(((1-math.sin(lat))/(1+math.sin(lat)))*(((1+e*math.sin(lat))/(1-e*math.sin(lat)))**e))
w1 = math.sqrt(1-(e**2)*((math.sin(LatS))**2))
w2 = math.sqrt(1-(e**2)*((math.sin(LatN))**2))
m1 = math.cos(LatS)/w1
m2 = math.cos(LatN)/w2
t0 = math.sqrt(((1-math.sin(LatO))/(1+math.sin(LatO)))*((1+e*math.sin(LatO))/(1-e*math.sin(LatO)))**e)
t1 = math.sqrt(((1-math.sin(LatS))/(1+math.sin(LatS)))*(((1+e*math.sin(LatS)))/(1-e*math.sin(LatS)))**e)
t2 = math.sqrt(((1-math.sin(LatN))/(1+math.sin(LatN)))*(((1+e*math.sin(LatN)))/(1-e*math.sin(LatN)))**e)
n = (math.log(m1)-math.log(m2))/((math.log(t1)-math.log(t2)))
F = m1/(n*t1**n)
Rb = a*F*t0**n
gamma = (lambdacm - math.degrees(long)) * n #convergency angle in decimal degrees
R = a*F*T**n
k = (R*n) / (a*M) #grid_factor
gamma = math.radians(gamma)
easting = (R * math.sin(gamma)) + Eo
northing = Rb - (R * math.cos(gamma)) + Nb
return(northing, easting, k, gamma)