Forum Discussion

samahiji's avatar
samahiji
Frequent Visitor
2 years ago
Solved

Latitude & Longitude calculation from UTM

Hi,

I have a table with utm coordinates that I need to convert them to latitude & longitude in order to be used in arcGIS.

The code below did will for latitude but not for longitude!!

Any idea?

 

let
    Source = ...,
   // Convert UTM to Lat/Long
   AddLatitude = Table.AddColumn(Source, "Latitude", each
       let
           Easting = [BOTTOMX],
           Northing = [BOTTOMY],
           UTMZoneNumber = 39,
           k0 = 0.9996,
           a = 6378137.0,
           e = 0.081819190842622,
           e1sq = 0.006739496742337,
           northernHemisphere = true,
           x = Easting - 500000.0,
           y = if northernHemisphere then Northing else Northing - 10000000.0,
           LongOrigin = (UTMZoneNumber - 1.0) * 6.0 - 180.0 + 3.0,
           M = y / k0,
           mu = M / (a * (1.0 - Number.Power(e, 2.0) / 4.0 - 3.0 * Number.Power(e, 4.0) / 64.0 - 5.0 * Number.Power(e, 6.0) / 256.0)),
           e1 = (1.0 - Number.Sqrt(1.0 - e * e)) / (1.0 + Number.Sqrt(1.0 - e * e)),
           phi1Rad = mu + (3.0 * e1 / 2.0 - 27.0 * Number.Power(e1, 3.0) / 32.0) * Number.Sin(2.0 * mu) + (21.0 * Number.Power(e1, 2.0) / 16.0 - 55.0 * Number.Power(e1, 4.0) / 32.0) * Number.Sin(4.0 * mu) + (151.0 * Number.Power(e1, 3.0) / 96.0) * Number.Sin(6.0 * mu),
           N1 = a / Number.Sqrt(1.0 - Number.Power(e * Number.Sin(phi1Rad), 2.0)),
           T1 = Number.Power(Number.Tan(phi1Rad), 2.0),
           C1 = e1sq * Number.Power(Number.Cos(phi1Rad), 2.0),
           R1 = a * (1.0 - e * e) / Number.Power(1.0 - Number.Power(e * Number.Sin(phi1Rad), 2.0), 1.5),
           D = x / (N1 * k0),
           Latitude = phi1Rad - (N1 * Number.Tan(phi1Rad) / R1) * (Number.Power(D, 2.0) / 2.0 - (5.0 + 3.0 * T1 + 10.0 * C1 - 4.0 * C1 * C1 - 9.0 * e1sq) * Number.Power(D, 4.0) / 24.0 + (61.0 + 90.0 * T1 + 298.0 * C1 + 45.0 * T1 * T1 - 252.0 * e1sq - 3.0 * C1 * C1) * Number.Power(D, 6.0) / 720.0),
           LatitudeDeg = Latitude * (180.0 / Number.PI)
       in
           LatitudeDeg
   ),
   AddLongitude = Table.AddColumn(AddLatitude, "Longitude", each let
   Easting = [BOTTOMX],
   Northing = [BOTTOMY],
   UTMZoneNumber = 39,
   k0 = 0.9996,
   a = 6378137.0,
   e = 0.081819190842622,
   e1sq = 0.006739496742337,
   northernHemisphere = true,
   x = Easting - 500000.0,
   y = if northernHemisphere then Northing else Northing - 10000000.0,
   LongOrigin = (UTMZoneNumber - 1.0) * 6.0 - 180.0 + 3.0,
   M = y / k0,
   mu = M / (a * (1.0 - Number.Power(e, 2.0) / 4.0 - 3.0 * Number.Power(e, 4.0) / 64.0 - 5.0 * Number.Power(e, 6.0) / 256.0)),
   e1 = (1.0 - Number.Sqrt(1.0 - e * e)) / (1.0 + Number.Sqrt(1.0 - e * e)),
   phi1Rad = mu + (3.0 * e1 / 2.0 - 27.0 * Number.Power(e1, 3.0) / 32.0) * Number.Sin(2.0 * mu) + (21.0 * Number.Power(e1, 2.0) / 16.0 - 55.0 * Number.Power(e1, 4.0) / 32.0) * Number.Sin(4.0 * mu) + (151.0 * Number.Power(e1, 3.0) / 96.0) * Number.Sin(6.0 * mu),
   N1 = a / Number.Sqrt(1.0 - Number.Power(e * Number.Sin(phi1Rad), 2.0)),
   T1 = Number.Power(Number.Tan(phi1Rad), 2.0),
   C1 = e1sq * Number.Power(Number.Cos(phi1Rad), 2.0),
   R1 = a * (1.0 - e * e) / Number.Power(1.0 - Number.Power(e * Number.Sin(phi1Rad), 2.0), 1.5),
   D = x / (N1 * k0),
   Longitude = LongOrigin + (D - (1.0 + 2.0 * T1 + C1) * Number.Power(D, 3.0) / 6.0 + (5.0 - 2.0 * C1 + 28.0 * T1 - 3.0 * Number.Power(T1, 2.0) + 8.0 * e1sq + 24.0 * Number.Power(T1, 3.0)) * Number.Power(D, 5.0) / 120.0) / Number.Cos(phi1Rad),
   LongitudeDeg = Longitude * (180.0 / Number.PI)
in
   LongitudeDeg)
in
    AddLongitude