Automatically Calculate Geometry of Points in Web Mercator

[BenWagner3]

I'm working in an Enterprise Geodatabase (Enterprise Version 11.3), and I have a point feature class that is in Iowa State Plane North. I currently have attribute rules that run on Insert and Update to calculate the Northing and Easting of the points. I have been asked to also include attribute fields called Latitude and Longitude and have those fields automatically calculate the Y & X values in Web Mercator. I don't believe that Arcade in an attribute rule can handle that process. I have tried creating a scheduled process to calculate the geometry, but so far my Python scripts have been generating values that aren't remotely close to being correct.

I would prefer not to have to manually run the Calculate Geometry tool, like I have the past couple of days.

Does anyone have pointers on how to write a Python script that would calculate the X & Y in Web Mercator?

Thanks!

[CameronRex1]

@BenWagner3 You can try the following (I am not entirely sure what Iowa State Plane North you are using, so I have assumed WKID: 3417). I have tested something similar in another state plane and was able to get correct results. (If you are using a different WKID, let me know and I can change it up.)

// === Projection parameters for NAD83 StatePlane Iowa North FIPS 1401 (US Feet) ===
var lat0 = 41.5;
var lon0 = -93.5;
var phi1 = 43.2666666666667;
var phi2 = 42.0666666666667;
var E0 = 4921250.0; // False easting (ft)
var N0 = 3280833.3333; // False northing (ft)

// GRS80 ellipsoid
var a = 6378137.0;
var f = 1 / 298.257222101;
var e = sqrt(2 * f - f * f);

// === Helper functions ===
function rad(d) { return d * PI / 180; }
function deg(r) { return r * 180 / PI; }
function m(phi) {
return cos(phi) / sqrt(1 - pow(e, 2) * pow(sin(phi), 2));
}
function t(phi) {
return tan(PI / 4 - phi / 2) / pow((1 - e * sin(phi)) / (1 + e * sin(phi)), e / 2);
}
function inverseLCC(x, y) {
var phi1_rad = rad(phi1);
var phi2_rad = rad(phi2);
var lat0_rad = rad(lat0);
var lon0_rad = rad(lon0);

var m1 = m(phi1_rad);
var m2 = m(phi2_rad);
var t0 = t(lat0_rad);
var t1 = t(phi1_rad);
var t2 = t(phi2_rad);

var n = (log(m1) - log(m2)) / (log(t1) - log(t2));
var F = m1 / (n * pow(t1, n));
var r0 = a * F * pow(t0, n);

var r = sqrt(pow(x - E0, 2) + pow(r0 - (y - N0), 2));
var t_p = pow(r / (a * F), 1 / n);

var chi = PI / 2 - 2 * atan(t_p);
var lat_rad = chi + (e*e/2 + 5*pow(e,4)/24 + pow(e,6)/12) * sin(2*chi)
+ (7*pow(e,4)/48 + 29*pow(e,6)/240) * sin(4*chi)
+ (7*pow(e,6)/120) * sin(6*chi);
var lat = deg(lat_rad);

var theta = atan((x - E0) / (r0 - (y - N0)));
var lon = lon0 + deg(theta / n);

return [lat, lon];
}

// === Get point coordinates from the feature ===
var x = $feature.x;
var y = $feature.y;

// === Convert to Lat/Lon
var latlon = inverseLCC(x, y);

// === Format output
"Lat: " + Text(latlon[0], "#.######") + ", Lon: " + Text(latlon[1], "#.######")

[MOBrien2]

Is there any chance you could do this with WKID 2260 (New York East US Feet)? It uses Transverse Mercator projection instead of Lambert Conic Conformal (2SP). I can copy all the details down from the projected coordinate system definition, but Transverse Mercator uses scale factor instead of standard parallels which means the phi variables wouldn't work the same. I wish Esri would come up with a set of scripts for the different projection methods to do this.

[CameronRex1]

@MOBrien2 Sorry I missed your post. I think the following code would work, but I have not fully tested.

// === Arcade Expression ===
// === NAD83 / New York East (ftUS) — EPSG:2260 ===
// Inverse Transverse Mercator to get Lat/Lon, then Web Mercator XY (optional).

// ----- Zone parameters (EPSG 2260) -----
var lat0_deg = 38.8333333333333;   // latitude of origin
var lon0_deg = -74.5;              // central meridian
var k0 = 0.9999;                   // scale factor at origin
var FE = 492125.0;                 // false easting (ftUS)
var FN = 0.0;                      // false northing (ftUS)

// ----- Ellipsoid (GRS80) -----
var a_m = 6378137.0;
var inv_f = 298.257222101;
var f = 1.0 / inv_f;
var e2 = 2*f - f*f;                // eccentricity squared
var e = Sqrt(e2);

// Use feet consistently (layer is ftUS)
var M_PER_FTUS = 0.3048006096012192;
var a_ft = a_m / M_PER_FTUS;

// Helpers
function rad(d){ return d * PI / 180.0; }
function deg(r){ return r * 180.0 / PI; }
function wrapLon(d){ return (d + 540) % 360 - 180; }

// --- Forward Web Mercator from lat/lon (deg) ---
function webMercatorXY(lat_deg, lon_deg){
  var maxLat = 85.0511287798066;
  var φ = rad(Max(-maxLat, Min(maxLat, lat_deg)));
  var λ = rad(lon_deg);
  var R = 6378137.0;
  return [ R * λ, R * Log(Tan(PI/4 + φ/2)) ];
}

// --- Inverse Transverse Mercator (Redfearn) ---
function inverseTM_2260(E, N){
  // Remove false origin
  var x = E - FE;
  var y = N - FN;

  // Footpoint latitude via meridional arc
  var A0 = 1 - e2/4 - 3*Pow(e2,2)/64 - 5*Pow(e2,3)/256;
  var A2 = 3/8 * (e2 + Pow(e2,2)/4 + 15*Pow(e2,3)/128);
  var A4 = 15/256 * (Pow(e2,2) + 3*Pow(e2,3)/4);
  var A6 = 35*Pow(e2,3)/3072;

  var M = y / k0;                 // meridional arc (ft)
  var mu = M / (a_ft * A0);

  // Series for footprint latitude φ1
  var e1 = (1 - Sqrt(1 - e2)) / (1 + Sqrt(1 - e2));
  var J1 = 3/2*e1 - 27/32*Pow(e1,3);
  var J2 = 21/16*Pow(e1,2) - 55/32*Pow(e1,4);
  var J3 = 151/96*Pow(e1,3);
  var J4 = 1097/512*Pow(e1,4);
  var phi1 = mu + J1*Sin(2*mu) + J2*Sin(4*mu) + J3*Sin(6*mu) + J4*Sin(8*mu);

  // Radii of curvature and helpers at φ1
  var sin1 = Sin(phi1);
  var cos1 = Cos(phi1);
  var tan1 = Tan(phi1);
  var e2p = e2 / (1 - e2);                    // second eccentricity squared
  var N1 = a_ft / Sqrt(1 - e2 * sin1*sin1);   // prime vertical radius (ft)
  var R1 = N1 * (1 - e2) / (1 - e2 * sin1*sin1); // meridional radius (ft)
  var T1 = tan1*tan1;
  var C1 = e2p * cos1*cos1;
  var D  = x / (N1 * k0);

  // Latitude
  var lat_rad = phi1
    - (N1 * tan1 / R1) * (
        (D*D)/2
        - Pow(D,4)/24 * (5 + 3*T1 + 10*C1 - 4*C1*C1 - 9*e2p)
        + Pow(D,6)/720 * (61 + 90*T1 + 298*C1 + 45*T1*T1 - 252*e2p - 3*C1*C1)
      );

  // Longitude
  var lon_rad = rad(lon0_deg)
    + ( D
        - Pow(D,3)/6 * (1 + 2*T1 + C1)
        + Pow(D,5)/120 * (5 + 28*T1 + 24*T1*T1 + 6*C1 + 8*C1*C1)
      ) / cos1;

  return [ deg(lat_rad), wrapLon(deg(lon_rad)) ];
}

// === main ===
var g = Geometry($feature);
if (IsEmpty(g)) return null;

var E = g.x;  // ftUS
var N = g.y;

var latlon = inverseTM_2260(E, N);
var lat = Round(latlon[0], 6);
var lon = Round(latlon[1], 6);

// For a popup:
"Lat: " + Text(lat, "#.######") + ", Lon: " + Text(lon, "#.######")

// For an Attribute Rule, return a dict instead:
// return { "Latitude": lat, "Longitude": lon };

[MOBrien2]

Thank you for responding to me with this, I'll try to test it later. Shortly after I made that reply I was eventually able to implement my own solution through ArcGIS Arcade, shown in this post: https://community.esri.com/t5/arcgis-pro-questions/is-there-a-way-to-automatically-update-latitude/m... 
which uses the following Arcade script:

//Transverse Mercator projected state plane coordinates converted to geodetic coordinates using these equations should be accurate to the milimeter or submilimeter. If modifying this script, test its accuracy before implementing.

// === Projection parameters for NAD83 State Plane New York East FIPS 3101 (US Survey Feet) (WKID 2260)===
//Change the variables in this section and others if using a different state plane or ellipsoid.

//US Survey foot compared to meters
var linear_unit = 0.3048006096012192; 
//Latitude of origin (decimal degrees) (38°50'N)
//var latitude_of_origin = 38.8333333333333;
//Commented out as this script does not rely on the latitude of origin
//Longitude of origin and Central Meridian (decimal degrees) (74°30'W)
var central_meridian = -74.5;
//Scale factor at natural origin/central meridian
var scale_factor_origin = 0.999900;
//False Easting (meters) (492125.0 in US Survey Feet)
var false_easting = 150000.0;
//False Northing (meters)
var false_northing = 0.0; 
//Meridional distance from the equator to latitude of origin, multiplied by the central meridian scale factor (written as S subscript o in documentation) (meters)
var meridional_distance_origin = 4299571.6693;

//GRS80 ellipsoid
var semimajor_axis = 6378137.0; 
//in meters
var semiminor_axis = 6356752.314140356; 
//in meters

//Helper functions
function degToRad(deg){return (deg*PI/180)};
function radToDeg(rad){return (rad*180/PI)};
function USsurveyFtToMeters(ft){return ft*linear_unit};

//Convert projection information into radians
var cmRad = degToRad(central_meridian);

//== Ellipsoid definitions ==
//Flattening of the ellipsoid
var f = (semimajor_axis - semiminor_axis) / semimajor_axis;
//Inverse flattening, should be 298.25722210088 or close to it
//var invFlat = 1 / f; //Commented out as no other expressions use it in this script
//First eccentricity squared, should be 0.0066943800229034 or close to it
var ecc_sq = 2*f-Pow(f,2);
//Second eccentricity squared
var secondEcc_sq = ecc_sq/(1-ecc_sq);

//Equations and constants adapted from https://geodesy.noaa.gov/library/pdfs/NOAA_Manual_NOS_NGS_0005.pdf
//REMEMBER TO CHECK IF LONGITUDE/X-COORDINATE FORMULA AND CENTRAL MERIDIAN USES POSITIVE WEST OR POSITIVE EAST!
//https://geodesy.noaa.gov/library/pdfs/Transverse_Mercator_Projection_Tables_New_York.pdf
//Projection constants
var r_s = 6367449.14577;//Radius of the rectifying sphere (GRS 80, meters)
var V_0 = 0.005022893948;
var V_2 = 0.000029370625;
var V_4 = 0.000000235059;
var V_6 = 0.000000002181;

function reverseTM(X,Y){
//Equations taken from and variables based on those used in chapter 3.2, Transverse Mercator Mapping Equations, of NOAA Manual NOS NGS 5, State Plane Coordinate System of 1983, linked above
//Linear units in these equations are the units of false northing/easting and ellipsoid (meters)
//All angles are in radians
var X_m = USsurveyFtToMeters(X);
var Y_m = USsurveyFtToMeters(Y);
var rectifying_latitude = (Y_m - false_northing + meridional_distance_origin)/(scale_factor_origin*r_s);
var footprint_lat = rectifying_latitude + (sin(rectifying_latitude)*cos(rectifying_latitude)) * (V_0 + Pow(cos(rectifying_latitude),2) * (V_2 + Pow(cos(rectifying_latitude),2) * (V_4 + V_6 * Pow(cos(rectifying_latitude),2))));
var footprint_R = scale_factor_origin * semimajor_axis / Pow(1 - ecc_sq * Pow(sin(footprint_lat),2),0.5);
var Q = (X_m - false_easting) / footprint_R;
var t_f = tan(footprint_lat); //tangent of footprint latitude
var Eta2_f = secondEcc_sq*Pow(cos(footprint_lat),2); //lowercase Eta squared, using footprint latitude
var B_2 = -0.5*t_f*(1+Eta2_f);
var B_4 = (-1/12)*(5+(3*Pow(t_f,2))+(Eta2_f*(1-(9*Pow(t_f,2))))-(4*Pow(Eta2_f,2)));
var B_6 = (1/360)*(61+(90*Pow(t_f,2))+(45*Pow(t_f,4))+(Eta2_f*(46-(252*Pow(t_f,2))-(90*Pow(t_f,4)))));
//latitude in radians; positive is north
var latitude_rad = footprint_lat+B_2*Pow(Q,2)*(Pow(Q,2)*(B_6*Pow(Q,2)+B_4)+1);
var B_3 = (-1/6)*(1+(2*Pow(t_f,2))+Eta2_f);
var B_5 = (1/120)*(5+(28*Pow(t_f,2))+(24*Pow(t_f,4))+(Eta2_f*(8*Pow(t_f,2)+6)));
var B_7 = (-1/5040)*(61+(662*Pow(t_f,2))+(1320*Pow(t_f,4))+(720*Pow(t_f,6)));
var L = Q*(Pow(Q,2)*(Pow(Q,2)*(B_7*Pow(Q,2)+B_5)+B_3)+1);
//longitude in radians, positive is east; longitude_rad should use subtraction instead of addition if central_meridian is positive.
var longitude_rad = cmRad+(L/cos(footprint_lat));

var Lat = radToDeg(latitude_rad);
var Lon = radToDeg(longitude_rad);

return [Lat,Lon];
}

// === Get point coordinates from the feature ===
var x = Geometry($feature).x;
var y = Geometry($feature).y;

//Convert to Lat/Lon
var latlon = reverseTM(x,y);

//Output for attribute rules:
return{
        "result":{
             "attributes":{
                   "LAT_CALCULATED": latlon[0],
                   "LON_CALCULATED": latlon[1] //Use the field names for latitude and longitude in the double quotes
                   }
                }
}

[BenWagner3]

I'm just now getting back to this project. The script you did is not quite working. Punching the returned coords into Google Maps takes me out to a spot in Wisconsin. I'm working with a sample parking sign. Here is the correct coordinates in Iowa State Plane North and lat/long:

Here is the lat/long that is being returned by the script:

Here is the spatial reference of the feature class:

[CameronRex1 · Accepted solution]

Wow, totally don't remember this post. But it appears that I found a working solution. In my quick test with your exact point, it works. (Note: I just did an edit to pull from the Geometry, instead of a calculated value in the attribute table).

Try This:

// === Arcade Expression ===
// Source SR: NAD83 / Iowa North (ftUS) — WKID 3417
// LCC(2SP) params from EPSG: lat0=41.5°, lon0=-93.5°, sp1=43°16'00", sp2=42°04'00"
// False Easting/Northing in US survey feet.

var FE_ft = 4921250.0;        // false easting (ftUS)
var FN_ft = 3280833.3333;     // false northing (ftUS)
var lat0_deg = 41.5;          // latitude of origin
var lon0_deg = -93.5;         // central meridian
var sp1_deg  = 43.2666666667; // standard parallel 1 (43°16')
var sp2_deg  = 42.0666666667; // standard parallel 2 (42°04')

// GRS80 ellipsoid (StatePlane NAD83 uses GRS80)
var a_m = 6378137.0;                 // semi-major (meters)
var f   = 1.0 / 298.257222101;       // flattening
var e   = Sqrt(2*f - f*f);           // eccentricity

// Units: layer is in US survey feet; convert 'a' to feet so math is consistent with x/y
var M_PER_FTUS = 0.3048006096012192;
var a_ft = a_m / M_PER_FTUS;   // semi-major in ftUS

// --- helpers ---
function rad(d){ return d * PI / 180.0; }
function deg(r){ return r * 180.0 / PI; }

function m_phi(phi){
  return Cos(phi) / Sqrt(1 - (e*e) * Pow(Sin(phi),2));
}
function t_phi(phi){
  return Tan(PI/4 - phi/2) / Pow((1 - e*Sin(phi)) / (1 + e*Sin(phi)), e/2);
}

// Inverse Lambert Conformal Conic (2SP) -> lat/lon (degrees)
function inverseLCC_3417(x_ft, y_ft){
  var phi1 = rad(sp1_deg);
  var phi2 = rad(sp2_deg);
  var phi0 = rad(lat0_deg);
  var lam0 = rad(lon0_deg);

  var m1 = m_phi(phi1);
  var m2 = m_phi(phi2);
  var t0 = t_phi(phi0);
  var t1 = t_phi(phi1);
  var t2 = t_phi(phi2);

  var n = (Log(m1) - Log(m2)) / (Log(t1) - Log(t2));
  var F = m1 / (n * Pow(t1, n));
  var r0 = a_ft * F * Pow(t0, n);

  var dx = x_ft - FE_ft;
  var dy = r0 - (y_ft - FN_ft);
  var r  = Sqrt(dx*dx + dy*dy);

  var tp = Pow(r / (a_ft * F), 1.0/n);

  // Snyder series: chi -> phi
  var chi = PI/2 - 2 * Atan(tp);
  var lat_rad =
      chi
    + (e*e/2 + 5*Pow(e,4)/24 + Pow(e,6)/12) * Sin(2*chi)
    + (7*Pow(e,4)/48 + 29*Pow(e,6)/240) * Sin(4*chi)
    + (7*Pow(e,6)/120) * Sin(6*chi);

  var theta = Atan(dx / dy); // dy won't be 0 except on a special arc
  var lon_rad = lam0 + theta / n;

  return [ deg(lat_rad), deg(lon_rad) ];
}

// Forward Web Mercator (Pseudo-Mercator) from lat/lon (degrees)
function webMercatorXY(lat_deg, lon_deg){
  // clamp latitude to the projection’s practical limit
  var maxLat = 85.0511287798066;
  var φ = rad(Max(-maxLat, Min(maxLat, lat_deg)));
  var λ = rad(lon_deg);
  var R = 6378137.0; // sphere radius used by EPSG:3857

  var x = R * λ;
  var y = R * Log(Tan(PI/4 + φ/2));
  return [x, y];
}

// === main ===
var pt = Geometry($feature)
var x_ft = pt.X;  // Iowa North feet
var y_ft = pt.Y;

var latlon = inverseLCC_3417(x_ft, y_ft);
var lat = Round(latlon[0], 6);
var lon = Round(latlon[1], 6);

// === Format output
"Lat: " + Text(lat, "#.######") + ", Lon: " + Text(lon, "#.######")

[BenWagner3]

Yeah, I started on this project several months ago, but then it got shelved for more pressing items and I'm just now getting back to it. This new script did the job! Thanks!