Arcade projectAs Geometry Function

[JacobFougstrupNicolajsen]

Working primarily with data in UTM projection, I often need ways to convert UTM coordinates to WGS to create deep links to various online services. 

Having to create duplicate datasets, or running python scripts as regular intervals is tedious and time consuming.

And arcade equivalent of the python projectAs method would solve almost all these issues.

[tleeming]

This would be great for calculating meaningful area values in Arcade attribute rules.

[OlivierBaudry]

+1 for this idea!

Working with a team of field workers creating and updating data continiously make this function highly necessary in order to have correct DD coordinate field from any system.

I'm currently looking for a reliable way to extract DD geometry from Lambert layer (EPSG:3812).

Do you know more about update of this idea @CraigWilliams ?

[MDB_GIS]

@TomNeer Thank you for your excellent work. I attempted to use this for to reverse engineer a way to do this in Georgia. It took a bit of work, because apparently Colorado's state plane projection is Lambert. Georgia's is a Transverse_Mercator, however. So the math is completely different. Just in case anyone else finds themselves in the same predicament here is my modified version of Tom's code that works for Georgia State Plane West. If you use a state plane coordinate system with a projection of Transverse_Mercator, you just need to update the values  under "Projection Parameters" with the values from your coordinate systems information. You can find it on ESPG's website or by looking at the projection properties inside arcgis pro. Hope this helps!

// Inverse Transverse Mercator for EPSG:2240 (NAD83 / Georgia West (ftUS))
// Returns latitude, longitude in decimal degrees

// Input projected coordinates (assumes geometry is in EPSG:2240)
var E = Geometry($feature).x; // Easting (ftUS)
var N = Geometry($feature).y; // Northing (ftUS)

// Projection parameters from EPSG:2240
var a = 6378137.0;                                // semimajor axis (meters)
var inv_f = 298.257222101;                        // inverse flattening
var f = 1 / inv_f;
var e2 = 2 * f - f * f;                           // eccentricity squared

var k0 = 0.9999;                                  // scale factor
var lon0 = -84.1666666666667;                     // central meridian (deg)
var lat0 = 30.0;                                  // latitude of origin (deg)
var falseE_ft = 2296583.333;                      // false easting (ftUS)
var falseN_ft = 0.0;                              // false northing (ftUS)
var ft_to_m = 0.3048006096012192;                 // US survey foot -> meter

// helpers
function degToRad(d) { return d * PI / 180.0; }
function radToDeg(r) { return r * 180.0 / PI; }

// convert inputs into meters and remove false easting/northing
var E_m = (E - falseE_ft) * ft_to_m;
var N_m = (N - falseN_ft) * ft_to_m;

// --- compute meridian arc for latitude of origin (M0) ---
function meridianArc(phi) {
  // phi in radians
  var n2 = e2;
  var a0 = 1 - (n2/4) - (3*Pow(n2,2)/64) - (5*Pow(n2,3)/256);
  var a2 = (3.0/8.0) * (n2 + (Pow(n2,2)/4) + (15*Pow(n2,3)/128));
  var a4 = (15.0/256.0) * (Pow(n2,2) + (3*Pow(n2,3)/4));
  var a6 = (35.0/3072.0) * Pow(n2,3);
  return a * (a0 * phi - a2 * Sin(2*phi) + a4 * Sin(4*phi) - a6 * Sin(6*phi));
}

var lat0R = degToRad(lat0);
var M0 = meridianArc(lat0R);

// --- compute M (meridian arc to footpoint) from N_m ---
var M = M0 + (N_m / k0);

// --- Compute footpoint latitude (phi1) from M using series ---
var mu = M / (a * (1 - e2/4 - 3*Pow(e2,2)/64 - 5*Pow(e2,3)/256));

// e1 (called e' or e1 in many references)
var e1 = (1.0 - Sqrt(1.0 - e2)) / (1.0 + Sqrt(1.0 - e2));

// series coefficients for footpoint latitude
var J1 = (3*e1/2) - (27*Pow(e1,3)/32);
var J2 = (21*Pow(e1,2)/16) - (55*Pow(e1,4)/32);
var J3 = (151*Pow(e1,3)/96);
var J4 = (1097*Pow(e1,4)/512);

var phi1 = mu + J1 * Sin(2*mu) + J2 * Sin(4*mu) + J3 * Sin(6*mu) + J4 * Sin(8*mu);

// --- compute parameters at phi1 ---
var sin1 = Sin(phi1);
var cos1 = Cos(phi1);
var tan1 = Tan(phi1);

var N1 = a / Sqrt(1 - e2 * Pow(sin1,2));                   // radius of curvature in prime vertical
var R1 = a * (1 - e2) / Pow(1 - e2 * Pow(sin1,2), 1.5);    // meridional radius of curvature
var T1 = Pow(tan1,2);
var C1 = (e2 / (1 - e2)) * Pow(cos1,2);                    // e'² * cos²(phi1)
var D = E_m / (N1 * k0);

// --- Series expansion to get latitude (phi) ---
var term1 = (D*D / 2.0);
var term2 = (5 + 3*T1 + 10*C1 - 4*Pow(C1,2) - 9*(e2/(1-e2))) * Pow(D,4) / 24.0;
var term3 = (61 + 90*T1 + 298*C1 + 45*Pow(T1,2) - 252*(e2/(1-e2)) - 3*Pow(C1,2)) * Pow(D,6) / 720.0;

var phi = phi1 - (N1 * tan1 / R1) * (term1 - term2 + term3);

// --- Series expansion to get longitude ---
var lon_term1 = D;
var lon_term2 = (1 + 2*T1 + C1) * Pow(D,3) / 6.0;
var lon_term3 = (5 - 2*C1 + 28*T1 - 3*Pow(C1,2) + 8*(e2/(1-e2)) + 24*Pow(T1,2)) * Pow(D,5) / 120.0;

var lon = degToRad(lon0) + (lon_term1 - lon_term2 + lon_term3) / cos1;

// convert to decimal degrees
var lat_dd = radToDeg(phi);
var lon_dd = radToDeg(lon);

// change me
return lat_dd

[BlakeTerhune]

@MDB_GIS, I'm attempting to use your code with NAD83 HARN / Arizona Central ft (EPSG WKID 2868), but the lon, lat values are still slightly off (mostly the longitude). Here's the parameters I changed. Can you tell where I'm going wrong? AI was no help for me on this one.

// Projection parameters from EPSG:2868
var a = 6378137.0;                                // semimajor axis (meters)
var inv_f = 298.257222101;                        // inverse flattening
var f = 1 / inv_f;
var e2 = 2 * f - f * f;                           // eccentricity squared

var k0 = 0.9999;                                  // scale factor
var lon0 = -111.9166666666667;                    // central meridian (deg)
var lat0 = 31.0;                                  // latitude of origin (deg)
var falseE_ft = 700000;                           // false easting (Intl Ft)
var falseN_ft = 0.0;                              // false northing (Intl Ft)
var ft_to_m = 0.3048006096012192;                 // US survey foot -> meter

[TomNeerDDS]

@BlakeTerhune  - Reviewing EPSG 2868, it appears to use an International Foot. Try changing your ft_to_m variable to "0.3048". 

With that said, there is a HARN definition in here that I am not sure I nor @MDB_GIS may have dealt with. If the change to the foot doesn't work, I would start looking at how to implement the HARN definition.

[DavidColey]

Hmm, I can tell you guys that last week I modified the projection parameters for WKID 2882:

Projected Coordinate System NAD 1983 HARN StatePlane Florida West FIPS 0902 (US Feet)
Projection Transverse Mercator
WKID 2882
Authority EPSG
Linear Unit US Survey Feet 0.3048006096012192
False Easting 656166.6666666665
False Northing 0.0
Central Meridian -82.0
Scale Factor 0.9999411764705882
Latitude Of Origin 24.33333333333333

and tested the expression in an ArcGIS Pro 3.6.2 field calculation expressions for lat and then for lon.

I then ran the Calculate Geometry tool and the results matched.

Both my layer and the pro map coordinate system for this test are in WKID 2882. 

I plan on setting up the expressions for the web map edit form and testing against a web mercator basemap this week in both the webmap and in Field Maps.

[BlakeTerhune]

Interesting. I had only been testing in ArcGIS Online. In ArcGIS Pro, the calculation (with ft_to_m = 0.3048) is correct. Has anyone else tried this in ArcGIS Online? Here's a comparison of what I'm seeing between the two apps (at the same location):

ArcGIS Pro:
-111.8669267268508, 33.31052991092958

ArcGIS Online:
-154.5911295109528, 33.42959248309001

[DavidColey]

My guess is is that the webmaps' basemap must also be in the same WKID, but I don't know yet

[MDB_GIS]

@BlakeTerhune 

This is the ESPG page for 2868. It confirms that it is Transverse Mercator, so the math part should be correct. @DavidColey is correct that your system is based on international feet while mine is based on survey feet. It is a very small difference but it would definitely add up. Change var ft_to_m to .3048 and see if that corrects the issue. 

Regarding using it in AGOL, I haven't tried that yet. My entire enterprise is set up in portal. One thing you can do to quickly check if this should work, however, is set a quick label expression. Use Geometry($feature).x or Geometry($feature).y to have it return the xy geometry of your points. Your Easting should be somewhere around 700,000 feet. Northing will be in the millions. If both x/y are in the millions, positive or negative, then your data is using web Mercator and this expression wouldn't work. 

[DavidHorwood]

@BlakeTerhune, it appears @MDB_GIS  and @DavidColey are on the right track. AGOL is published in Web Mercator (3857 or 102100). If you take the Web Mercator coordinates of the the first geographic coordinates (ArcGIS Pro) and enter them as 2868 coordinates you will get the second set of coordinates (AGOL).

If you are just trying to get latitude and longitude from Web Mercator coordinates, it would be a different calculation (but should be easier, spherical rather than elliptical).

One thing of note: it appears 2868 implements a geographic transformation from NAD83 to WGS84.

TOWGS84[-0.991,1.9072,0.5129,-1.25033E-07,-4.6785E-08,-5.6529E-08,0]]

This appears to correspond to NAD_1983_HARN_To_WGS_1984_3 (1901) which is essentially NAD 1983 to ITRF 1996. It will result in a displacement of between one to 1.5 meters (3 to 5 feet) depending where you are in North America. This will complicate the calculation somewhat, since you will be returning WGS84 latitude and longitude which will need to be transformed to NAD83 latitude and longitude.

It would be really helpful if Esri implemented projectAs in Arcade. When implementing parcel fabrics for large jurisdictions across multiple zones (UTM or local projections), we regularly need to convert between systems on the fly, particularly for labelling and validation.