# Untitled

unknown
python
10 months ago
3.0 kB
3
Indexable
Never
```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
```