Arcade Script For Converting From State Plane to Decimal Degrees

[RPGIS]

Hi,

I found this code that I copied from copilot and modified for my specific needs and works somewhat but needs some major corrections.

// Constants for NAD83 Georgia West (FIPS 1002)
var A = 6378137.0; // Semi-major axis for NAD83
var F = 1 / 298.257222101; // Flattening for NAD83
var E2 = (2 * F) - (F * F); // Eccentricity squared
var centralMeridian = -84.16666666666667 * (PI / 180); // Central meridian for Georgia West (degrees to radians)
var latitudeOfOrigin = 30.0 // Latitude of origin (degrees to radians)
var scaleFactor = 0.9999; // Scale factor
var falseEasting = 2296583.333333 * 0.3048006096012192; // False easting in meters (1 US foot = 0.3048006096012192 meters)
var falseNorthing = 0.0; // False northing in meters

// Function to convert State Plane (NAD83) to Latitude and Longitude
function statePlaneToLatLong(x, y) {
    x *= 0.3048006096012192; // Convert x from US Feet to meters
    y *= 0.3048006096012192; // Convert y from US Feet to meters
    x -= falseEasting;
    y -= falseNorthing;

    var m = y / scaleFactor;
    var mu = m / (A * (1 - (E2 / 4) - (3 * E2 * E2 / 64) - (5 * E2 * E2 * E2 / 256)));

    var e1 = (1 - sqrt(1 - E2)) / (1 + sqrt(1 - E2));
    
    var phi1 = mu + (3 * e1 / 2 - 27 * pow(e1, 3) / 32) * sin(2 * mu)
        + (21 * pow(e1, 2) / 16 - 55 * pow(e1, 4) / 32) * sin(4 * mu)
        + (151 * pow(e1, 3) / 96) * sin(6 * mu)
        + (1097 * pow(e1, 4) / 512) * sin(8 * mu);

    var sinPhi1 = sin(phi1);
    var cosPhi1 = cos(phi1);
    var t1 = tan(phi1) * tan(phi1);
    var c1 = E2 * cosPhi1 * cosPhi1 / (1 - E2);
    var r1 = A * (1 - E2) / pow(1 - E2 * sinPhi1 * sinPhi1, 1.5);
    var n1 = A / sqrt(1 - E2 * sinPhi1 * sinPhi1);
    var d = x / (n1 * scaleFactor);

    // Calculate latitude and longitude in radians
    var latitude = phi1 - (n1 * tan(phi1) / r1)
        * (d * d / 2 - (5 + 3 * t1 + 10 * c1 - 4 * c1 * c1 - 9 * E2) * pow(d, 4) / 24
        + (61 + 90 * t1 + 298 * c1 + 45 * t1 * t1 - 252 * E2 - 3 * c1 * c1) * pow(d, 6) / 720);

    var longitude = centralMeridian + (d - (1 + 2 * t1 + c1) * pow(d, 3) / 6
        + (5 - 2 * c1 + 28 * t1 - 3 * c1 * c1 + 8 * E2 + 24 * t1 * t1) * pow(d, 5) / 120) / cosPhi1;

    // Convert latitude and longitude to degrees
    var latitudeDeg = latitude * (180 / PI)*10;
    var longitudeDeg = longitude * (180 / PI);

    // Format results as degrees North and degrees West
    return {
        latitude: round( abs(latitudeDeg) , 6),
        longitude: round( abs(longitudeDeg) , 6)
    };
}

// Example usage with provided coordinates
var x = Geometry($feature).x; // X coordinate in US feet
var y = Geometry($feature).y; // Y coordinate in US feet

var result = statePlaneToLatLong(x, y);
console(result); // Output: {latitude: "xy.xyzxyz° N", longitude: "xy.xyzxyz° W"}
return result.latitude

[LindsayRaabe_FPCWA]

Can I suggest adding the purpose of the code to the post heading to make it more searchable? I believe it's a Lat/Long conversion code, which I know a lot of people would be looking for. Also maybe add some tags to improve search results. 

[RPGIS]

Thanks @LindsayRaabe_FPCWA for the suggestion. I went ahead and made the changes, and I am hoping that this can be used to create a function for attribute rules. I think there might be something simpler than this, but I will need to explore further.

[DavidColey]

Hello @RPGIS  - would you be able to share how you used copilot to derive your expression?  I ask because when I try to parametrize for Florida State Plane West HARN, the only differences are the centralMeridian, latitutudeOfOrigin and the scale factor.  

var centralMeridian = -82.0 * (PI / 180); // Central meridian for NAD83 HARN Florida State Plane West (degrees to radians)
var latitudeOfOrigin = 24.33333333333333 // Latitude of origin (degrees to radians)
var scaleFactor = 0.9999411764705882; // Scale factor

Other than the scale factor for HARN carrying more precision than the NAD 83 above I cannot seem to find where HARN would cause my conversions to have a 0.1 degree shift north in latitude and 0.05 shift east in longitude.

For example, when I plug in the expressions for my respective lat and lon attributes:

Pro: 27.078110470, -82.43048022

Arcade: 27.505790, -82.383999

Its interesting too because the console returns different values from the actual calculation:

{"latitude":27.851679,"longitude":-82.409632}

Any insights you can provide are appreciated, (I just don't have the trig to keep up) -

Thanks,

David

Link to a similar post with an Arcade expression that works for State Plane North / South coordinate systems (Lambert Conformal Conic)

https://community.esri.com/t5/arcgis-pro-ideas/arcade-projectas-geometry-function/idi-p/1171382/page...

@TomNeer 

[RPGIS]

I noticed that there some issues with the console values as well that I am trying to breakdown the math but this code was copied, like you said, using copilot but it is incorrect.

When I ran some tests it worked in some instances and not in others. I also noticed that sometimes it would be off by 4 degrees in some spots but I am trying to figure out where that issue is in the code.

[MDB_GIS]

There are a few issues that I've noticed. I actually just created basically the exact same version of this script using another geonet users post as a starting point, only their data was in Colorado state plane west. But I've managed to get it to work for Georgia. Here are the issues that I've noticed:

 1. Your latitude of origin is unused. It is needed to calculate the correct meridian arc.

//unused
latitudeOfOrigin = 30.0

2. Unless I'm missing something, you are multiplying your latitude by 10 for no reason which would just make your latitude 10 times too large in every scenario. 

var latitudeDeg = latitude * (180 / PI)*10;

3. You convert x,y to meters first, and then subtract the false easting, which was already converted to meters. It's a bit cleaner and easier to keep them in ft-us and then convert to meters at the end like so:

var ft_to_m = 0.3048006096012192;                 // US survey foot -> metervar E_m = (E - falseE_ft) * ft_to_m;
var N_m = (N - falseN_ft) * ft_to_m;

Here is my version which works for Georgia state plane west. Let me know if you have any questions. 

// 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

[RPGIS]

Thank you very much @MDB_GIS,

I will be sure to test this out and see how this works to align with our CAD coordinate system since that has been somewhat of a boondoggle. 

[MDB_GIS]

Great! Let me know how it goes!