Commit 6e1ba703 authored by Richard Didier's avatar Richard Didier
Browse files

fixes for specific cases (pole, ellipsoid, sphere) : aea, aeqd : ellipsoid... 

fixes for specific cases (pole, ellipsoid, sphere) : aea, aeqd : ellipsoid case added cea : ellipsoid and sphere cases rewriting omerc : use of no_rot, no_off, lon_1 and lon_2 parameters stere : pole case added, use of lat_ts cass : pole case added. See ticket #80
parent b15f00a4

lib/projCode/aea.js

+3 −12
Original line number Diff line number Diff line
@@ -99,19 +99,10 @@ Proj4js.Proj.aea = { 
      theta = Math.atan2(con * p.x, con * p.y);

    }

    con = rh1 * this.ns0 / this.a;

    qs = (this.c - con * con) / this.ns0;

    if (this.e3 >= 1e-10) {

      con = 1 - .5 * (1.0 -this.es) * Math.log((1.0 - this.e3) / (1.0 + this.e3))/this.e3;

      if (Math.abs(Math.abs(con) - Math.abs(qs)) > .0000000001 ) {

          lat = this.phi1z(this.e3,qs);

      } else {

          if (qs >= 0) {

             lat = .5 * Proj4js.common.PI;

          } else {

             lat = -.5 * Proj4js.common.PI;

          }

      }

    if (this.sphere) {

      lat = Math.asin((this.c-con*con)/(2.0*this.ns0));

    } else {

      qs = (this.c - con * con) / this.ns0;

      lat = this.phi1z(this.e3,qs);

    }

lib/projCode/aeqd.js

+149 −51
Original line number Diff line number Diff line
@@ -8,32 +8,81 @@ Proj4js.Proj.aeqd = { 
  forward: function(p) {

    var lon=p.x;

    var lat=p.y;

    var ksp;



    var sinphi=Math.sin(p.y);

    var cosphi=Math.cos(p.y); 

    var dlon = Proj4js.common.adjust_lon(lon - this.long0);

    var coslon = Math.cos(dlon);

    var g = this.sin_p12 * sinphi + this.cos_p12 * cosphi * coslon;

    if (Math.abs(Math.abs(g) - 1.0) < Proj4js.common.EPSLN) {

       ksp = 1.0;

       if (g < 0.0) {

         Proj4js.reportError("aeqd:Fwd:PointError");

         return;

    

    if (this.sphere){

	if (Math.abs(this.sin_p12-1.0)<=Proj4js.common.EPSLN){

		//North Pole case

		p.x=this.x0 + this.a * (Proj4js.common.HALF_PI-lat) *  Math.sin(dlon);

		p.y=this.y0 - this.a * (Proj4js.common.HALF_PI-lat) *  Math.cos(dlon);

		return p;

	} else if (Math.abs(this.sin_p12+1.0)<=Proj4js.common.EPSLN){

		//South Pole case

		p.x=this.x0 + this.a * (Proj4js.common.HALF_PI+lat) *  Math.sin(dlon);

		p.y=this.y0 + this.a * (Proj4js.common.HALF_PI+lat) *  Math.cos(dlon);

		return p;

	} else {

		//default case

		var cos_c=this.sin_p12*sinphi+this.cos_p12*cosphi*Math.cos(dlon);

		var c = Math.acos(cos_c);

		var kp = c/Math.sin(c);

		p.x=this.x0 + this.a*kp*cosphi*Math.sin(dlon);

		p.y=this.y0 + this.a*kp*(this.cos_p12*sinphi-this.sin_p12*cosphi*Math.cos(dlon));

		return p;

	}

    } else {

       var z = Math.acos(g);

       ksp = z/Math.sin(z);

	var e0 = Proj4js.common.e0fn(this.es);

	var e1 = Proj4js.common.e1fn(this.es);

	var e2 = Proj4js.common.e2fn(this.es);

	var e3 = Proj4js.common.e3fn(this.es);

	if (Math.abs(this.sin_p12-1.0)<=Proj4js.common.EPSLN){

		//North Pole case

		var Mlp = this.a*Proj4js.common.mlfn(e0,e1,e2,e3,Proj4js.common.HALF_PI);

		var Ml = this.a*Proj4js.common.mlfn(e0,e1,e2,e3,lat);

		p.x = this.x0 + (Mlp-Ml)*Math.sin(dlon);

		p.y = this.y0 - (Mlp-Ml)*Math.cos(dlon);

		return p;

	} else if (Math.abs(this.sin_p12+1.0)<=Proj4js.common.EPSLN){

		//South Pole case

		var Mlp = this.a*Proj4js.common.mlfn(e0,e1,e2,e3,Proj4js.common.HALF_PI);

		var Ml = this.a*Proj4js.common.mlfn(e0,e1,e2,e3,lat);

		p.x = this.x0 + (Mlp+Ml)*Math.sin(dlon);

		p.y = this.y0 + (Mlp+Ml)*Math.cos(dlon);

		return p;

	} else {

		//Default case

		var tanphi=sinphi/cosphi;

		var Nl1 = Proj4js.common.gN(this.a,this.e, this.sin_p12);

		var Nl = Proj4js.common.gN(this.a, this.e, sinphi);

		var psi = Math.atan((1.0-this.es)*tanphi+this.es*Nl1*this.sin_p12/(Nl*cosphi));

		var Az = Math.atan2(Math.sin(dlon),this.cos_p12*Math.tan(psi)-this.sin_p12*Math.cos(dlon));

		var s;

		if (Az==0) {

			s=Math.asin(this.cos_p12*Math.sin(psi)-this.sin_p12*Math.cos(psi));

		} else if (Math.abs(Math.abs(Az)-Proj4js.common.PI)<=Proj4js.common.EPSLN){

			s=-Math.asin(this.cos_p12*Math.sin(psi)-this.sin_p12*Math.cos(psi));

		} else {

			s=Math.asin(Math.sin(dlon)*Math.cos(psi)/Math.sin(Az));

		}

    p.x = this.x0 + this.a * ksp * cosphi * Math.sin(dlon);

    p.y = this.y0 + this.a * ksp * (this.cos_p12 * sinphi - this.sin_p12 * cosphi * coslon);

		var G = this.e*this.sin_p12/Math.sqrt(1.0-this.es);

		var H = this.e*this.cos_p12*Math.cos(Az)/Math.sqrt(1.0-this.es);

		var Hs = H*H;

		var c = Nl1*s*(1.0-s*s*Hs*(1.0-Hs)/6.0+s*s*s/8.0*G*H*(1.0-2.0*Hs)+s*s*s*s/120.0*(Hs*(4.0-7.0*Hs)-3.0*G*G*(1.0-7.0*Hs))-s*s*s*s*s/48.0*G*H);

		p.x=this.x0+c*Math.sin(Az);

		p.y=this.y0+c*Math.cos(Az);

		return p;

	}

    }

    

   

  },



  inverse: function(p){

    p.x -= this.x0;

    p.y -= this.y0;



	if (this.sphere){

		var rh = Math.sqrt(p.x * p.x + p.y *p.y);

		if (rh > (2.0 * Proj4js.common.HALF_PI * this.a)) {

			Proj4js.reportError("aeqdInvDataError");
@@ -58,18 +107,67 @@ Proj4js.Proj.aeqd = { 
					lon = Proj4js.common.adjust_lon(this.long0 - Math.atan2(-p.x , p.y));

				}

			} else {

        con = cosz - this.sin_p12 * Math.sin(lat);

				/*con = cosz - this.sin_p12 * Math.sin(lat);

				if ((Math.abs(con) < Proj4js.common.EPSLN) && (Math.abs(p.x) < Proj4js.common.EPSLN)) {

					//no-op, just keep the lon value as is

				} else {

					var temp = Math.atan2((p.x * sinz * this.cos_p12), (con * rh));

					lon = Proj4js.common.adjust_lon(this.long0 + Math.atan2((p.x * sinz * this.cos_p12), (con * rh)));

				}*/

				lon=Proj4js.common.adjust_lon(this.long0+Math.atan2(p.x*sinz,rh*this.cos_p12*cosz-p.y*this.sin_p12*sinz));

			}

		}



		p.x = lon;

		p.y = lat;

		return p;

	}

	else {

		var e0 = Proj4js.common.e0fn(this.es);

		var e1 = Proj4js.common.e1fn(this.es);

		var e2 = Proj4js.common.e2fn(this.es);

		var e3 = Proj4js.common.e3fn(this.es);

		if (Math.abs(this.sin_p12-1.0)<=Proj4js.common.EPSLN){

			//North pole case

			var Mlp = this.a*Proj4js.common.mlfn(e0,e1,e2,e3,Proj4js.common.HALF_PI);

			var rh = Math.sqrt(p.x*p.x+p.y*p.y);

			var M = Mlp-rh;

			var lat = Proj4js.common.imlfn(M/this.a,e0, e1,e2,e3);

			var lon = Proj4js.common.adjust_lon(this.long0+Math.atan2(p.x,-1.0*p.y));

			p.x=lon,

			p.y=lat;

			return p;

		} else if (Math.abs(this.sin_p12+1.0)<=Proj4js.common.EPSLN){

			//South pole case

			var Mlp = this.a*Proj4js.common.mlfn(e0,e1,e2,e3,Proj4js.common.HALF_PI);

			var rh = Math.sqrt(p.x*p.x+p.y*p.y);

			var M = rh-Mlp;

			

			var lat = Proj4js.common.imlfn(M/this.a,e0, e1,e2,e3);

			var lon = Proj4js.common.adjust_lon(this.long0+Math.atan2(p.x,p.y));

			p.x=lon,

			p.y=lat;

			return p;

		} else {

			//default case

			var rh = Math.sqrt(p.x*p.x+p.y*p.y);

			var Az = Math.atan2(p.x,p.y);

			var N1 = Proj4js.common.gN(this.a, this.e, this.sin_p12);

			var cosAz = Math.cos(Az);

			var tmp = this.e*this.cos_p12*cosAz;

			var A = -tmp*tmp/(1.0 - this.es);

			var B=3.0*this.es*(1.0-A) * this.sin_p12*this.cos_p12*cosAz/(1.0-this.es);

			var D = rh/N1;

			var Ee = D-A*(1.0+A)*Math.pow(D,3.0)/6.0-B*(1+3.0*A)*Math.pow(D,4.0)/24.0;

			var F = 1.0-A*Ee*Ee/2.0-D*Ee*Ee*Ee/6.0;

			var psi = Math.asin(this.sin_p12*Math.cos(Ee)+this.cos_p12*Math.sin(Ee)*cosAz);

			var lon = Proj4js.common.adjust_lon(this.long0+Math.asin(Math.sin(Az)*Math.sin(Ee)/Math.cos(psi)));

			var lat = Math.atan((1.0-this.es*F*this.sin_p12/Math.sin(psi))*Math.tan(psi)/(1.0-this.es));

			p.x=lon;

			p.y=lat;

			return p;

		}

	}

    

  } 

};

lib/projCode/cass.js

+49 −41
Original line number Diff line number Diff line
@@ -27,16 +27,14 @@ ALGORITHM REFERENCES 
Proj4js.Proj.cass = {

  init : function() {

    if (!this.sphere) {

      this.en = Proj4js.common.pj_enfn(this.es)

      this.m0 = Proj4js.common.pj_mlfn(this.lat0, Math.sin(this.lat0), Math.cos(this.lat0), this.en);

      this.e0 = Proj4js.common.e0fn(this.es);

      this.e1 = Proj4js.common.e1fn(this.es);

      this.e2 = Proj4js.common.e2fn(this.es);

      this.e3 = Proj4js.common.e3fn(this.es);

      this.ml0 = this.a*Proj4js.common.mlfn(this.e0,this.e1,this.e2,this.e3,this.lat0);

    }

  },



  C1:	.16666666666666666666,

  C2:	.00833333333333333333,

  C3:	.04166666666666666666,

  C4:	.33333333333333333333,

  C5:	.06666666666666666666,





/* Cassini forward equations--mapping lat,long to x,y
@@ -51,25 +49,27 @@ Proj4js.Proj.cass = { 
    lam = Proj4js.common.adjust_lon(lam - this.long0);

    

    if (this.sphere) {

      x = Math.asin(Math.cos(phi) * Math.sin(lam));

      y = Math.atan2(Math.tan(phi) , Math.cos(lam)) - this.phi0;

      x = this.a*Math.asin(Math.cos(phi) * Math.sin(lam));

      y = this.a*(Math.atan2(Math.tan(phi) , Math.cos(lam)) - this.lat0);

    } else {

        //ellipsoid

      this.n = Math.sin(phi);

      this.c = Math.cos(phi);

      y = Proj4js.common.pj_mlfn(phi, this.n, this.c, this.en);

      this.n = 1./Math.sqrt(1. - this.es * this.n * this.n);

      this.tn = Math.tan(phi); 

      this.t = this.tn * this.tn;

      this.a1 = lam * this.c;

      this.c *= this.es * this.c / (1 - this.es);

      this.a2 = this.a1 * this.a1;

      x = this.n * this.a1 * (1. - this.a2 * this.t * (this.C1 - (8. - this.t + 8. * this.c) * this.a2 * this.C2));

      y -= this.m0 - this.n * this.tn * this.a2 * (.5 + (5. - this.t + 6. * this.c) * this.a2 * this.C3);

      var sinphi = Math.sin(phi);

      var cosphi = Math.cos(phi);

      var nl = Proj4js.common.gN(this.a,this.e,sinphi);

      var tl = Math.tan(phi)*Math.tan(phi);

      var al = lam*Math.cos(phi);

      var asq = al*al;

      var cl = this.es*cosphi*cosphi/(1.0-this.es);

      var ml = this.a*Proj4js.common.mlfn(this.e0,this.e1,this.e2,this.e3,phi);

      

      x = nl*al*(1.0-asq*tl*(1.0/6.0-(8.0-tl+8.0*cl)*asq/120.0));

      y = ml-this.ml0 + nl*sinphi/cosphi*asq*(0.5+(5.0-tl+6.0*cl)*asq/24.0);

      

      

    }

    

    p.x = this.a*x + this.x0;

    p.y = this.a*y + this.y0;

    p.x = x + this.x0;

    p.y = y + this.y0;

    return p;

  },//cassFwd()
@@ -83,26 +83,34 @@ Proj4js.Proj.cass = { 
    var phi, lam;



    if (this.sphere) {

      this.dd = y + this.lat0;

      phi = Math.asin(Math.sin(this.dd) * Math.cos(x));

      lam = Math.atan2(Math.tan(x), Math.cos(this.dd));

      var dd = y + this.lat0;

      phi = Math.asin(Math.sin(dd) * Math.cos(x));

      lam = Math.atan2(Math.tan(x), Math.cos(dd));

    } else {

      /* ellipsoid */

      var ph1 = Proj4js.common.pj_inv_mlfn(this.m0 + y, this.es, this.en);

      this.tn = Math.tan(ph1); 

      this.t = this.tn * this.tn;

      this.n = Math.sin(ph1);

      this.r = 1. / (1. - this.es * this.n * this.n);

      this.n = Math.sqrt(this.r);

      this.r *= (1. - this.es) * this.n;

      this.dd = x / this.n;

      this.d2 = this.dd * this.dd;

      phi = ph1 - (this.n * this.tn / this.r) * this.d2 * (.5 - (1. + 3. * this.t) * this.d2 * this.C3);

      lam = this.dd * (1. + this.t * this.d2 * (-this.C4 + (1. + 3. * this.t) * this.d2 * this.C5)) / Math.cos(ph1);

      var ml1 = this.ml0/this.a + y;

      var phi1 = Proj4js.common.imlfn(ml1, this.e0,this.e1,this.e2,this.e3);

      if (Math.abs(Math.abs(phi1)-Proj4js.common.HALF_PI)<=Proj4js.common.EPSLN){

	p.x=this.long0;

	p.y=Proj4js.common.HALF_PI;

	if (y<0.0){p.y*=-1.0;}

	return p;

      }

      var nl1 = Proj4js.common.gN(this.a,this.e, Math.sin(phi1));

      

      var rl1 = nl1*nl1*nl1/this.a/this.a*(1.0-this.es);

      var tl1 = Math.pow(Math.tan(phi1),2.0);

      var dl = x*this.a/nl1;

      var dsq=dl*dl;

      phi = phi1-nl1*Math.tan(phi1)/rl1*dl*dl*(0.5-(1.0+3.0*tl1)*dl*dl/24.0);

      lam = dl*(1.0-dsq*(tl1/3.0+(1.0+3.0*tl1)*tl1*dsq/15.0))/Math.cos(phi1);

      

    } 

    p.x = Proj4js.common.adjust_lon(this.long0+lam);

    p.y = phi;

    

    p.x=Proj4js.common.adjust_lon(lam+this.long0);

    p.y=Proj4js.common.adjust_lat(phi);

    return p;

    

   }//cassInv()



}
 
};

lib/projCode/cea.js

+20 −15
Original line number Diff line number Diff line
@@ -39,6 +39,9 @@ Proj4js.Proj.cea = { 
  -------------------------------------------*/

  init: function() {

    //no-op

    if (!this.sphere){

	    this.k0 = Proj4js.common.msfnz(this.e, Math.sin(this.lat_ts), Math.cos(this.lat_ts));

    }

  },
@@ -47,21 +50,18 @@ Proj4js.Proj.cea = { 
  forward: function(p) {

    var lon=p.x;

    var lat=p.y;

    var x,y;

    /* Forward equations

      -----------------*/

    var dlon = Proj4js.common.adjust_lon(lon -this.long0);

    var x = this.x0 + this.a * dlon * Math.cos(this.lat_ts);

    var y = this.y0 + this.a * Math.sin(lat) / Math.cos(this.lat_ts);

   /* Elliptical Forward Transform

      Not implemented due to a lack of a matchign inverse function

    {

      var Sin_Lat = Math.sin(lat);

      var Rn = this.a * (Math.sqrt(1.0e0 - this.es * Sin_Lat * Sin_Lat ));

    if (this.sphere){

	x = this.x0 + this.a * dlon * Math.cos(this.lat_ts);

      y = this.y0 + Rn * Math.sin(lat) / Math.cos(this.lat_ts);

	y = this.y0 + this.a * Math.sin(lat) / Math.cos(this.lat_ts);

    } else {

	var qs = Proj4js.common.qsfnz(this.e,Math.sin(lat));

	x = this.x0 + this.a*this.k0*dlon;

	y = this.y0 + this.a*qs*0.5/this.k0;

    }

   */





    p.x=x;

    p.y=y;
@@ -73,10 +73,15 @@ Proj4js.Proj.cea = { 
  inverse: function(p) {

    p.x -= this.x0;

    p.y -= this.y0;



    var lon = Proj4js.common.adjust_lon( this.long0 + (p.x / this.a) / Math.cos(this.lat_ts) );



    var lat = Math.asin( (p.y/this.a) * Math.cos(this.lat_ts) );

    var lon, lat;

    

    if (this.sphere){

	lon = Proj4js.common.adjust_lon( this.long0 + (p.x / this.a) / Math.cos(this.lat_ts) );

        lat = Math.asin( (p.y/this.a) * Math.cos(this.lat_ts) );

    } else {

	lat=Proj4js.common.iqsfnz(this.e,2.0*p.y*this.k0/this.a);

	lon = Proj4js.common.adjust_lon( this.long0 + p.x/(this.a*this.k0));

    }



    p.x=lon;

    p.y=lat;

lib/projCode/eqdc.js

+61 −69
Original line number Diff line number Diff line
@@ -32,8 +32,14 @@ Proj4js.Proj.eqdc = { 


    /* Place parameters in static storage for common use

      -------------------------------------------------*/



    if(!this.mode) this.mode=0;//chosen default mode

	// Standard Parallels cannot be equal and on opposite sides of the equator

	if (Math.abs(this.lat1 + this.lat2) < Proj4js.common.EPSLN) {

		Proj4js.common.reportError("eqdc:init: Equal Latitudes");

		return;

	}

	if (this.lat2 == null) {

		this.lat2=this.lat1;

	}

	this.temp = this.b / this.a;

	this.es = 1.0 - Math.pow(this.temp,2);

	this.e = Math.sqrt(this.es);
@@ -48,25 +54,15 @@ Proj4js.Proj.eqdc = { 
	this.ms1 = Proj4js.common.msfnz(this.e,this.sinphi,this.cosphi);

	this.ml1 = Proj4js.common.mlfn(this.e0, this.e1, this.e2,this.e3, this.lat1);



    /* format B

    ---------*/

    if (this.mode != 0) {

      if (Math.abs(this.lat1 + this.lat2) < Proj4js.common.EPSLN) {

	if (Math.abs(this.lat1 - this.lat2) < Proj4js.common.EPSLN) {

		this.ns = this.sinphi;

		Proj4js.reportError("eqdc:Init:EqualLatitudes");

            //return(81);

       }

        } else {

		this.sinphi=Math.sin(this.lat2);

		this.cosphi=Math.cos(this.lat2); 



		this.ms2 = Proj4js.common.msfnz(this.e,this.sinphi,this.cosphi);

		this.ml2 = Proj4js.common.mlfn(this.e0, this.e1, this.e2, this.e3, this.lat2);

       if (Math.abs(this.lat1 - this.lat2) >= Proj4js.common.EPSLN) {

		this.ns = (this.ms1 - this.ms2) / (this.ml2 - this.ml1);

       } else {

          this.ns = this.sinphi;

       }

    } else {

      this.ns = this.sinphi;

	}

	this.g = this.ml1 + this.ms1/this.ns;

	this.ml0 = Proj4js.common.mlfn(this.e0, this.e1,this. e2, this.e3, this.lat0);
@@ -79,13 +75,17 @@ Proj4js.Proj.eqdc = { 
  forward: function(p) {

    var lon=p.x;

    var lat=p.y;

    var rh1;



    /* Forward equations

      -----------------*/

    if (this.sphere){

	rh1 = this.a *(this.g - lat);

    } else {

	var ml = Proj4js.common.mlfn(this.e0, this.e1, this.e2, this.e3, lat);

    var rh1 = this.a * (this.g - ml);

	rh1 = this.a * (this.g - ml);

    }

    var theta = this.ns * Proj4js.common.adjust_lon(lon - this.long0);



    var x = this.x0  + rh1 * Math.sin(theta);

    var y = this.y0 + this.rh - rh1 * Math.cos(theta);

    p.x=x;
@@ -98,7 +98,7 @@ Proj4js.Proj.eqdc = { 
  inverse: function(p) {

    p.x -= this.x0;

    p.y  = this.rh - p.y + this.y0;

    var con, rh1;

    var con, rh1, lat, lon;

    if (this.ns >= 0) {

       rh1 = Math.sqrt(p.x *p.x + p.y * p.y); 

       con = 1.0;
@@ -107,34 +107,26 @@ Proj4js.Proj.eqdc = { 
       con = -1.0;

    }

    var theta = 0.0;

    if (rh1 != 0.0) theta = Math.atan2(con *p.x, con *p.y);

    var ml = this.g - rh1 /this.a;

    var lat = this.phi3z(ml,this.e0,this.e1,this.e2,this.e3);

    var lon = Proj4js.common.adjust_lon(this.long0 + theta / this.ns);

    if (rh1 != 0.0) {theta = Math.atan2(con *p.x, con *p.y);}

    

    if (this.sphere){

	lon=Proj4js.common.adjust_lon(this.long0+theta/this.ns);

	lat=Proj4js.common.adjust_lat(this.g-rh1/this.a);

	p.x=lon;

	p.y=lat;

	return p;

    } else {

	var ml = this.g - rh1 /this.a;

	lat = Proj4js.common.imlfn(ml,this.e0,this.e1,this.e2,this.e3);

	lon = Proj4js.common.adjust_lon(this.long0 + theta / this.ns);

	p.x=lon;

	p.y=lat;  

	return p;

    },

    

/* Function to compute latitude, phi3, for the inverse of the Equidistant

   Conic projection.

-----------------------------------------------------------------*/

  phi3z: function(ml,e0,e1,e2,e3) {

    var phi;

    var dphi;



    phi = ml;

    for (var i = 0; i < 15; i++) {

      dphi = (ml + e1 * Math.sin(2.0 * phi) - e2 * Math.sin(4.0 * phi) + e3 * Math.sin(6.0 * phi))/ e0 - phi;

      phi += dphi;

      if (Math.abs(dphi) <= .0000000001) {

        return phi;

      }

    }

    Proj4js.reportError("PHI3Z-CONV:Latitude failed to converge after 15 iterations");

    return null;

    

    }

    





    

};
 
 No newline at end of file