Skip to content

Instantly share code, notes, and snippets.

@JedIV
Last active September 14, 2018 14:03
Show Gist options
  • Save JedIV/07682f123443676d77a30079a6fcfe2b to your computer and use it in GitHub Desktop.
Save JedIV/07682f123443676d77a30079a6fcfe2b to your computer and use it in GitHub Desktop.
CRS to Lat/Lon with only base packages
import math
def utmToLatLng(utmZone, mean_easting, mean_northing, northernHemisphere):
utmZone = int(utmZone)
mean_easting = float(mean_easting)
mean_northing = float(mean_northing)
northernHemisphere = int(northernHemisphere)
if not northernHemisphere:
mean_northing = 10000000 - mean_northing
a = 6378137
e = 0.081819191
e1sq = 0.006739497
k0 = 0.9996
arc = mean_northing / k0
mu = arc / (a * (1 - math.pow(e, 2) / 4.0 - 3 * math.pow(e, 4) / 64.0 - 5 * math.pow(e, 6) / 256.0))
ei = (1 - math.pow((1 - e * e), (1 / 2.0))) / (1 + math.pow((1 - e * e), (1 / 2.0)))
ca = 3 * ei / 2 - 27 * math.pow(ei, 3) / 32.0
cb = 21 * math.pow(ei, 2) / 16 - 55 * math.pow(ei, 4) / 32
cc = 151 * math.pow(ei, 3) / 96
cd = 1097 * math.pow(ei, 4) / 512
phi1 = mu + ca * math.sin(2 * mu) + cb * math.sin(4 * mu) + cc * math.sin(6 * mu) + cd * math.sin(8 * mu)
n0 = a / math.pow((1 - math.pow((e * math.sin(phi1)), 2)), (1 / 2.0))
r0 = a * (1 - e * e) / math.pow((1 - math.pow((e * math.sin(phi1)), 2)), (3 / 2.0))
fact1 = n0 * math.tan(phi1) / r0
a1 = 500000 - mean_easting
dd0 = a1 / (n0 * k0)
fact2 = dd0 * dd0 / 2
t0 = math.pow(math.tan(phi1), 2)
Q0 = e1sq * math.pow(math.cos(phi1), 2)
fact3 = (5 + 3 * t0 + 10 * Q0 - 4 * Q0 * Q0 - 9 * e1sq) * math.pow(dd0, 4) / 24
fact4 = (61 + 90 * t0 + 298 * Q0 + 45 * t0 * t0 - 252 * e1sq - 3 * Q0 * Q0) * math.pow(dd0, 6) / 720
lof1 = a1 / (n0 * k0)
lof2 = (1 + 2 * t0 + Q0) * math.pow(dd0, 3) / 6.0
lof3 = (5 - 2 * Q0 + 28 * t0 - 3 * math.pow(Q0, 2) + 8 * e1sq + 24 * math.pow(t0, 2)) * math.pow(dd0, 5) / 120
_a2 = (lof1 - lof2 + lof3) / math.cos(phi1)
_a3 = _a2 * 180 / math.pi
latitude = 180 * (phi1 - fact1 * (fact2 + fact3 + fact4)) / math.pi
if not northernHemisphere:
latitude = -latitude
longitude = ((utmZone > 0) and (6 * utmZone - 183.0) or 3.0) - _a3
return (latitude, longitude)
def process(row):
# In 'cell' mode, the process function must return
# a single cell value for each row,
# which will be affected to a new column.
# The 'row' argument is a dictionary of columns of the row
lat_lon = utmToLatLng(row["utmZone"], row["mean_easting"], row["mean_northing"], row["northernHemisphere"])
return lat_lon
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment