Untitled
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