var HALF_PI = Math.PI/2;
|
var PJD_3PARAM = 1;
|
var PJD_7PARAM = 2;
|
var PJD_GRIDSHIFT = 3;
|
var PJD_WGS84 = 4; // WGS84 or equivalent
|
var PJD_NODATUM = 5; // WGS84 or equivalent
|
var SEC_TO_RAD = 4.84813681109535993589914102357e-6;
|
var AD_C = 1.0026000;
|
var COS_67P5 = 0.38268343236508977;
|
var datum = function(proj) {
|
if (!(this instanceof datum)) {
|
return new datum(proj);
|
}
|
this.datum_type = PJD_WGS84; //default setting
|
if (!proj) {
|
return;
|
}
|
if (proj.datumCode && proj.datumCode === 'none') {
|
this.datum_type = PJD_NODATUM;
|
}
|
|
if (proj.datum_params) {
|
this.datum_params = proj.datum_params.map(parseFloat);
|
if (this.datum_params[0] !== 0 || this.datum_params[1] !== 0 || this.datum_params[2] !== 0) {
|
this.datum_type = PJD_3PARAM;
|
}
|
if (this.datum_params.length > 3) {
|
if (this.datum_params[3] !== 0 || this.datum_params[4] !== 0 || this.datum_params[5] !== 0 || this.datum_params[6] !== 0) {
|
this.datum_type = PJD_7PARAM;
|
this.datum_params[3] *= SEC_TO_RAD;
|
this.datum_params[4] *= SEC_TO_RAD;
|
this.datum_params[5] *= SEC_TO_RAD;
|
this.datum_params[6] = (this.datum_params[6] / 1000000.0) + 1.0;
|
}
|
}
|
}
|
|
// DGR 2011-03-21 : nadgrids support
|
this.datum_type = proj.grids ? PJD_GRIDSHIFT : this.datum_type;
|
|
this.a = proj.a; //datum object also uses these values
|
this.b = proj.b;
|
this.es = proj.es;
|
this.ep2 = proj.ep2;
|
if (this.datum_type === PJD_GRIDSHIFT) {
|
this.grids = proj.grids;
|
}
|
};
|
datum.prototype = {
|
|
|
/****************************************************************/
|
// cs_compare_datums()
|
// Returns TRUE if the two datums match, otherwise FALSE.
|
compare_datums: function(dest) {
|
if (this.datum_type !== dest.datum_type) {
|
return false; // false, datums are not equal
|
}
|
else if (this.a !== dest.a || Math.abs(this.es - dest.es) > 0.000000000050) {
|
// the tolerence for es is to ensure that GRS80 and WGS84
|
// are considered identical
|
return false;
|
}
|
else if (this.datum_type === PJD_3PARAM) {
|
return (this.datum_params[0] === dest.datum_params[0] && this.datum_params[1] === dest.datum_params[1] && this.datum_params[2] === dest.datum_params[2]);
|
}
|
else if (this.datum_type === PJD_7PARAM) {
|
return (this.datum_params[0] === dest.datum_params[0] && this.datum_params[1] === dest.datum_params[1] && this.datum_params[2] === dest.datum_params[2] && this.datum_params[3] === dest.datum_params[3] && this.datum_params[4] === dest.datum_params[4] && this.datum_params[5] === dest.datum_params[5] && this.datum_params[6] === dest.datum_params[6]);
|
}
|
else if (this.datum_type === PJD_GRIDSHIFT || dest.datum_type === PJD_GRIDSHIFT) {
|
//alert("ERROR: Grid shift transformations are not implemented.");
|
//return false
|
//DGR 2012-07-29 lazy ...
|
return this.nadgrids === dest.nadgrids;
|
}
|
else {
|
return true; // datums are equal
|
}
|
}, // cs_compare_datums()
|
|
/*
|
* The function Convert_Geodetic_To_Geocentric converts geodetic coordinates
|
* (latitude, longitude, and height) to geocentric coordinates (X, Y, Z),
|
* according to the current ellipsoid parameters.
|
*
|
* Latitude : Geodetic latitude in radians (input)
|
* Longitude : Geodetic longitude in radians (input)
|
* Height : Geodetic height, in meters (input)
|
* X : Calculated Geocentric X coordinate, in meters (output)
|
* Y : Calculated Geocentric Y coordinate, in meters (output)
|
* Z : Calculated Geocentric Z coordinate, in meters (output)
|
*
|
*/
|
geodetic_to_geocentric: function(p) {
|
var Longitude = p.x;
|
var Latitude = p.y;
|
var Height = p.z ? p.z : 0; //Z value not always supplied
|
var X; // output
|
var Y;
|
var Z;
|
|
var Error_Code = 0; // GEOCENT_NO_ERROR;
|
var Rn; /* Earth radius at location */
|
var Sin_Lat; /* Math.sin(Latitude) */
|
var Sin2_Lat; /* Square of Math.sin(Latitude) */
|
var Cos_Lat; /* Math.cos(Latitude) */
|
|
/*
|
** Don't blow up if Latitude is just a little out of the value
|
** range as it may just be a rounding issue. Also removed longitude
|
** test, it should be wrapped by Math.cos() and Math.sin(). NFW for PROJ.4, Sep/2001.
|
*/
|
if (Latitude < -HALF_PI && Latitude > -1.001 * HALF_PI) {
|
Latitude = -HALF_PI;
|
}
|
else if (Latitude > HALF_PI && Latitude < 1.001 * HALF_PI) {
|
Latitude = HALF_PI;
|
}
|
else if ((Latitude < -HALF_PI) || (Latitude > HALF_PI)) {
|
/* Latitude out of range */
|
//..reportError('geocent:lat out of range:' + Latitude);
|
return null;
|
}
|
|
if (Longitude > Math.PI) {
|
Longitude -= (2 * Math.PI);
|
}
|
Sin_Lat = Math.sin(Latitude);
|
Cos_Lat = Math.cos(Latitude);
|
Sin2_Lat = Sin_Lat * Sin_Lat;
|
Rn = this.a / (Math.sqrt(1.0e0 - this.es * Sin2_Lat));
|
X = (Rn + Height) * Cos_Lat * Math.cos(Longitude);
|
Y = (Rn + Height) * Cos_Lat * Math.sin(Longitude);
|
Z = ((Rn * (1 - this.es)) + Height) * Sin_Lat;
|
|
p.x = X;
|
p.y = Y;
|
p.z = Z;
|
return Error_Code;
|
}, // cs_geodetic_to_geocentric()
|
|
|
geocentric_to_geodetic: function(p) {
|
/* local defintions and variables */
|
/* end-criterium of loop, accuracy of sin(Latitude) */
|
var genau = 1e-12;
|
var genau2 = (genau * genau);
|
var maxiter = 30;
|
|
var P; /* distance between semi-minor axis and location */
|
var RR; /* distance between center and location */
|
var CT; /* sin of geocentric latitude */
|
var ST; /* cos of geocentric latitude */
|
var RX;
|
var RK;
|
var RN; /* Earth radius at location */
|
var CPHI0; /* cos of start or old geodetic latitude in iterations */
|
var SPHI0; /* sin of start or old geodetic latitude in iterations */
|
var CPHI; /* cos of searched geodetic latitude */
|
var SPHI; /* sin of searched geodetic latitude */
|
var SDPHI; /* end-criterium: addition-theorem of sin(Latitude(iter)-Latitude(iter-1)) */
|
var At_Pole; /* indicates location is in polar region */
|
var iter; /* # of continous iteration, max. 30 is always enough (s.a.) */
|
|
var X = p.x;
|
var Y = p.y;
|
var Z = p.z ? p.z : 0.0; //Z value not always supplied
|
var Longitude;
|
var Latitude;
|
var Height;
|
|
At_Pole = false;
|
P = Math.sqrt(X * X + Y * Y);
|
RR = Math.sqrt(X * X + Y * Y + Z * Z);
|
|
/* special cases for latitude and longitude */
|
if (P / this.a < genau) {
|
|
/* special case, if P=0. (X=0., Y=0.) */
|
At_Pole = true;
|
Longitude = 0.0;
|
|
/* if (X,Y,Z)=(0.,0.,0.) then Height becomes semi-minor axis
|
* of ellipsoid (=center of mass), Latitude becomes PI/2 */
|
if (RR / this.a < genau) {
|
Latitude = HALF_PI;
|
Height = -this.b;
|
return;
|
}
|
}
|
else {
|
/* ellipsoidal (geodetic) longitude
|
* interval: -PI < Longitude <= +PI */
|
Longitude = Math.atan2(Y, X);
|
}
|
|
/* --------------------------------------------------------------
|
* Following iterative algorithm was developped by
|
* "Institut for Erdmessung", University of Hannover, July 1988.
|
* Internet: www.ife.uni-hannover.de
|
* Iterative computation of CPHI,SPHI and Height.
|
* Iteration of CPHI and SPHI to 10**-12 radian resp.
|
* 2*10**-7 arcsec.
|
* --------------------------------------------------------------
|
*/
|
CT = Z / RR;
|
ST = P / RR;
|
RX = 1.0 / Math.sqrt(1.0 - this.es * (2.0 - this.es) * ST * ST);
|
CPHI0 = ST * (1.0 - this.es) * RX;
|
SPHI0 = CT * RX;
|
iter = 0;
|
|
/* loop to find sin(Latitude) resp. Latitude
|
* until |sin(Latitude(iter)-Latitude(iter-1))| < genau */
|
do {
|
iter++;
|
RN = this.a / Math.sqrt(1.0 - this.es * SPHI0 * SPHI0);
|
|
/* ellipsoidal (geodetic) height */
|
Height = P * CPHI0 + Z * SPHI0 - RN * (1.0 - this.es * SPHI0 * SPHI0);
|
|
RK = this.es * RN / (RN + Height);
|
RX = 1.0 / Math.sqrt(1.0 - RK * (2.0 - RK) * ST * ST);
|
CPHI = ST * (1.0 - RK) * RX;
|
SPHI = CT * RX;
|
SDPHI = SPHI * CPHI0 - CPHI * SPHI0;
|
CPHI0 = CPHI;
|
SPHI0 = SPHI;
|
}
|
while (SDPHI * SDPHI > genau2 && iter < maxiter);
|
|
/* ellipsoidal (geodetic) latitude */
|
Latitude = Math.atan(SPHI / Math.abs(CPHI));
|
|
p.x = Longitude;
|
p.y = Latitude;
|
p.z = Height;
|
return p;
|
}, // cs_geocentric_to_geodetic()
|
|
/** Convert_Geocentric_To_Geodetic
|
* The method used here is derived from 'An Improved Algorithm for
|
* Geocentric to Geodetic Coordinate Conversion', by Ralph Toms, Feb 1996
|
*/
|
geocentric_to_geodetic_noniter: function(p) {
|
var X = p.x;
|
var Y = p.y;
|
var Z = p.z ? p.z : 0; //Z value not always supplied
|
var Longitude;
|
var Latitude;
|
var Height;
|
|
var W; /* distance from Z axis */
|
var W2; /* square of distance from Z axis */
|
var T0; /* initial estimate of vertical component */
|
var T1; /* corrected estimate of vertical component */
|
var S0; /* initial estimate of horizontal component */
|
var S1; /* corrected estimate of horizontal component */
|
var Sin_B0; /* Math.sin(B0), B0 is estimate of Bowring aux variable */
|
var Sin3_B0; /* cube of Math.sin(B0) */
|
var Cos_B0; /* Math.cos(B0) */
|
var Sin_p1; /* Math.sin(phi1), phi1 is estimated latitude */
|
var Cos_p1; /* Math.cos(phi1) */
|
var Rn; /* Earth radius at location */
|
var Sum; /* numerator of Math.cos(phi1) */
|
var At_Pole; /* indicates location is in polar region */
|
|
X = parseFloat(X); // cast from string to float
|
Y = parseFloat(Y);
|
Z = parseFloat(Z);
|
|
At_Pole = false;
|
if (X !== 0.0) {
|
Longitude = Math.atan2(Y, X);
|
}
|
else {
|
if (Y > 0) {
|
Longitude = HALF_PI;
|
}
|
else if (Y < 0) {
|
Longitude = -HALF_PI;
|
}
|
else {
|
At_Pole = true;
|
Longitude = 0.0;
|
if (Z > 0.0) { /* north pole */
|
Latitude = HALF_PI;
|
}
|
else if (Z < 0.0) { /* south pole */
|
Latitude = -HALF_PI;
|
}
|
else { /* center of earth */
|
Latitude = HALF_PI;
|
Height = -this.b;
|
return;
|
}
|
}
|
}
|
W2 = X * X + Y * Y;
|
W = Math.sqrt(W2);
|
T0 = Z * AD_C;
|
S0 = Math.sqrt(T0 * T0 + W2);
|
Sin_B0 = T0 / S0;
|
Cos_B0 = W / S0;
|
Sin3_B0 = Sin_B0 * Sin_B0 * Sin_B0;
|
T1 = Z + this.b * this.ep2 * Sin3_B0;
|
Sum = W - this.a * this.es * Cos_B0 * Cos_B0 * Cos_B0;
|
S1 = Math.sqrt(T1 * T1 + Sum * Sum);
|
Sin_p1 = T1 / S1;
|
Cos_p1 = Sum / S1;
|
Rn = this.a / Math.sqrt(1.0 - this.es * Sin_p1 * Sin_p1);
|
if (Cos_p1 >= COS_67P5) {
|
Height = W / Cos_p1 - Rn;
|
}
|
else if (Cos_p1 <= -COS_67P5) {
|
Height = W / -Cos_p1 - Rn;
|
}
|
else {
|
Height = Z / Sin_p1 + Rn * (this.es - 1.0);
|
}
|
if (At_Pole === false) {
|
Latitude = Math.atan(Sin_p1 / Cos_p1);
|
}
|
|
p.x = Longitude;
|
p.y = Latitude;
|
p.z = Height;
|
return p;
|
}, // geocentric_to_geodetic_noniter()
|
|
/****************************************************************/
|
// pj_geocentic_to_wgs84( p )
|
// p = point to transform in geocentric coordinates (x,y,z)
|
geocentric_to_wgs84: function(p) {
|
|
if (this.datum_type === PJD_3PARAM) {
|
// if( x[io] === HUGE_VAL )
|
// continue;
|
p.x += this.datum_params[0];
|
p.y += this.datum_params[1];
|
p.z += this.datum_params[2];
|
|
}
|
else if (this.datum_type === PJD_7PARAM) {
|
var Dx_BF = this.datum_params[0];
|
var Dy_BF = this.datum_params[1];
|
var Dz_BF = this.datum_params[2];
|
var Rx_BF = this.datum_params[3];
|
var Ry_BF = this.datum_params[4];
|
var Rz_BF = this.datum_params[5];
|
var M_BF = this.datum_params[6];
|
// if( x[io] === HUGE_VAL )
|
// continue;
|
var x_out = M_BF * (p.x - Rz_BF * p.y + Ry_BF * p.z) + Dx_BF;
|
var y_out = M_BF * (Rz_BF * p.x + p.y - Rx_BF * p.z) + Dy_BF;
|
var z_out = M_BF * (-Ry_BF * p.x + Rx_BF * p.y + p.z) + Dz_BF;
|
p.x = x_out;
|
p.y = y_out;
|
p.z = z_out;
|
}
|
}, // cs_geocentric_to_wgs84
|
|
/****************************************************************/
|
// pj_geocentic_from_wgs84()
|
// coordinate system definition,
|
// point to transform in geocentric coordinates (x,y,z)
|
geocentric_from_wgs84: function(p) {
|
|
if (this.datum_type === PJD_3PARAM) {
|
//if( x[io] === HUGE_VAL )
|
// continue;
|
p.x -= this.datum_params[0];
|
p.y -= this.datum_params[1];
|
p.z -= this.datum_params[2];
|
|
}
|
else if (this.datum_type === PJD_7PARAM) {
|
var Dx_BF = this.datum_params[0];
|
var Dy_BF = this.datum_params[1];
|
var Dz_BF = this.datum_params[2];
|
var Rx_BF = this.datum_params[3];
|
var Ry_BF = this.datum_params[4];
|
var Rz_BF = this.datum_params[5];
|
var M_BF = this.datum_params[6];
|
var x_tmp = (p.x - Dx_BF) / M_BF;
|
var y_tmp = (p.y - Dy_BF) / M_BF;
|
var z_tmp = (p.z - Dz_BF) / M_BF;
|
//if( x[io] === HUGE_VAL )
|
// continue;
|
|
p.x = x_tmp + Rz_BF * y_tmp - Ry_BF * z_tmp;
|
p.y = -Rz_BF * x_tmp + y_tmp + Rx_BF * z_tmp;
|
p.z = Ry_BF * x_tmp - Rx_BF * y_tmp + z_tmp;
|
} //cs_geocentric_from_wgs84()
|
}
|
};
|
|
/** point object, nothing fancy, just allows values to be
|
passed back and forth by reference rather than by value.
|
Other point classes may be used as long as they have
|
x and y properties, which will get modified in the transform method.
|
*/
|
module.exports = datum;
|
