mirror of
https://github.com/proj4js/proj4js.git
synced 2026-01-25 16:26:11 +00:00
478 lines
13 KiB
JavaScript
478 lines
13 KiB
JavaScript
exports.PI = 3.141592653589793238; //Math.PI;
|
|
exports.HALF_PI = 1.570796326794896619; //Math.PI*0.5;
|
|
exports.TWO_PI = 6.283185307179586477; //Math.PI*2;
|
|
exports.FORTPI = 0.78539816339744833;
|
|
exports.R2D = 57.29577951308232088;
|
|
exports.D2R = 0.01745329251994329577;
|
|
exports.SEC_TO_RAD = 4.84813681109535993589914102357e-6;
|
|
/* SEC_TO_RAD = Pi/180/3600 */
|
|
exports.EPSLN = 1.0e-10;
|
|
exports.MAX_ITER = 20;
|
|
// following constants from geocent.c
|
|
exports.COS_67P5 = 0.38268343236508977;
|
|
/* cosine of 67.5 degrees */
|
|
exports.AD_C = 1.0026000;
|
|
/* Toms region 1 constant */
|
|
|
|
/* datum_type values */
|
|
exports.PJD_UNKNOWN = 0;
|
|
exports.PJD_3PARAM = 1;
|
|
exports.PJD_7PARAM = 2;
|
|
exports.PJD_GRIDSHIFT = 3;
|
|
exports.PJD_WGS84 = 4; // WGS84 or equivalent
|
|
exports.PJD_NODATUM = 5; // WGS84 or equivalent
|
|
exports.SRS_WGS84_SEMIMAJOR = 6378137; // only used in grid shift transforms
|
|
exports.SRS_WGS84_ESQUARED = 0.006694379990141316; //DGR: 2012-07-29
|
|
|
|
// ellipoid pj_set_ell.c
|
|
exports.SIXTH = 0.1666666666666666667;
|
|
/* 1/6 */
|
|
exports.RA4 = 0.04722222222222222222;
|
|
/* 17/360 */
|
|
exports.RA6 = 0.02215608465608465608;
|
|
/* 67/3024 */
|
|
exports.RV4 = 0.06944444444444444444;
|
|
/* 5/72 */
|
|
exports.RV6 = 0.04243827160493827160;
|
|
/* 55/1296 */
|
|
|
|
// Function to compute the constant small m which is the radius of
|
|
// a parallel of latitude; phi; divided by the semimajor axis.
|
|
// -----------------------------------------------------------------
|
|
exports.msfnz = function(eccent, sinphi, cosphi) {
|
|
var con = eccent * sinphi;
|
|
return cosphi / (Math.sqrt(1 - con * con));
|
|
};
|
|
|
|
// Function to compute the constant small t for use in the forward
|
|
// computations in the Lambert Conformal Conic and the Polar
|
|
// Stereographic projections.
|
|
// -----------------------------------------------------------------
|
|
exports.tsfnz = function(eccent, phi, sinphi) {
|
|
var con = eccent * sinphi;
|
|
var com = 0.5 * eccent;
|
|
con = Math.pow(((1 - con) / (1 + con)), com);
|
|
return (Math.tan(0.5 * (this.HALF_PI - phi)) / con);
|
|
};
|
|
|
|
// Function to compute the latitude angle; phi2; for the inverse of the
|
|
// Lambert Conformal Conic and Polar Stereographic projections.
|
|
// ----------------------------------------------------------------
|
|
exports.phi2z = function(eccent, ts) {
|
|
var eccnth = 0.5 * eccent;
|
|
var con, dphi;
|
|
var phi = this.HALF_PI - 2 * Math.atan(ts);
|
|
for (var i = 0; i <= 15; i++) {
|
|
con = eccent * Math.sin(phi);
|
|
dphi = this.HALF_PI - 2 * Math.atan(ts * (Math.pow(((1 - con) / (1 + con)), eccnth))) - phi;
|
|
phi += dphi;
|
|
if (Math.abs(dphi) <= 0.0000000001) {
|
|
return phi;
|
|
}
|
|
}
|
|
//console.log("phi2z has NoConvergence");
|
|
return -9999;
|
|
};
|
|
|
|
/* Function to compute constant small q which is the radius of a
|
|
parallel of latitude, phi, divided by the semimajor axis.
|
|
------------------------------------------------------------*/
|
|
exports.qsfnz = function(eccent, sinphi) {
|
|
var con;
|
|
if (eccent > 1.0e-7) {
|
|
con = eccent * sinphi;
|
|
return ((1 - eccent * eccent) * (sinphi / (1 - con * con) - (0.5 / eccent) * Math.log((1 - con) / (1 + con))));
|
|
}
|
|
else {
|
|
return (2 * sinphi);
|
|
}
|
|
};
|
|
|
|
/* Function to compute the inverse of qsfnz
|
|
------------------------------------------------------------*/
|
|
exports.iqsfnz = function(eccent, q) {
|
|
var temp = 1 - (1 - eccent * eccent) / (2 * eccent) * Math.log((1 - eccent) / (1 + eccent));
|
|
if (Math.abs(Math.abs(q) - temp) < 1.0E-6) {
|
|
if (q < 0) {
|
|
return (-1 * exports.HALF_PI);
|
|
}
|
|
else {
|
|
return exports.HALF_PI;
|
|
}
|
|
}
|
|
//var phi = 0.5* q/(1-eccent*eccent);
|
|
var phi = Math.asin(0.5 * q);
|
|
var dphi;
|
|
var sin_phi;
|
|
var cos_phi;
|
|
var con;
|
|
for (var i = 0; i < 30; i++) {
|
|
sin_phi = Math.sin(phi);
|
|
cos_phi = Math.cos(phi);
|
|
con = eccent * sin_phi;
|
|
dphi = Math.pow(1 - con * con, 2) / (2 * cos_phi) * (q / (1 - eccent * eccent) - sin_phi / (1 - con * con) + 0.5 / eccent * Math.log((1 - con) / (1 + con)));
|
|
phi += dphi;
|
|
if (Math.abs(dphi) <= 0.0000000001) {
|
|
return phi;
|
|
}
|
|
}
|
|
|
|
//console.log("IQSFN-CONV:Latitude failed to converge after 30 iterations");
|
|
return NaN;
|
|
};
|
|
|
|
/* Function to eliminate roundoff errors in asin
|
|
----------------------------------------------*/
|
|
exports.asinz = function(x) {
|
|
if (Math.abs(x) > 1) {
|
|
x = (x > 1) ? 1 : -1;
|
|
}
|
|
return Math.asin(x);
|
|
};
|
|
|
|
// following functions from gctpc cproj.c for transverse mercator projections
|
|
exports.e0fn = function(x) {
|
|
return (1 - 0.25 * x * (1 + x / 16 * (3 + 1.25 * x)));
|
|
};
|
|
exports.e1fn = function(x) {
|
|
return (0.375 * x * (1 + 0.25 * x * (1 + 0.46875 * x)));
|
|
};
|
|
exports.e2fn = function(x) {
|
|
return (0.05859375 * x * x * (1 + 0.75 * x));
|
|
};
|
|
exports.e3fn = function(x) {
|
|
return (x * x * x * (35 / 3072));
|
|
};
|
|
exports.mlfn = function(e0, e1, e2, e3, phi) {
|
|
return (e0 * phi - e1 * Math.sin(2 * phi) + e2 * Math.sin(4 * phi) - e3 * Math.sin(6 * phi));
|
|
};
|
|
exports.imlfn = function(ml, e0, e1, e2, e3) {
|
|
var phi;
|
|
var dphi;
|
|
|
|
phi = ml / e0;
|
|
for (var i = 0; i < 15; i++) {
|
|
dphi = (ml - (e0 * phi - e1 * Math.sin(2 * phi) + e2 * Math.sin(4 * phi) - e3 * Math.sin(6 * phi))) / (e0 - 2 * e1 * Math.cos(2 * phi) + 4 * e2 * Math.cos(4 * phi) - 6 * e3 * Math.cos(6 * phi));
|
|
phi += dphi;
|
|
if (Math.abs(dphi) <= 0.0000000001) {
|
|
return phi;
|
|
}
|
|
}
|
|
|
|
//..reportError("IMLFN-CONV:Latitude failed to converge after 15 iterations");
|
|
return NaN;
|
|
};
|
|
|
|
exports.srat = function(esinp, exp) {
|
|
return (Math.pow((1 - esinp) / (1 + esinp), exp));
|
|
};
|
|
|
|
// Function to return the sign of an argument
|
|
exports.sign = function(x) {
|
|
if (x < 0) {
|
|
return (-1);
|
|
}
|
|
else {
|
|
return (1);
|
|
}
|
|
};
|
|
|
|
// Function to adjust longitude to -180 to 180; input in radians
|
|
exports.adjust_lon = function(x) {
|
|
x = (Math.abs(x) < this.PI) ? x : (x - (this.sign(x) * this.TWO_PI));
|
|
return x;
|
|
};
|
|
|
|
// IGNF - DGR : algorithms used by IGN France
|
|
|
|
// Function to adjust latitude to -90 to 90; input in radians
|
|
exports.adjust_lat = function(x) {
|
|
x = (Math.abs(x) < this.HALF_PI) ? x : (x - (this.sign(x) * this.PI));
|
|
return x;
|
|
};
|
|
|
|
// Latitude Isometrique - close to tsfnz ...
|
|
exports.latiso = function(eccent, phi, sinphi) {
|
|
if (Math.abs(phi) > this.HALF_PI) {
|
|
return Number.NaN;
|
|
}
|
|
if (phi === this.HALF_PI) {
|
|
return Number.POSITIVE_INFINITY;
|
|
}
|
|
if (phi === -1 * this.HALF_PI) {
|
|
return Number.NEGATIVE_INFINITY;
|
|
}
|
|
|
|
var con = eccent * sinphi;
|
|
return Math.log(Math.tan((this.HALF_PI + phi) / 2)) + eccent * Math.log((1 - con) / (1 + con)) / 2;
|
|
};
|
|
|
|
exports.fL = function(x, L) {
|
|
return 2 * Math.atan(x * Math.exp(L)) - this.HALF_PI;
|
|
};
|
|
|
|
// Inverse Latitude Isometrique - close to ph2z
|
|
exports.invlatiso = function(eccent, ts) {
|
|
var phi = this.fL(1, ts);
|
|
var Iphi = 0;
|
|
var con = 0;
|
|
do {
|
|
Iphi = phi;
|
|
con = eccent * Math.sin(Iphi);
|
|
phi = this.fL(Math.exp(eccent * Math.log((1 + con) / (1 - con)) / 2), ts);
|
|
} while (Math.abs(phi - Iphi) > 1.0e-12);
|
|
return phi;
|
|
};
|
|
|
|
// Needed for Gauss Schreiber
|
|
// Original: Denis Makarov (info@binarythings.com)
|
|
// Web Site: http://www.binarythings.com
|
|
exports.sinh = function(x) {
|
|
var r = Math.exp(x);
|
|
r = (r - 1 / r) / 2;
|
|
return r;
|
|
};
|
|
|
|
exports.cosh = function(x) {
|
|
var r = Math.exp(x);
|
|
r = (r + 1 / r) / 2;
|
|
return r;
|
|
};
|
|
|
|
exports.tanh = function(x) {
|
|
var r = Math.exp(x);
|
|
r = (r - 1 / r) / (r + 1 / r);
|
|
return r;
|
|
};
|
|
|
|
exports.asinh = function(x) {
|
|
var s = (x >= 0 ? 1 : -1);
|
|
return s * (Math.log(Math.abs(x) + Math.sqrt(x * x + 1)));
|
|
};
|
|
|
|
exports.acosh = function(x) {
|
|
return 2 * Math.log(Math.sqrt((x + 1) / 2) + Math.sqrt((x - 1) / 2));
|
|
};
|
|
|
|
exports.atanh = function(x) {
|
|
return Math.log((x - 1) / (x + 1)) / 2;
|
|
};
|
|
|
|
// Grande Normale
|
|
exports.gN = function(a, e, sinphi) {
|
|
var temp = e * sinphi;
|
|
return a / Math.sqrt(1 - temp * temp);
|
|
};
|
|
|
|
//code from the PROJ.4 pj_mlfn.c file; this may be useful for other projections
|
|
exports.pj_enfn = function(es) {
|
|
var en = [];
|
|
en[0] = this.C00 - es * (this.C02 + es * (this.C04 + es * (this.C06 + es * this.C08)));
|
|
en[1] = es * (this.C22 - es * (this.C04 + es * (this.C06 + es * this.C08)));
|
|
var t = es * es;
|
|
en[2] = t * (this.C44 - es * (this.C46 + es * this.C48));
|
|
t *= es;
|
|
en[3] = t * (this.C66 - es * this.C68);
|
|
en[4] = t * es * this.C88;
|
|
return en;
|
|
};
|
|
|
|
exports.pj_mlfn = function(phi, sphi, cphi, en) {
|
|
cphi *= sphi;
|
|
sphi *= sphi;
|
|
return (en[0] * phi - cphi * (en[1] + sphi * (en[2] + sphi * (en[3] + sphi * en[4]))));
|
|
};
|
|
|
|
exports.pj_inv_mlfn = function(arg, es, en) {
|
|
var k = 1 / (1 - es);
|
|
var phi = arg;
|
|
for (var i = exports.MAX_ITER; i; --i) { /* rarely goes over 2 iterations */
|
|
var s = Math.sin(phi);
|
|
var t = 1 - es * s * s;
|
|
//t = this.pj_mlfn(phi, s, Math.cos(phi), en) - arg;
|
|
//phi -= t * (t * Math.sqrt(t)) * k;
|
|
t = (this.pj_mlfn(phi, s, Math.cos(phi), en) - arg) * (t * Math.sqrt(t)) * k;
|
|
phi -= t;
|
|
if (Math.abs(t) < exports.EPSLN) {
|
|
return phi;
|
|
}
|
|
}
|
|
//..reportError("cass:pj_inv_mlfn: Convergence error");
|
|
return phi;
|
|
};
|
|
exports.nadInterBreakout = function(indx, frct, letter, number, ct) {
|
|
var inx;
|
|
if (indx[letter] < 0) {
|
|
if (!(indx[letter] === -1 && frct[letter] > 0.99999999999)) {
|
|
return false;
|
|
}
|
|
indx[letter]++;
|
|
frct[letter] = 0;
|
|
}
|
|
else {
|
|
inx = indx[letter] + 1;
|
|
if (inx >= ct.lim[number]) {
|
|
if (!(inx === ct.lim[number] && frct[letter] < 1e-11)) {
|
|
return false;
|
|
}
|
|
if (letter === 'x') {
|
|
indx[letter]--;
|
|
}
|
|
else {
|
|
indx[letter]++;
|
|
}
|
|
frct[letter] = 1;
|
|
}
|
|
}
|
|
return [indx, frct];
|
|
};
|
|
/**
|
|
* Determine correction values
|
|
* source: nad_intr.c (DGR: 2012-07-29)
|
|
*/
|
|
exports.nad_intr = function(pin, ct) {
|
|
// force computation by decreasing by 1e-7 to be as closed as possible
|
|
// from computation under C:C++ by leveraging rounding problems ...
|
|
var t = {
|
|
x: (pin.x - 1e-7) / ct.del[0],
|
|
y: (pin.y - 1e-7) / ct.del[1]
|
|
};
|
|
var indx = {
|
|
x: Math.floor(t.x),
|
|
y: Math.floor(t.y)
|
|
};
|
|
var frct = {
|
|
x: t.x - 1 * indx.x,
|
|
y: t.y - 1 * indx.y
|
|
};
|
|
var val = {
|
|
x: Number.NaN,
|
|
y: Number.NaN
|
|
};
|
|
|
|
|
|
var temp = exports.nadInterBreakout(indx, frct, 'x', 0, ct);
|
|
if (temp) {
|
|
indx = temp[0];
|
|
frct = temp[1];
|
|
}
|
|
else {
|
|
return val;
|
|
}
|
|
temp = exports.nadInterBreakout(indx, frct, 'y', 1, ct);
|
|
if (temp) {
|
|
indx = temp[0];
|
|
frct = temp[1];
|
|
}
|
|
else {
|
|
return val;
|
|
}
|
|
var inx = (indx.y * ct.lim[0]) + indx.x;
|
|
var f00 = {
|
|
x: ct.cvs[inx][0],
|
|
y: ct.cvs[inx][1]
|
|
};
|
|
inx++;
|
|
var f10 = {
|
|
x: ct.cvs[inx][0],
|
|
y: ct.cvs[inx][1]
|
|
};
|
|
inx += ct.lim[0];
|
|
var f11 = {
|
|
x: ct.cvs[inx][0],
|
|
y: ct.cvs[inx][1]
|
|
};
|
|
inx--;
|
|
var f01 = {
|
|
x: ct.cvs[inx][0],
|
|
y: ct.cvs[inx][1]
|
|
};
|
|
var m11 = frct.x * frct.y,
|
|
m10 = frct.x * (1 - frct.y),
|
|
m00 = (1 - frct.x) * (1 - frct.y),
|
|
m01 = (1 - frct.x) * frct.y;
|
|
val.x = (m00 * f00.x + m10 * f10.x + m01 * f01.x + m11 * f11.x);
|
|
val.y = (m00 * f00.y + m10 * f10.y + m01 * f01.y + m11 * f11.y);
|
|
return val;
|
|
};
|
|
|
|
/**
|
|
* Correct value
|
|
* source: nad_cvt.c (DGR: 2012-07-29)
|
|
*/
|
|
exports.inverseNadCvt = function(t, val, tb, ct) {
|
|
if (isNaN(t.x)) {
|
|
return val;
|
|
}
|
|
t.x = tb.x + t.x;
|
|
t.y = tb.y - t.y;
|
|
var i = 9,
|
|
tol = 1e-12;
|
|
var dif, del;
|
|
do {
|
|
del = exports.nad_intr(t, ct);
|
|
if (isNaN(del.x)) {
|
|
this.reportError("Inverse grid shift iteration failed, presumably at grid edge. Using first approximation.");
|
|
break;
|
|
}
|
|
dif = {
|
|
"x": t.x - del.x - tb.x,
|
|
"y": t.y + del.y - tb.y
|
|
};
|
|
t.x -= dif.x;
|
|
t.y -= dif.y;
|
|
} while (i-- && Math.abs(dif.x) > tol && Math.abs(dif.y) > tol);
|
|
if (i < 0) {
|
|
this.reportError("Inverse grid shift iterator failed to converge.");
|
|
return val;
|
|
}
|
|
val.x = exports.adjust_lon(t.x + ct.ll[0]);
|
|
val.y = t.y + ct.ll[1];
|
|
return val;
|
|
};
|
|
exports.nad_cvt = function(pin, inverse, ct) {
|
|
var val = {
|
|
"x": Number.NaN,
|
|
"y": Number.NaN
|
|
};
|
|
if (isNaN(pin.x)) {
|
|
return val;
|
|
}
|
|
var tb = {
|
|
"x": pin.x,
|
|
"y": pin.y
|
|
};
|
|
tb.x -= ct.ll[0];
|
|
tb.y -= ct.ll[1];
|
|
tb.x = exports.adjust_lon(tb.x - exports.PI) + exports.PI;
|
|
var t = exports.nad_intr(tb, ct);
|
|
if (inverse) {
|
|
return exports.inverseNadCvt(t, val, tb, ct);
|
|
}
|
|
else {
|
|
if (!isNaN(t.x)) {
|
|
val.x = pin.x - t.x;
|
|
val.y = pin.y + t.y;
|
|
}
|
|
}
|
|
return val;
|
|
};
|
|
|
|
/* meridinal distance for ellipsoid and inverse
|
|
** 8th degree - accurate to < 1e-5 meters when used in conjuction
|
|
** with typical major axis values.
|
|
** Inverse determines phi to EPS (1e-11) radians, about 1e-6 seconds.
|
|
*/
|
|
exports.C00 = 1;
|
|
exports.C02 = 0.25;
|
|
exports.C04 = 0.046875;
|
|
exports.C06 = 0.01953125;
|
|
exports.C08 = 0.01068115234375;
|
|
exports.C22 = 0.75;
|
|
exports.C44 = 0.46875;
|
|
exports.C46 = 0.01302083333333333333;
|
|
exports.C48 = 0.00712076822916666666;
|
|
exports.C66 = 0.36458333333333333333;
|
|
exports.C68 = 0.00569661458333333333;
|
|
exports.C88 = 0.3076171875;
|