/*******************************************************************************
|
NAME LAMBERT AZIMUTHAL EQUAL-AREA
|
|
PURPOSE: Transforms input longitude and latitude to Easting and
|
Northing for the Lambert Azimuthal Equal-Area projection. The
|
longitude and latitude must be in radians. The Easting
|
and Northing values will be returned in meters.
|
|
PROGRAMMER DATE
|
---------- ----
|
D. Steinwand, EROS March, 1991
|
|
This function was adapted from the Lambert Azimuthal Equal Area projection
|
code (FORTRAN) in the General Cartographic Transformation Package software
|
which is available from the U.S. Geological Survey National Mapping Division.
|
|
ALGORITHM REFERENCES
|
|
1. "New Equal-Area Map Projections for Noncircular Regions", John P. Snyder,
|
The American Cartographer, Vol 15, No. 4, October 1988, pp. 341-355.
|
|
2. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
|
Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
|
State Government Printing Office, Washington D.C., 1987.
|
|
3. "Software Documentation for GCTP General Cartographic Transformation
|
Package", U.S. Geological Survey National Mapping Division, May 1982.
|
*******************************************************************************/
|
|
Proj4js.Proj.laea = {
|
S_POLE: 1,
|
N_POLE: 2,
|
EQUIT: 3,
|
OBLIQ: 4,
|
|
|
/* Initialize the Lambert Azimuthal Equal Area projection
|
------------------------------------------------------*/
|
init: function() {
|
var t = Math.abs(this.lat0);
|
if (Math.abs(t - Proj4js.common.HALF_PI) < Proj4js.common.EPSLN) {
|
this.mode = this.lat0 < 0. ? this.S_POLE : this.N_POLE;
|
} else if (Math.abs(t) < Proj4js.common.EPSLN) {
|
this.mode = this.EQUIT;
|
} else {
|
this.mode = this.OBLIQ;
|
}
|
if (this.es > 0) {
|
var sinphi;
|
|
this.qp = Proj4js.common.qsfnz(this.e, 1.0);
|
this.mmf = .5 / (1. - this.es);
|
this.apa = this.authset(this.es);
|
switch (this.mode) {
|
case this.N_POLE:
|
case this.S_POLE:
|
this.dd = 1.;
|
break;
|
case this.EQUIT:
|
this.rq = Math.sqrt(.5 * this.qp);
|
this.dd = 1. / this.rq;
|
this.xmf = 1.;
|
this.ymf = .5 * this.qp;
|
break;
|
case this.OBLIQ:
|
this.rq = Math.sqrt(.5 * this.qp);
|
sinphi = Math.sin(this.lat0);
|
this.sinb1 = Proj4js.common.qsfnz(this.e, sinphi) / this.qp;
|
this.cosb1 = Math.sqrt(1. - this.sinb1 * this.sinb1);
|
this.dd = Math.cos(this.lat0) / (Math.sqrt(1. - this.es * sinphi * sinphi) * this.rq * this.cosb1);
|
this.ymf = (this.xmf = this.rq) / this.dd;
|
this.xmf *= this.dd;
|
break;
|
}
|
} else {
|
if (this.mode == this.OBLIQ) {
|
this.sinph0 = Math.sin(this.lat0);
|
this.cosph0 = Math.cos(this.lat0);
|
}
|
}
|
},
|
|
/* Lambert Azimuthal Equal Area forward equations--mapping lat,long to x,y
|
-----------------------------------------------------------------------*/
|
forward: function(p) {
|
|
/* Forward equations
|
-----------------*/
|
var x,y;
|
var lam=p.x;
|
var phi=p.y;
|
lam = Proj4js.common.adjust_lon(lam - this.long0);
|
|
if (this.sphere) {
|
var coslam, cosphi, sinphi;
|
|
sinphi = Math.sin(phi);
|
cosphi = Math.cos(phi);
|
coslam = Math.cos(lam);
|
switch (this.mode) {
|
case this.OBLIQ:
|
case this.EQUIT:
|
y = (this.mode == this.EQUIT) ? 1. + cosphi * coslam : 1. + this.sinph0 * sinphi + this.cosph0 * cosphi * coslam;
|
if (y <= Proj4js.common.EPSLN) {
|
Proj4js.reportError("laea:fwd:y less than eps");
|
return null;
|
}
|
y = Math.sqrt(2. / y);
|
x = y * cosphi * Math.sin(lam);
|
y *= (this.mode == this.EQUIT) ? sinphi : this.cosph0 * sinphi - this.sinph0 * cosphi * coslam;
|
break;
|
case this.N_POLE:
|
coslam = -coslam;
|
case this.S_POLE:
|
if (Math.abs(phi + this.phi0) < Proj4js.common.EPSLN) {
|
Proj4js.reportError("laea:fwd:phi < eps");
|
return null;
|
}
|
y = Proj4js.common.FORTPI - phi * .5;
|
y = 2. * ((this.mode == this.S_POLE) ? Math.cos(y) : Math.sin(y));
|
x = y * Math.sin(lam);
|
y *= coslam;
|
break;
|
}
|
} else {
|
var coslam, sinlam, sinphi, q, sinb=0.0, cosb=0.0, b=0.0;
|
|
coslam = Math.cos(lam);
|
sinlam = Math.sin(lam);
|
sinphi = Math.sin(phi);
|
q = Proj4js.common.qsfnz(this.e, sinphi);
|
if (this.mode == this.OBLIQ || this.mode == this.EQUIT) {
|
sinb = q / this.qp;
|
cosb = Math.sqrt(1. - sinb * sinb);
|
}
|
switch (this.mode) {
|
case this.OBLIQ:
|
b = 1. + this.sinb1 * sinb + this.cosb1 * cosb * coslam;
|
break;
|
case this.EQUIT:
|
b = 1. + cosb * coslam;
|
break;
|
case this.N_POLE:
|
b = Proj4js.common.HALF_PI + phi;
|
q = this.qp - q;
|
break;
|
case this.S_POLE:
|
b = phi - Proj4js.common.HALF_PI;
|
q = this.qp + q;
|
break;
|
}
|
if (Math.abs(b) < Proj4js.common.EPSLN) {
|
Proj4js.reportError("laea:fwd:b < eps");
|
return null;
|
}
|
switch (this.mode) {
|
case this.OBLIQ:
|
case this.EQUIT:
|
b = Math.sqrt(2. / b);
|
if (this.mode == this.OBLIQ) {
|
y = this.ymf * b * (this.cosb1 * sinb - this.sinb1 * cosb * coslam);
|
} else {
|
y = (b = Math.sqrt(2. / (1. + cosb * coslam))) * sinb * this.ymf;
|
}
|
x = this.xmf * b * cosb * sinlam;
|
break;
|
case this.N_POLE:
|
case this.S_POLE:
|
if (q >= 0.) {
|
x = (b = Math.sqrt(q)) * sinlam;
|
y = coslam * ((this.mode == this.S_POLE) ? b : -b);
|
} else {
|
x = y = 0.;
|
}
|
break;
|
}
|
}
|
|
//v 1.0
|
/*
|
var sin_lat=Math.sin(lat);
|
var cos_lat=Math.cos(lat);
|
|
var sin_delta_lon=Math.sin(delta_lon);
|
var cos_delta_lon=Math.cos(delta_lon);
|
|
var g =this.sin_lat_o * sin_lat +this.cos_lat_o * cos_lat * cos_delta_lon;
|
if (g == -1.0) {
|
Proj4js.reportError("laea:fwd:Point projects to a circle of radius "+ 2.0 * R);
|
return null;
|
}
|
var ksp = this.a * Math.sqrt(2.0 / (1.0 + g));
|
var x = ksp * cos_lat * sin_delta_lon + this.x0;
|
var y = ksp * (this.cos_lat_o * sin_lat - this.sin_lat_o * cos_lat * cos_delta_lon) + this.y0;
|
*/
|
p.x = this.a*x + this.x0;
|
p.y = this.a*y + this.y0;
|
return p;
|
},//lamazFwd()
|
|
/* Inverse equations
|
-----------------*/
|
inverse: function(p) {
|
p.x -= this.x0;
|
p.y -= this.y0;
|
var x = p.x/this.a;
|
var y = p.y/this.a;
|
var lam, phi;
|
|
if (this.sphere) {
|
var cosz=0.0, rh, sinz=0.0;
|
|
rh = Math.sqrt(x*x + y*y);
|
phi = rh * .5;
|
if (phi > 1.) {
|
Proj4js.reportError("laea:Inv:DataError");
|
return null;
|
}
|
phi = 2. * Math.asin(phi);
|
if (this.mode == this.OBLIQ || this.mode == this.EQUIT) {
|
sinz = Math.sin(phi);
|
cosz = Math.cos(phi);
|
}
|
switch (this.mode) {
|
case this.EQUIT:
|
phi = (Math.abs(rh) <= Proj4js.common.EPSLN) ? 0. : Math.asin(y * sinz / rh);
|
x *= sinz;
|
y = cosz * rh;
|
break;
|
case this.OBLIQ:
|
phi = (Math.abs(rh) <= Proj4js.common.EPSLN) ? this.phi0 : Math.asin(cosz * this.sinph0 + y * sinz * this.cosph0 / rh);
|
x *= sinz * this.cosph0;
|
y = (cosz - Math.sin(phi) * this.sinph0) * rh;
|
break;
|
case this.N_POLE:
|
y = -y;
|
phi = Proj4js.common.HALF_PI - phi;
|
break;
|
case this.S_POLE:
|
phi -= Proj4js.common.HALF_PI;
|
break;
|
}
|
lam = (y == 0. && (this.mode == this.EQUIT || this.mode == this.OBLIQ)) ? 0. : Math.atan2(x, y);
|
} else {
|
var cCe, sCe, q, rho, ab=0.0;
|
|
switch (this.mode) {
|
case this.EQUIT:
|
case this.OBLIQ:
|
x /= this.dd;
|
y *= this.dd;
|
rho = Math.sqrt(x*x + y*y);
|
if (rho < Proj4js.common.EPSLN) {
|
p.x = 0.;
|
p.y = this.phi0;
|
return p;
|
}
|
sCe = 2. * Math.asin(.5 * rho / this.rq);
|
cCe = Math.cos(sCe);
|
x *= (sCe = Math.sin(sCe));
|
if (this.mode == this.OBLIQ) {
|
ab = cCe * this.sinb1 + y * sCe * this.cosb1 / rho
|
q = this.qp * ab;
|
y = rho * this.cosb1 * cCe - y * this.sinb1 * sCe;
|
} else {
|
ab = y * sCe / rho;
|
q = this.qp * ab;
|
y = rho * cCe;
|
}
|
break;
|
case this.N_POLE:
|
y = -y;
|
case this.S_POLE:
|
q = (x * x + y * y);
|
if (!q ) {
|
p.x = 0.;
|
p.y = this.phi0;
|
return p;
|
}
|
/*
|
q = this.qp - q;
|
*/
|
ab = 1. - q / this.qp;
|
if (this.mode == this.S_POLE) {
|
ab = - ab;
|
}
|
break;
|
}
|
lam = Math.atan2(x, y);
|
phi = this.authlat(Math.asin(ab), this.apa);
|
}
|
|
/*
|
var Rh = Math.Math.sqrt(p.x *p.x +p.y * p.y);
|
var temp = Rh / (2.0 * this.a);
|
|
if (temp > 1) {
|
Proj4js.reportError("laea:Inv:DataError");
|
return null;
|
}
|
|
var z = 2.0 * Proj4js.common.asinz(temp);
|
var sin_z=Math.sin(z);
|
var cos_z=Math.cos(z);
|
|
var lon =this.long0;
|
if (Math.abs(Rh) > Proj4js.common.EPSLN) {
|
var lat = Proj4js.common.asinz(this.sin_lat_o * cos_z +this. cos_lat_o * sin_z *p.y / Rh);
|
var temp =Math.abs(this.lat0) - Proj4js.common.HALF_PI;
|
if (Math.abs(temp) > Proj4js.common.EPSLN) {
|
temp = cos_z -this.sin_lat_o * Math.sin(lat);
|
if(temp!=0.0) lon=Proj4js.common.adjust_lon(this.long0+Math.atan2(p.x*sin_z*this.cos_lat_o,temp*Rh));
|
} else if (this.lat0 < 0.0) {
|
lon = Proj4js.common.adjust_lon(this.long0 - Math.atan2(-p.x,p.y));
|
} else {
|
lon = Proj4js.common.adjust_lon(this.long0 + Math.atan2(p.x, -p.y));
|
}
|
} else {
|
lat = this.lat0;
|
}
|
*/
|
//return(OK);
|
p.x = Proj4js.common.adjust_lon(this.long0+lam);
|
p.y = phi;
|
return p;
|
},//lamazInv()
|
|
/* determine latitude from authalic latitude */
|
P00: .33333333333333333333,
|
P01: .17222222222222222222,
|
P02: .10257936507936507936,
|
P10: .06388888888888888888,
|
P11: .06640211640211640211,
|
P20: .01641501294219154443,
|
|
authset: function(es) {
|
var t;
|
var APA = new Array();
|
APA[0] = es * this.P00;
|
t = es * es;
|
APA[0] += t * this.P01;
|
APA[1] = t * this.P10;
|
t *= es;
|
APA[0] += t * this.P02;
|
APA[1] += t * this.P11;
|
APA[2] = t * this.P20;
|
return APA;
|
},
|
|
authlat: function(beta, APA) {
|
var t = beta+beta;
|
return(beta + APA[0] * Math.sin(t) + APA[1] * Math.sin(t+t) + APA[2] * Math.sin(t+t+t));
|
}
|
|
};
|
