closes #34: add support of proj=laea (76dd0cda) · Commits · REDMIC Project / Third party / Proj4js

lib/proj4js.js

+2 −3

Original line number	Diff line number	Diff line
		@@ -843,7 +843,7 @@ Proj4js.common = {
		/* Function to compute constant small q which is the radius of a
		parallel of latitude, phi, divided by the semimajor axis.
		------------------------------------------------------------*/
		qsfnz : function(eccent,sinphi,cosphi) {
		qsfnz : function(eccent,sinphi) {
		var con;
		if (eccent > 1.0e-7) {
		con = eccent * sinphi;
		@@ -891,8 +891,7 @@ Proj4js.common = {
		},

		// Latitude Isometrique - close to tsfnz ...
		latiso : function(eccent, phi, sinphi)
		{
		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.0this.HALF_PI) return -1.0Number.POSITIVE_INFINITY;

lib/projCode/laea.js

+256 −11

Original line number	Diff line number	Diff line
		@@ -28,13 +28,56 @@ ALGORITHM REFERENCES
		*******************************************************************************/

		Proj4js.Proj.laea = {
		S_POLE: 1,
		N_POLE: 2,
		EQUIT: 3,
		OBLIQ: 4,


		/* Initialize the Lambert Azimuthal Equal Area projection
		------------------------------------------------------*/
		init: function() {
		this.sin_lat_o=Math.sin(this.lat0);
		this.cos_lat_o=Math.cos(this.lat0);
		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
		@@ -43,11 +86,97 @@ Proj4js.Proj.laea = {

		/* Forward equations
		-----------------*/
		var lon=p.x;
		var lat=p.y;
		var delta_lon = Proj4js.common.adjust_lon(lon - this.long0);
		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.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);

		@@ -62,8 +191,9 @@ Proj4js.Proj.laea = {
		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 = x;
		p.y = y;
		*/
		p.x = this.a*x + this.x0;
		p.y = this.a*y + this.y0;
		return p;
		},//lamazFwd()

		@@ -72,8 +202,94 @@ Proj4js.Proj.laea = {
		inverse: function(p) {
		p.x -= this.x0;
		p.y -= this.y0;
		var x = p.x/this.a;
		var y = p.y/this.a;

		if (this.sphere) {
		var cosz=0.0, rh, sinz=0.0;

		var Rh = Math.sqrt(p.x p.x +p.y p.y);
		rh = Math.sqrt(xx + yy);
		var 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 * sinph0 + y * sinz * cosph0 / rh);
		x = sinz cosph0;
		y = (cosz - Math.sin(phi) * 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. : 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(xx + yy);
		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) {
		@@ -100,11 +316,40 @@ Proj4js.Proj.laea = {
		} else {
		lat = this.lat0;
		}
		*/
		//return(OK);
		p.x = lon;
		p.y = lat;
		p.x = Proj4js.common.adjust_lon(this.long0+lam);
		p.y = phi;
		return p;
		}//lamazInv()
		},//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));
		}

		};

test/testdata.js

+54 −30

Original line number	Diff line number	Diff line
		// a set of points in map XY and Lon/Lat that are supposed to correspond between
		// forward and invers transforms
		Proj4js.defs["EPSG:54003"] = "+proj=mill +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +R_A +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
		Proj4js.defs["EPSG:54029"] = "+proj=vandg +lon_0=0 +x_0=0 +y_0=0 +R_A +ellps=WGS84 +datum=WGS84 +units=m +no_defs";
		Proj4js.defs["EPSG:2303X"] = "+proj=utm +zone=30 +ellps=intl +units=m +towgs84=-157.89,-17.16,-78.41,2.118,2.697,-1.434,-1.1097046576093785 +no_defs ";
		Proj4js.defs["EPSG:3035"] = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs";

		Proj4js.testPoints = [
		{code: 'EPSG:3035',
		xy: [4388138.60, 3321736.46],
		ll: [11.0, 53.0]
		},
		{code: 'EPSG:23030',
		xy: [168035.13,4199884.83,-216.62],
		ll: [-6.77432123185356, 37.88456231505968]
		},
		{code: 'EPSG:29100',
		xy: [5110899.06,10552971.81,-22.99],
		ll: [-53.0, 5.0,0.0]
		},
		{code: 'EPSG:27700',
		xy: [343733.14, 612144.53, -51.89],
		ll: [-2.89, 55.4, 0]
		},
		{code: 'EPSG:27492',
		xy: [25260.493584, -9579.245052],
		ll: [-7.84, 39.58]
		},
		{code: 'EPSG:3411',
		xy: [1070076.44,-4635010.27,-136.63],
		ll: [-32, 48, 0]
		},
		{code: 'EPSG:2403',
		xy: [27500000.00, 4198690.08, -109.02],
		ll: [81, 37.92, 0]
		},
		{code: 'EPSG:21781',
		xy: [660389.52, 185731.63, -49.23],
		ll: [8.23, 46.82, 0]
		},
		{code: 'EPSG:27563',
		xy: [653704.865208, 176887.660037],
		ll: [3.005, 43.89]
		},
		{code: 'EPSG:54029',
		xy: [1094702.50,6496789.74,-6468.21],
		ll: [11.0, 53.0,0.0]
		},
		{code: 'EPSG:54003',
		xy: [1223145.57,6491218.13,-6468.21],
		ll: [11.0, 53.0]
		},
		{code: 'EPSG:3573',
		xy: [2923052.02009, 1054885.46559],
		ll: [9.84375, 61.875]
		},
		}/*,
		{code: 'EPSG:54009',
		xy: [-10617602.79013849,4108337.84708608,0.00000000 ],
		ll: [-119,34,0]
		@@ -15,18 +63,10 @@ Proj4js.testPoints = [
		xy: [2547685.01212,5699155.7345],
		ll: [6.685,51.425]
		},
		{code: 'EPSG:54029',
		xy: [1094702.50,6496789.74,-6468.21],
		ll: [11.0, 53.0,0.0]
		},
		{code: 'EPSG:54008',
		xy: [738509.49,5874620.38],
		ll: [11.0, 53.0]
		},
		{code: 'EPSG:29100',
		xy: [5110899.06,10552971.81,-22.99],
		ll: [-53.0, 5.0,0.0] /**
		},
		{code: 'EPSG:2057',
		xy: [-11608322.26,18282612.23,-281.67],
		ll: [-53.0, 5.0,0.0]
		@@ -35,10 +75,6 @@ Proj4js.testPoints = [
		xy: [804759.21,6164983.82,-13598.03],
		ll: [11.0, 53.0, 0.0]
		},
		{code: 'EPSG:54003',
		xy: [1223145.57,6491218.13,-6468.21],
		ll: [11.0, 53.0]
		},
		{code: 'EPSG:3035',
		xy: [4388138.60, 3321736.46],
		ll: [11.0, 53.0]
		@@ -55,18 +91,14 @@ Proj4js.testPoints = [
		xy: [931625.91, 789252.65],
		ll: [-127.0, 52.11]
		},
		{code: 'EPSG:27492',
		xy: [25260.493584, -9579.245052],
		ll: [-7.84, 39.58]
		},
		{code: 'EPSG:27700',
		xy: [343733.14, 612144.53, -51.89],
		ll: [-2.89, 55.4, 0]
		},
		{code: 'EPSG:32615',
		xy: [500000, 4649776.22482],
		ll: [-93, 42]
		},
		{code: 'EPSG:26916',
		xy: [5110899.06,10552971.81,-22.99],
		ll: [-86.6056, 34.5790,0.0]
		},
		{code: 'EPSG:32612',
		xy: [383357.429537, 6684599.06392],
		ll: [-113.109375, 60.28125]
		@@ -79,10 +111,6 @@ Proj4js.testPoints = [
		xy: [450055.70, 5262356.33],
		ll: [13.33333333333, 47.5]
		},
		{code: 'EPSG:27563',
		xy: [653704.865208, 176887.660037],
		ll: [3.005, 43.89]
		},
		{code: 'EPSG:2736',
		xy: [603933.40, 7677505.64],
		ll: [34.0, -21.0]
		@@ -102,9 +130,5 @@ Proj4js.testPoints = [
		{code: 'EPSG:28992',
		xy: [148312.15, 457804.79, 698.48],
		ll: [5.29, 52.11]
		},
		{code: 'EPSG:2403',
		xy: [27500000.00, 4198690.08, -109.02],
		ll: [81, 37.92, 0] ****/
		}
		}*/
		];
		No newline at end of file