Untitled

mail@pastecode.io avatar
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
Leave a Comment