Untitled
unknown
python
2 years ago
3.0 kB
11
Indexable
import math
def lat_long_to_oig(lat, long):
olat = lat * math.pi / 180
olon = long * math.pi / 180
from_esq = 0.0066943800042608068
a = 6378137.0
e = math.sqrt(from_esq)
b = 6356752.3141
slat1 = math.sin(olat)
clat1 = math.cos(olat)
clat1sq = math.cos(olat) ** 2
tanlat1sq = slat1 * slat1 / clat1sq
e2 = from_esq
e4 = from_esq ** 2
e6 = from_esq ** 3
eg2 = (e * a / b) ** 2
l1 = 1 - e2 / 4 - 3 * e4 / 64 - 5 * e6 / 256
l2 = 3 * e2 / 8 + 3 * e4 / 32 + 45 * e6 / 1024
l3 = 15 * e4 / 256 + 45 * e6 / 1024
l4 = 35 * e6 / 3072
M = a * (l1 * olat - l2 * math.sin(2 * olat) + l3 * math.sin(4 * olat) - l4 * math.sin(6 * olat))
nu = a / math.sqrt(1 - (e * slat1) * (e * slat1))
p = olon - 0.61443473225468925
k0 = 1.0000067
K1 = M * k0
K2 = k0 * nu * slat1 * clat1 / 2
K3 = (k0 * nu * slat1 * clat1 * clat1sq / 24) * (5 - tanlat1sq + 9 * eg2 * clat1sq + 4 * eg2 * eg2 * clat1sq * clat1sq)
K4 = k0 * nu * clat1
K5 = (k0 * nu * clat1 * clat1sq / 6) * (1 - tanlat1sq + eg2 * clat1 * clat1)
X = round(K4 * p + K5 * p * p * p + 219529.584, 0) - 50103
Y = (round(K1 + K2 * p * p + K3 * p * p * p * p - 2885516.9488, 0) - 216 + 500000) % 1000000
return X, Y
def utm_to_lat_long(x, y):
E = x / 1000
N = y / 1000
a = 6378.137
f = 1 / 298.257223563
N0 = 0
k0 = 0.9996
E0 = 500
n = f / (2 - f)
A = (a / (n + 1)) * (1 + (n ** 2) / 4 + (n ** 4) / 64) # ...
b1 = n / 2 - 2 * (n ** 2) / 3 + 37 * (n ** 3) / 96
b2 = (n ** 2) / 48 + (n ** 3) / 15
b3 = 17 * (n ** 3) / 480
d1 = 2 * n - 2 * (n ** 2) / 3 - 2 * (n ** 3)
d2 = 7 * (n ** 2) / 3 - 8 * (n ** 3) / 5
d3 = 56 * (n ** 3) / 15
zone = 36
ep = (N - N0) / (k0 * A)
et = (E - E0) / (k0 * A)
ept = ep - (b1 * math.sin(2 * 1 * ep) * math.cosh(2 * 1 * et) +
b2 * math.sin(2 * 2 * ep) * math.cosh(2 * 2 * et) +
b3 * math.sin(2 * 3 * ep) * math.cosh(2 * 3 * et))
ett = et - (b1 * math.cos(2 * 1 * ep) * math.sinh(2 * 1 * et) +
b2 * math.cos(2 * 2 * ep) * math.sinh(2 * 2 * et) +
b3 * math.cos(2 * 3 * ep) * math.sinh(2 * 3 * et))
X = math.asin(math.sin(ept) / math.cosh(ett))
lat = (X + (d1 * math.sin(2 * 1 * X) +
d2 * math.sin(2 * 2 * X) +
d3 * math.sin(2 * 3 * X))) * 180 / math.pi
l0 = zone * 6 - 183
long = l0 + math.atan(math.sinh(ett) / math.cos(ept)) * 180 / math.pi
return lat, long
def utm_to_oig(utm_x, utm_y):
utm_y = utm_y % 1000000
if str(utm_y)[0] != "3":
utm_y += 3000000
lat, long = utm_to_lat_long(utm_x, utm_y)
return lat_long_to_oig(lat, long)
def nig_to_oig(oig_x, oig_y):
return oig_x - 50000, (oig_y + 500000) % 1000000
def oig_to_nig(nig_x, nig_y):
return nig_x + 50000, (nig_y + 500000) % 1000000
Editor is loading...
Leave a Comment