nitsrik

Convert Decimal Degrees to State Plane (Javascript)

Blog Post created by nitsrik on Jul 8, 2019

In case anybody needs it.

 

// Code to convert from decimal degrees to EPSG:2227 (feet)
// Code based on asp code from Gerry Daumiller, Montana State Library 8-23-02 email. ASP code on web at http://nris.state.mt.us/gis/projection/projection.html. Translated to javascript by Jeff Miller, Jefferson County, WA.
// The formulae this program is based on are from "Map Projections, A Working Manual" by John P. Snyder, U.S. GeoLogical Survey - https://pubs.er.usgs.gov/publication/pp1395
// Professional Paper 1395, 1987, pages 295-298
// Set up the coordinate system parameters. Need to edit these based on state plane zone.

 

function convertSP(uX,uY) {
// uX = lambda, uY = phi

feet = 3.28084;
angRad = 0.01745329252; //number of radians in a degree

 

// Constants - https://spatialreference.org/ref/epsg/2227/
a = 6378137 * feet; //20925604.47416666667; //major radius of ellipsoid, map units (NAD 83)
ec = Math.sqrt((2 / 298.257222101) - (1 / Math.pow(298.257222101, 2))); //0.081819191; //eccentricity of ellipsoid (NAD 83)
uY = uY * angRad;
uX = uX * angRad;
pi4 = Math.PI/4.0 //3.141592653582 / 4; //Pi / 4
phi0 = 36.5 * angRad; //latitude of origin
phi1 = 37.066666666666667 * angRad; //latitude of first standard parallel
phi2 = 38.43333333333 * angRad; //latitude of second standard parallel
lambda0 = -120.5 * angRad; //central meridian
x0 = 6561666.6666667; //False easting of central meridian, map units
y0 = 1640416.6666667; // False northing

 

// Calculate the coordinate system constants. pages 107, 108, 296 in Map Projects, A Working Manual
with (Math) {
m1 = cos(phi1) / sqrt(1 - (pow(ec,2)) * pow(sin(phi1),2));
m2 = cos(phi2) / sqrt(1 - (pow(ec,2)) * pow(sin(phi2),2));

t = tan(pi4 - (uY/2)) / pow(((1 - (ec * sin(uY))) / (1 + (ec * sin(uY)))), ec / 2);
t0 = tan(pi4 - (phi0/2)) / pow(((1 - (ec * sin(phi0))) / (1 + (ec * sin(phi0)))), ec / 2);
t1 = tan(pi4 - (phi1/2)) / pow(((1 - (ec * sin(phi1))) / (1 + (ec * sin(phi1)))), ec / 2);
t2 = tan(pi4 - (phi2/2)) / pow(((1 - (ec * sin(phi2))) / (1 + (ec * sin(phi2)))), ec / 2);

n = (log(m1) - log(m2)) / (log(t1) - log(t2));
f = m1 / (n * pow(t1,n));
rho0 = a * f * pow(t0,n);

rho = a * f * pow(t, n);
theta = n * (uX - lambda0);
easting = (rho * sin(theta)) + x0;
northing = rho0 - (rho * cos(theta)) + y0;

alert(easting + " | " + northing);
}
}

Outcomes