Covert point geometry to Latitude and Longitude

[RyanBohan]

Increase functionality for Geometry($feature) to covert point geometry to Latitude and Longitude

This is possible with a bit of arcade, would be great to have a ToLatitude(geometry($feature).x) and ToLongitude(geometry($feature).y)

Arcade Functionality from XanderBakker

var x = Geometry($feature).x

var y = Geometry($feature).y

// Converts XY point from Spherical Mercator EPSG:900913 to lat/lon in WGS84 Datum

// Fuente: http://www.maptiler.org/google-maps-coordinates-tile-bounds-projection/

var originShift = 2.0 * PI * 6378137.0 / 2.0;

var lon = (x / originShift) * 180.0;

var lat = (y / originShift) * 180.0;

lat = 180.0 / PI * (2.0 * Atan( Exp( lat * PI / 180.0)) - PI / 2.0);

return [lat, lon];

[AbiDhakal]

Yes, from any projection to Lat/lon.

[CyrusCaughey]

I created a Transverse Mercator version of the formulas from this post. If parameters are changed appropriately for each projection, this should hopefully work for all Transverse Mercator projections. I was working in State Plane NAD83(HARN)/New Mexico Central, so my units were in feet, but the formulas are made for meters.

//Call feature geometry coordinates in feet and convert to meters
var easting = Geometry($feature).x * 0.3048
var northing = Geometry($feature).y * 0.3048

//Define supplementary hyperbolic trig functions
function Asinh(x) {
    return log(x + Sqrt(1 + Pow(x, 2)))
}

function Atanh(x) {
      return (log(1+x) - log(1-x))/2
}

function Sinh (x) {
    return (Exp(x) - Exp((-1*x)))/2
}

function Cosh (x) {
    return (Exp(x) + Exp((-1*x)))/2
}

function Tanh (x) {
    return (Exp(x) - Exp((-1*x)))/(Exp(x) + Exp((-1*x)))
}

//Function to convert Transverse Mercator Coordinates to WGS84 LatLon, derived from formulas in https://www.iogp.org/wp-content/uploads/2019/09/373-07-02.pdf starting on page 52
//Find parameters for each projection at epsg.io
//From the OGC Well Known Text definition at the bottom of the page, the semimajor axis and flattening values will be found in  SPHEROID[crs,semimajor_axis,flattening]
function convertTransverseMercator(x, y, semimajor_axis, flattening, latitude_of_origin, central_meridian, scale_factor, false_easting, false_northing) {
  var a = semimajor_axis
  var f = 1 / flattening
  var b = a*(1 - f)
  var e = Sqrt(1 - Pow((b/a), 2))
  var phi_0 = latitude_of_origin
  var phi_0_rad = phi_0 * (PI/180)
  var lambda_0 = central_meridian
  var lambda_0_rad = lambda_0 * (PI/180)
  var k_0 = scale_factor
  var FE = false_easting
  var FN = false_northing
  var n = f/(2-f)
  var B_ = (a/(1+n)) * (1+(Pow(n, 2)/4) + (Pow(n, 4)/64))

  var h_1 = (n/2) - ((2/3)*Pow(n, 2)) + ((5/16)*Pow(n, 3)) + ((41/180)*Pow(n, 4))
  var h_2 = ((13/48)*Pow(n, 2)) - ((3/5)*Pow(n, 3)) + ((557/1440)*Pow(n, 4))
  var h_3 = ((61/240)*Pow(n, 3)) - ((103/140)*Pow(n, 4))
  var h_4 = (49561/161280)*Pow(n, 4)
  
  var Q_0 = Asinh(Tan(phi_0_rad)) - (e * Atanh((e * Sin(phi_0_rad))))
  var epsilon_0_0 = Atan(Sinh(Q_0))
  var epsilon_0_1 = h_1 * Sin(2*epsilon_0_0)
  var epsilon_0_2 = h_2 * Sin(4*epsilon_0_0)
  var epsilon_0_3 = h_3 * Sin(6*epsilon_0_0)
  var epsilon_0_4 = h_4 * Sin(8*epsilon_0_0)
  var epsilon_0 = epsilon_0_0 + epsilon_0_1 + epsilon_0_2 + epsilon_0_3 + epsilon_0_4
  var M_0 = B_ * epsilon_0

  var h_1_rev = (n/2) - ((2/3)*Pow(n, 2)) + ((37/96)*Pow(n, 3)) - ((1/360)*Pow(n, 4))
  var h_2_rev = ((1/48)*Pow(n, 2)) + ((1/15)*Pow(n, 3)) - ((437/1440)*Pow(n, 4))
  var h_3_rev = ((17/480)*Pow(n, 3)) - ((37/840)*Pow(n, 4))
  var h_4_rev = (4397/161280) * Pow(n, 4)

  var eta_rev = (x - FE)/(B_*k_0)
  var epsilon_rev = ((y - FN) + (k_0 * M_0)) / (B_ * k_0)

  var epsilon_1_rev = h_1_rev * Sin(2 * epsilon_rev) * Cosh(2 * eta_rev)
  var epsilon_2_rev = h_2_rev * Sin(4 * epsilon_rev) * Cosh(4 * eta_rev)
  var epsilon_3_rev = h_3_rev * Sin(6 * epsilon_rev) * Cosh(6 * eta_rev)
  var epsilon_4_rev = h_4_rev * Sin(8 * epsilon_rev) * Cosh(8 * eta_rev)
  var epsilon_0_rev = epsilon_rev - (epsilon_1_rev + epsilon_2_rev + epsilon_3_rev + epsilon_4_rev)

  var eta_1_rev = h_1_rev * Cos(2 * epsilon_rev) * Sinh(2 * eta_rev)
  var eta_2_rev = h_2_rev * Cos(4 * epsilon_rev) * Sinh(4 * eta_rev)
  var eta_3_rev = h_3_rev * Cos(6 * epsilon_rev) * Sinh(6 * eta_rev)
  var eta_4_rev = h_4_rev * Cos(8 * epsilon_rev) * Sinh(8 * eta_rev)
  var eta_0_rev = eta_rev - (eta_1_rev + eta_2_rev + eta_3_rev + eta_4_rev)

  var beta_rev = Asin(Sin(epsilon_0_rev)/Cosh(eta_0_rev))
  var Q_rev = Asinh(Tan(beta_rev))
  var Q_rev_1 = Q_rev + (e * Atanh((e * Tanh(Q_rev))))
  var Q_rev_2 = Q_rev + (e * Atanh((e * Tanh(Q_rev_1))))
  var Q_rev_3 = Q_rev + (e * Atanh((e * Tanh(Q_rev_2))))
  var Q_rev_4 = Q_rev + (e * Atanh((e * Tanh(Q_rev_3))))
  //Changes in Q often became insignificant after the 4th iteration, but you are welcome to iterate further

  var phi_rad = Atan(Sinh(Q_rev_4))
  var phi = phi_rad * (180/PI)
  var lambda_rad = lambda_0_rad + Asin(Tanh(eta_0_rev) / Cos(beta_rev))
  var lambda = lambda_rad * (180/PI)
  
  return {
          'Latitude' : phi,
          'Longitude' : lambda,
  }
  }
//Example use converting from NAD83(HARN) / New Mexico Central to WGS84
convertTransverseMercator(easting, northing, 6378137, 298.257222101, 31, -106.25, 0.99992833, 500000, 0)

However, this formula is not 100% perfect. When working with the New Mexico Central data, I found that the converted coordinates became slightly less accurate the farther they were from the latitude of origin and the central meridian. I used a set of test data with Feature Geometry in State Plane and fields with accurate LatLon coordinates that I ran my arcade code on to see the difference in coordinate calculations. Since the points in my data were very far north of the projection's latitude of origin, the differences in Longitude were insignificant, while the Latitude calculations were off by ~0.000124 degrees.

I ran a simple linear regression between each feature's new calculated latitude and its distance from the more accurate conversion's latitude, which returned a nearly exact linear relationship as the latitude moved away from the latitude of origin.

Lat_diff = 0.0000303 * Latitude - 0.00094

R^2 = 0.9999988042

I added this code in at line 90 to correct  my phi values and returned new_phi as my latitude.

var new_phi = ((0.0000303*phi) - 0.00094) + phi

 My new latitude coordinates were then only off by ~0.000001 degrees from the previous converted coordinates. I am not sure how spatially sound this data manipulation method is for all projections, but it worked to correct the calculations on my data. I would suggest anyone attempting to use the larger code block functions with their own projection to have test data with LatLon coordinates already calculated to identify any linear patterns of difference and then modify their code accordingly.