Math for converting coordinates

[DuncanNisbett]

I know that I can use the geometryService to convert between 2 different coordinate systems, but I'd like to reduce the number of calls made to the server to update coordinates. I'm trying to convert from wkid 32149 (NAD_1983_StatePlane_Washington_South_FIPS_4602) to Latitude/Longitude 4326. Does anyone know where I can find the math so I can create my own function to convert the coordinates? I've tried my Google Kung Fu, but come up dry.

[JeffPace]

I know that I can use the geometryService to convert between 2 different coordinate systems, but I'd like to reduce the number of calls made to the server to update coordinates. I'm trying to convert from wkid 32149 (NAD_1983_StatePlane_Washington_South_FIPS_4602) to Latitude/Longitude 4326. Does anyone know where I can find the math so I can create my own function to convert the coordinates? I've tried my Google Kung Fu, but come up dry.

is washington south transverse mercator or lambert?

/* 
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */

dojo.provide("org.mymanatee.common.latlong");

dojo.require("dojo.string");

org.mymanatee.common.latlong.convertSPtoLL_LLC = function (uX,uY) {
var a = 20925604.48;     //major radius of ellipsoid, map units (NAD 83)
var ec = 0.08181905782;    //eccentricity of ellipsoid (NAD 83)
var angRad = 0.01745329252;   //number of radians in a degree
var pi4 = Math.PI / 4;   //Pi / 4
var p0 = 24.333333 * angRad;    //latitude of origin
var p1 = 24.333333 * angRad;    //latitude of first standard parallel
var p2 = 31.5 * angRad;   //latitude of second standard parallel
var m0 = -82.000000 * angRad;   //central meridian
var x0 = 656166.666667;   //False easting of central meridian, map units

// Calculate the coordinate system constants.
with (Math) {
var m1 = cos(p1) / sqrt(1 - (pow(ec,2)) * pow(sin(p1),2));
var m2 = cos(p2) / sqrt(1 - (pow(ec,2)) * pow(sin(p2),2));
var t0 = tan(pi4 - (p0 / 2));
var t1 = tan(pi4 - (p1 / 2));
var t2 = tan(pi4 - (p2 / 2));
t0 = t0 / pow(((1 - (ec * (sin(p0)))) / (1 + (ec * (sin(p0))))),ec/2);
t1 = t1 / pow(((1 - (ec * (sin(p1)))) / (1 + (ec * (sin(p1))))),ec/2);
t2 = t2 / pow(((1 - (ec * (sin(p2)))) / (1 + (ec * (sin(p2))))),ec/2);
var n = log(m1 / m2) / log(t1 / t2);
var f = m1 / (n * pow(t1,n));
var rho0 = a * f * pow(t0,n);

// Convert the coordinate to Latitude/Longitude.

// Calculate the Longitude.
var uX = uX - x0;
var pi2 = pi4 * 2;

var rho = sqrt(pow(uX,2) + pow((rho0 - uY),2));
var theta = atan(uX / (rho0 - uY));
var txy = pow((rho / (a * f)),(1 / n));
var lon = (theta / n) + m0;
var uX = uX + x0;

// Estimate the Latitude
var lat0 = pi2 - (2 * atan(txy));

// Substitute the estimate into the iterative calculation that
// converges on the correct Latitude value.
var part1 = (1 - (ec * sin(lat0))) / (1 + (ec * sin(lat0)));
var lat1 = pi2 - (2 * atan(txy * pow(part1,(ec/2))));

while ((abs(lat1 - lat0)) > 0.000000002) {
  lat0 = lat1;
  part1 = (1 - (ec * sin(lat0))) / (1 + (ec * sin(lat0)));
  lat1 = pi2 - (2 * atan(txy * pow(part1,(ec/2))));
  }

// Convert from radians to degrees.
var Lat = lat1 / angRad;
var Lon = lon / angRad;

//alert(Lat);
//alert(Lon);

//Round the latitude and longitude
//lat = (CLng(lat * 100000)) / 100000;
//lon = (CLng(lon * 100000)) / 100000;

Lat = Lat.toPrecision(7);
Lon = Lon.toPrecision(8);

return [Lat,Lon];

}
};
//End of Code for converting map units to decimal degrees.


//COnversion to LL from State Plane for Florida West (Transverse mercator, not Lamberic
org.mymanatee.common.latlong.convertSPtoLL_TM = function (uX,uY) {

    //code completely borrowed from Dr Bill Hazelton's HP-33S calculator
    // http://homepage.mac.com/nwjh/HP-33S/ 10-Feb-10
 var    metconv = 1200/3937; //exact survey feet to meters
 var    uXm = uX*metconv;
 var    uYm = uY*metconv;

  var   a = 6378137.000000000000000000; // semimajor axis from ESRI
 var    b =  6356752.314140356100000000; // semiminor axis from ESRI
 var    e2 = ((a*a)-(b*b))/(a*a);//ellipsoid
 var    p0 = 24.333333;    //latitude of origin
 var    m0 = -82.000000;   //central meridian in degrees
 var    x0 = 656166.666667;  //False easting of central meridian, map units
 var    y0 = 0; //False Northing
 var    k0 = 0.999941;    //central scale factor k0
 var    x0m = x0*metconv;
 var    y0m = y0*metconv;

 var    p0rad = p0*Math.PI/180;
 var    m0rad = m0*Math.PI/180;
 var    e = Math.sqrt(e2);


 var    b = a * Math.sqrt(1-e2);
 var    n = (a - b)/(a + b);
 var    G = a*(1-n)*(1-(n*n))*(1+2.25*(n*n)+225/64*(n*n*n*n))*Math.PI/180;
 var    A0 = 1-e2/4-3*(e2*e2)/64-5*(e2*e2*e2)/256;
 var    A2 = 3/8*(e2+(e2*e2)/4+15*(e2*e2*e2)/128);
var     A4 = 15/256*((e2*e2)+0.75*(e2*e2*e2));
var     A6 = 35*(e2*e2*e2)/3072;
var     lat0 = p0rad;
 var    m0atlat0 = a*(A0*lat0-A2*Math.sin(2*lat0)+A4*Math.sin(4*lat0)-A6*Math.sin(6*lat0))
 var    NminusNo = uYm - y0m;

var     m = m0atlat0+NminusNo/k0;
var     s = (m/G)*Math.PI/180;
 var    fpl1 = (1.5*n-27*(n*n*n)/32)*Math.sin(2*s);
  var   fpl2 = (21*(n*n)/16-55*(n*n*n*n)/32)*Math.sin(4*s);
 var    fpl3 = 151*(n*n*n)/96*Math.sin(6*s);
 var    fpl4 = 1097*(n*n*n*n)/512*Math.sin(8*s);
  var   footptlat = (s+fpl1+fpl2+fpl3+fpl4);
  var   rho = a*(1-e2)/Math.pow((1-e2*Math.sin(footptlat)*Math.sin(footptlat)),1.5);
 var    nu = a/(Math.sqrt(1-e2*Math.sin(footptlat)*Math.sin(footptlat)));
 var    psi = nu/rho;
  var   t = Math.tan(footptlat);
 var    E = uXm-x0m;
 var    x = E/k0/nu;

 var    lat1 = x*E/2;
 var    lat2 = (x*x*x)*E/24*(-4*psi*psi+9*psi*(1-t*t)+12*t*t);
  var   lat3 = (x*x*x*x*x)*E/720*(8*(psi*psi*psi*psi)*(11-24*(t*t))-12*(psi*psi*psi)*(21-71*(t*t))+15*(psi*psi)*(15-98*(t*t)+15*(t*t*t*t))+180*psi*(5*(t*t)-3*(t*t*t*t))+360*(t*t*t*t));
 var    lat4 = (x*x*x*x*x*x*x)*E/40320*(1385+3633*(t*t)+4095*(t*t*t*t)+1575*(t*t*t*t*t*t));
 var    lat = footptlat-(t/k0/rho)*(lat1-lat2+lat3-lat4);
 var    latdd = lat*180/Math.PI;

 var    w1 = x/Math.cos(footptlat);
 var    w2 = Math.pow(-1*x,3)/6*(psi+2*t*t)/Math.cos(footptlat);
 var    w3 = Math.pow(x,5)/120/Math.cos(footptlat)*(-4*(psi*psi*psi)*(1-6*(t*t))+(psi*psi)*(9-68*(t*t))+72*(psi)*(t*t)+24*(t*t*t*t));
 var    w4 = Math.pow(-1*x,7)/5040/Math.cos(footptlat)*(61+662*(t*t)+1320*(t*t*t*t)+720*(t*t*t*t*t*t));
 var    w = w1+w2+w3+w4;
 var    lon = (-1*m0)*Math.PI/180-w;
 var    londd = (-1*lon)*180/Math.PI;


//    alert(londd);

    return [latdd,londd];
    };

org.mymanatee.common.latlong.convertLLtoSP_TM = function (lat,lon) {
   var metconv = 1200/3937; //exact survey feet to meters

   var a = 6378137.000000000000000000; // semimajor axis from ESRI
   var b =  6356752.314140356100000000; // semiminor axis from ESRI
   var e2 = ((a*a)-(b*b))/(a*a);//ellipsoid
   var p0 = 24.333333;    //latitude of origin
   var m0 = -82.000000;   //central meridian in degrees
   var x0 = 656166.666667;  //False easting of central meridian, map units
   var y0 = 0; //False Northing
   var k0 = 0.999941;    //central scale factor k0
   var e = Math.sqrt(e2);
    
   var lat0rad_lat0=p0*Math.PI/180;
   var latrad_lat=lat*Math.PI/180;
   var rho_lat=a*(1-e2)/Math.pow((1-e2*Math.sin(latrad_lat)*Math.sin(latrad_lat)),1.5);
 //   rho_lat0=a*(1-e2)/Math.pow((1-e2*Math.sin(lat0rad_lat0)*Math.sin(lat0rad_lat0)),1.5); not used
 var   nu_lat=a/(Math.sqrt(1-e2*Math.sin(latrad_lat)*Math.sin(latrad_lat)));
 //   nu_lat0=a/(Math.sqrt(1-e2*Math.sin(lat0rad_lat0)*Math.sin(lat0rad_lat0))); not used
  var  omega_lat=((lon)*Math.PI/180)-((m0)*Math.PI/180); //provide negative lon so no negative multiplier
 //   omega_lat0=""; not used
 var   psi_lat=nu_lat/rho_lat;
 //   psi_lat0=""; not used
  var  t_lat=Math.tan(latrad_lat);
  //  t_lat0=""; not used
  var  a0_lat=1-e2/4-3*e2*e2/64-5*e2*e2*e2/256;
  var  a2_lat=3/8*(e2+e2*e2/4+15*e2*e2*e2/128);
  var  a4_lat=15/256*(e2*e2+0.75*e2*e2*e2);
  var  a6_lat=35*e2*e2*e2/3072;
  var  mlat_lat=a*(a0_lat*latrad_lat-a2_lat*Math.sin(2*latrad_lat)+a4_lat*Math.sin(4*latrad_lat)-a6_lat*Math.sin(6*latrad_lat));
  var  mlat_lat0=a*(a0_lat*lat0rad_lat0-a2_lat*Math.sin(2*lat0rad_lat0)+a4_lat*Math.sin(4*lat0rad_lat0)-a6_lat*Math.sin(6*lat0rad_lat0));
 var   e1_lat=nu_lat*omega_lat*Math.cos(latrad_lat);
 var   e2_lat=nu_lat*omega_lat*omega_lat*omega_lat/6*Math.cos(latrad_lat)*Math.cos(latrad_lat)*Math.cos(latrad_lat)*(psi_lat-t_lat*t_lat);
 var   e3_lat=nu_lat*Math.pow(omega_lat,5)/120*Math.pow(Math.cos(latrad_lat),5)*(4*psi_lat*psi_lat*psi_lat*(1-6*t_lat*t_lat)+psi_lat*psi_lat*(1+8*t_lat*t_lat)-psi_lat*2*t_lat*t_lat+Math.pow(t_lat,4));
  var  e4_lat=nu_lat*Math.pow(omega_lat,7)/5040*Math.pow(Math.cos(latrad_lat),7)*(61-479*t_lat*t_lat+179*Math.pow(t_lat,4)-Math.pow(t_lat,6));
 var   e_lat=k0*(e1_lat+e2_lat+e3_lat+e4_lat);
 var   n1_lat=(omega_lat*omega_lat)/2*Math.cos(latrad_lat);
 var   n1_lat0=n1_lat*nu_lat*Math.sin(latrad_lat);
 var   n2_lat=(Math.pow(omega_lat,4)/24)*(Math.pow(Math.cos(latrad_lat),3))*(4*psi_lat*psi_lat+psi_lat-t_lat*t_lat);
 var   n2_lat0=n2_lat*nu_lat*Math.sin(latrad_lat);
 var   n3_lat=(Math.pow(omega_lat,6)/720)*(Math.pow(Math.cos(latrad_lat),5))*(8*Math.pow(psi_lat,4)*(11-24*t_lat*t_lat)-28*Math.pow(psi_lat,3)*(1-6*t_lat*t_lat)+psi_lat*psi_lat*(1-32*t_lat*t_lat)-2*psi_lat*t_lat*t_lat+Math.pow(t_lat,4));
  var  n3_lat0=n3_lat*nu_lat*Math.sin(latrad_lat);
 var   n4_lat=(Math.pow(omega_lat,8)/40320)*(Math.pow(Math.cos(latrad_lat),7))*(1385-3111*t_lat*t_lat+543*Math.pow(t_lat,4)-Math.pow(t_lat,6));
 var   n4_lat0=n4_lat*nu_lat*Math.sin(latrad_lat);
 var   n_lat=k0*(mlat_lat+nu_lat*Math.sin(latrad_lat)*(n1_lat+n2_lat+n3_lat+n4_lat));
 var   n_lat0=k0*mlat_lat0;
 var   no=n_lat-n_lat0;
  var  x=e_lat+x0*metconv;
 var   y=(n_lat-n_lat0)+y0*metconv;
//    alert(n_lat);
//    alert("x="+x+"y="+y);
  //  alert("a0="+a0_lat+"a2="+a2_lat+"a4="+a4_lat+"a6="+a6_lat);
 var   xft=x/metconv;
 var   yft=y/metconv;
 //   alert("x="+xft+"y="+yft);
  return [xft,yft];
    }

[JeffPace]

also

    org.mymanatee.common.latlong.convertWebMercAuxtoLL = function (uX,uY) {
  //code completely borrowed from http://www.frosties.com/index.php?option=com_mojo&Itemid=45&p=7
   // lat and long backwards on link for some reason
   //alert("IN CONVERSION");
  var x = uX;
   // alert("IN CONVERSION2");
  var  y = uY;
   
 var num3 = x / 6378137.0;
 var num4 = num3 * 57.295779513082323;
 //alert("IN CONVERSION3");
 var num5 = Math.floor(((num4 + 180.0) / 360.0));
 var num6 = num4 - (num5 * 360.0);
 //alert("IN CONVERSION4");
var  num7 = 1.5707963267948966 - (2.0 * Math.atan(Math.exp((-1.0 * y) / 6378137.0)));
var lon  = num6;
var lat = num7*57.295779513082323;



   // alert([lat,lon]);


    return [lat,lon];
    };
    
        org.mymanatee.common.latlong.convertLLtoWebMercAux = function (lat,lon) {
  //code completely borrowed from Dr Bill Hazelton's HP-33S calculator
    // http://homepage.mac.com/nwjh/HP-33S/ 10-Feb-10

    //code completely borrowed from Dr Bill Hazelton's HP-33S calculator
    // http://homepage.mac.com/nwjh/HP-33S/ 10-Feb-10
var num = lon * 0.017453292519943295;
var x = 6378137.0 * num;
var a = lat * 0.017453292519943295;

var mercatorX = x;
var mercatorY = 3189068.5*Math.log((1.0 + Math.sin(a))/(1.0 - Math.sin(a)));

//    alert(londd);

    return [mercatorX,mercatorY];
    };

[DuncanNisbett]

is washington south transverse mercator or lambert?

It's Lambert Conformal Conic

[MelitaKennedy]

One of the original sources is John Snyder's Map Projections: A Working Manual. It also contains the state plane parameter values (NAD83 changes from NAD27 are listed in an appendix).

If you use Jeff Pace's code as a template, please note that the "map units" of the first section are in US survey feet, not meters.

Melita

[JeffPace]

that is correct

[DuncanNisbett]

One of the original sources is John Snyder's Map Projections: A Working Manual. It also contains the state plane parameter values (NAD83 changes from NAD27 are listed in an appendix).

If you use Jeff Pace's code as a template, please note that the "map units" of the first section are in US survey feet, not meters.

Melita

I admit to feeling a bit overwhelmed, this doesn't appear to be as straightforward a process as I had hoped. The PDF has a lot of information and am unable to process until I can really sit down and read it.

I took a look at the 5 functions that Jeff Pace included and no single function performs an accurate or even close conversion. Do I need to use a combination of these functions? If so, what combination should I be using?

[DavidKitten]

I stumbled upon this old post while working on a similar project and was able to get Jeff's code to work for the first function(convertSPtoLL_LLC) by adding in the false northing variable and Y coordinate correction.

New lines are in red.

1. org.mymanatee.common.latlong.convertSPtoLL_LLC = function (uX,uY) {  
2. var a = 20925604.48;     //major radius of ellipsoid, map units (NAD 83)  
3. var ec = 0.08181905782;    //eccentricity of ellipsoid (NAD 83)  
4. var angRad = 0.01745329252;   //number of radians in a degree  
5. var pi4 = Math.PI / 4;   //Pi / 4  
6. var p0 = 24.333333 * angRad;    //latitude of origin  
7. var p1 = 24.333333 * angRad;    //latitude of first standard parallel  
8. var p2 = 31.5 * angRad;   //latitude of second standard parallel  
9. var m0 = -82.000000 * angRad;   //central meridian  
10. var x0 = 656166.666667;   //False easting of central meridian, map units  
11. var y0 = value;       //False northing of central meridian, map units  
12. // Calculate the coordinate system constants.  
13. with (Math) {  
14. var m1 = cos(p1) / sqrt(1 - (pow(ec,2)) * pow(sin(p1),2));  
15. var m2 = cos(p2) / sqrt(1 - (pow(ec,2)) * pow(sin(p2),2));  
16. var t0 = tan(pi4 - (p0 / 2));  
17. var t1 = tan(pi4 - (p1 / 2));  
18. var t2 = tan(pi4 - (p2 / 2));  
19. t0 = t0 / pow(((1 - (ec * (sin(p0)))) / (1 + (ec * (sin(p0))))),ec/2);  
20. t1 = t1 / pow(((1 - (ec * (sin(p1)))) / (1 + (ec * (sin(p1))))),ec/2);  
21. t2 = t2 / pow(((1 - (ec * (sin(p2)))) / (1 + (ec * (sin(p2))))),ec/2);  
22. var n = log(m1 / m2) / log(t1 / t2);  
23. var f = m1 / (n * pow(t1,n));  
24. var rho0 = a * f * pow(t0,n);  
25.   
26. // Convert the coordinate to Latitude/Longitude.  
27.   
28. // Calculate the Longitude.  
29. var uX = uX - x0; 
30. var uY = uY - y0; 
31. var pi2 = pi4 * 2;  
32.   
33. var rho = sqrt(pow(uX,2) + pow((rho0 - uY),2));  
34. var theta = atan(uX / (rho0 - uY));  
35. var txy = pow((rho / (a * f)),(1 / n));  
36. var lon = (theta / n) + m0;  
37. var uX = uX + x0;  
38.   
39. // Estimate the Latitude  
40. var lat0 = pi2 - (2 * atan(txy));  
41.   
42. // Substitute the estimate into the iterative calculation that  
43. // converges on the correct Latitude value.  
44. var part1 = (1 - (ec * sin(lat0))) / (1 + (ec * sin(lat0)));  
45. var lat1 = pi2 - (2 * atan(txy * pow(part1,(ec/2))));  
46.   
47. while ((abs(lat1 - lat0)) > 0.000000002) {  
48.   lat0 = lat1;  
49.   part1 = (1 - (ec * sin(lat0))) / (1 + (ec * sin(lat0)));  
50.   lat1 = pi2 - (2 * atan(txy * pow(part1,(ec/2))));  
51.   }  
52.   
53. // Convert from radians to degrees.  
54. var Lat = lat1 / angRad;  
55. var Lon = lon / angRad;  
56.   
57. //alert(Lat);  
58. //alert(Lon);  
59.   
60. //Round the latitude and longitude  
61. //lat = (CLng(lat * 100000)) / 100000;  
62. //lon = (CLng(lon * 100000)) / 100000;  
63.   
64. Lat = Lat.toPrecision(7);  
65. Lon = Lon.toPrecision(8);  
66.   
67. return [Lat,Lon];  
68.   
69. }  
70. };  
71. //End of Code for converting map units to decimal degrees.