/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/stere.php |
---|
New file |
0,0 → 1,303 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
// Initialize the Stereographic projection |
class Proj4phpProjStere { |
protected $TOL = 1.e-8; |
protected $NITER = 8; |
protected $CONV = 1.e-10; |
protected $S_POLE = 0; |
protected $N_POLE = 1; |
protected $OBLIQ = 2; |
protected $EQUIT = 3; |
/** |
* |
* @param type $phit |
* @param type $sinphi |
* @param type $eccen |
* @return type |
*/ |
public function ssfn_( $phit, $sinphi, $eccen ) { |
$sinphi *= $eccen; |
return (tan( .5 * (Proj4php::$common->HALF_PI + $phit) ) * pow( (1. - $sinphi) / (1. + $sinphi), .5 * $eccen )); |
} |
/** |
* |
*/ |
public function init() { |
$this->phits = $this->lat_ts ? $this->lat_ts : Proj4php::$common->HALF_PI; |
$t = abs( $this->lat0 ); |
if( (abs( $t ) - Proj4php::$common->HALF_PI) < Proj4php::$common->EPSLN ) { |
$this->mode = $this->lat0 < 0. ? $this->S_POLE : $this->N_POLE; |
} else { |
$this->mode = $t > Proj4php::$common->EPSLN ? $this->OBLIQ : $this->EQUIT; |
} |
$this->phits = abs( $this->phits ); |
if( $this->es ) { |
#$X; |
switch( $this->mode ) { |
case $this->N_POLE: |
case $this->S_POLE: |
if( abs( $this->phits - Proj4php::$common->HALF_PI ) < Proj4php::$common->EPSLN ) { |
$this->akm1 = 2. * $this->k0 / sqrt( pow( 1 + $this->e, 1 + $this->e ) * pow( 1 - $this->e, 1 - $this->e ) ); |
} else { |
$t = sin( $this->phits ); |
$this->akm1 = cos( $this->phits ) / Proj4php::$common->tsfnz( $this->e, $this->phits, $t ); |
$t *= $this->e; |
$this->akm1 /= sqrt( 1. - $t * $t ); |
} |
break; |
case $this->EQUIT: |
$this->akm1 = 2. * $this->k0; |
break; |
case $this->OBLIQ: |
$t = sin( $this->lat0 ); |
$X = 2. * atan( $this->ssfn_( $this->lat0, $t, $this->e ) ) - Proj4php::$common->HALF_PI; |
$t *= $this->e; |
$this->akm1 = 2. * $this->k0 * cos( $this->lat0 ) / sqrt( 1. - $t * $t ); |
$this->sinX1 = sin( $X ); |
$this->cosX1 = cos( $X ); |
break; |
} |
} else { |
switch( $this->mode ) { |
case $this->OBLIQ: |
$this->sinph0 = sin( $this->lat0 ); |
$this->cosph0 = cos( $this->lat0 ); |
case $this->EQUIT: |
$this->akm1 = 2. * $this->k0; |
break; |
case $this->S_POLE: |
case $this->N_POLE: |
$this->akm1 = abs( $this->phits - Proj4php::$common->HALF_PI ) >= Proj4php::$common->EPSLN ? |
cos( $this->phits ) / tan( Proj4php::$common->FORTPI - .5 * $this->phits ) : |
2. * $this->k0; |
break; |
} |
} |
} |
/** |
* Stereographic forward equations--mapping lat,long to x,y |
* |
* @param type $p |
* @return type |
*/ |
public function forward( $p ) { |
$lon = $p->x; |
$lon = Proj4php::$common->adjust_lon( $lon - $this->long0 ); |
$lat = $p->y; |
#$x; |
#$y; |
if( $this->sphere ) { |
/* |
$sinphi; |
$cosphi; |
$coslam; |
$sinlam; |
*/ |
$sinphi = sin( $lat ); |
$cosphi = cos( $lat ); |
$coslam = cos( $lon ); |
$sinlam = sin( $lon ); |
switch( $this->mode ) { |
case $this->EQUIT: |
$y = 1. + $cosphi * $coslam; |
if( y <= Proj4php::$common->EPSLN ) { |
Proj4php::reportError("stere:forward:Equit"); |
} |
$y = $this->akm1 / $y; |
$x = $y * $cosphi * $sinlam; |
$y *= $sinphi; |
break; |
case $this->OBLIQ: |
$y = 1. + $this->sinph0 * $sinphi + $this->cosph0 * $cosphi * $coslam; |
if( $y <= Proj4php::$common->EPSLN ) { |
Proj4php::reportError("stere:forward:Obliq"); |
} |
$y = $this->akm1 / $y; |
$x = $y * $cosphi * $sinlam; |
$y *= $this->cosph0 * $sinphi - $this->sinph0 * $cosphi * $coslam; |
break; |
case $this->N_POLE: |
$coslam = -$coslam; |
$lat = -$lat; |
//Note no break here so it conitnues through S_POLE |
case $this->S_POLE: |
if( abs( $lat - Proj4php::$common->HALF_PI ) < $this->TOL ) { |
Proj4php::reportError("stere:forward:S_POLE"); |
} |
$y = $this->akm1 * tan( Proj4php::$common->FORTPI + .5 * $lat ); |
$x = $sinlam * $y; |
$y *= $coslam; |
break; |
} |
} else { |
$coslam = cos( $lon ); |
$sinlam = sin( $lon ); |
$sinphi = sin( $lat ); |
if( $this->mode == $this->OBLIQ || $this->mode == $this->EQUIT ) { |
$Xt = 2. * atan( $this->ssfn_( $lat, $sinphi, $this->e ) ); |
$sinX = sin( $Xt - Proj4php::$common->HALF_PI ); |
$cosX = cos( $Xt ); |
} |
switch( $this->mode ) { |
case $this->OBLIQ: |
$A = $this->akm1 / ($this->cosX1 * (1. + $this->sinX1 * $sinX + $this->cosX1 * $cosX * $coslam)); |
$y = $A * ($this->cosX1 * $sinX - $this->sinX1 * $cosX * $coslam); |
$x = $A * $cosX; |
break; |
case $this->EQUIT: |
$A = 2. * $this->akm1 / (1. + $cosX * $coslam); |
$y = $A * $sinX; |
$x = $A * $cosX; |
break; |
case $this->S_POLE: |
$lat = -$lat; |
$coslam = - $coslam; |
$sinphi = -$sinphi; |
case $this->N_POLE: |
$x = $this->akm1 * Proj4php::$common->tsfnz( $this->e, $lat, $sinphi ); |
$y = - $x * $coslam; |
break; |
} |
$x = $x * $sinlam; |
} |
$p->x = $x * $this->a + $this->x0; |
$p->y = $y * $this->a + $this->y0; |
return $p; |
} |
/** |
* Stereographic inverse equations--mapping x,y to lat/long |
* |
* @param type $p |
* @return type |
*/ |
public function inverse( $p ) { |
$x = ($p->x - $this->x0) / $this->a; /* descale and de-offset */ |
$y = ($p->y - $this->y0) / $this->a; |
/* |
$lon; |
$lat; |
$cosphi; |
$sinphi; |
$rho; |
$tp = 0.0; |
$phi_l = 0.0; |
$i; |
*/ |
$halfe = 0.0; |
$pi2 = 0.0; |
if( $this->sphere ) { |
/* |
$c; |
$rh; |
$sinc; |
$cosc; |
*/ |
$rh = sqrt( $x * $x + $y * $y ); |
$c = 2. * atan( $rh / $this->akm1 ); |
$sinc = sin( $c ); |
$cosc = cos( $c ); |
$lon = 0.; |
switch( $this->mode ) { |
case $this->EQUIT: |
if( abs( $rh ) <= Proj4php::$common->EPSLN ) { |
$lat = 0.; |
} else { |
$lat = asin( $y * $sinc / $rh ); |
} |
if( $cosc != 0. || $x != 0. ) |
$lon = atan2( $x * $sinc, $cosc * $rh ); |
break; |
case $this->OBLIQ: |
if( abs( $rh ) <= Proj4php::$common->EPSLN ) { |
$lat = $this->phi0; |
} else { |
$lat = asin( $cosc * $this->sinph0 + $y * $sinc * $this->cosph0 / $rh ); |
} |
$c = $cosc - $this->sinph0 * sin( $lat ); |
if( $c != 0. || $x != 0. ) { |
$lon = atan2( $x * $sinc * $this->cosph0, $c * $rh ); |
} |
break; |
case $this->N_POLE: |
$y = -$y; |
case $this->S_POLE: |
if( abs( $rh ) <= Proj4php::$common->EPSLN ) { |
$lat = $this->phi0; |
} else { |
$lat = asin( $this->mode == $this->S_POLE ? -$cosc : $cosc ); |
} |
$lon = ($x == 0. && $y == 0.) ? 0. : atan2( $x, $y ); |
break; |
} |
$p->x = Proj4php::$common->adjust_lon( $lon + $this->long0 ); |
$p->y = $lat; |
} else { |
$rho = sqrt( $x * $x + $y * $y ); |
switch( $this->mode ) { |
case $this->OBLIQ: |
case $this->EQUIT: |
$tp = 2. * atan2( $rho * $this->cosX1, $this->akm1 ); |
$cosphi = cos( $tp ); |
$sinphi = sin( $tp ); |
if( $rho == 0.0 ) { |
$phi_l = asin( $cosphi * $this->sinX1 ); |
} else { |
$phi_l = asin( $cosphi * $this->sinX1 + ($y * $sinphi * $this->cosX1 / $rho) ); |
} |
$tp = tan( .5 * (Proj4php::$common->HALF_PI + $phi_l) ); |
$x *= $sinphi; |
$y = $rho * $this->cosX1 * $cosphi - $y * $this->sinX1 * $sinphi; |
$pi2 = Proj4php::$common->HALF_PI; |
$halfe = .5 * $this->e; |
break; |
case $this->N_POLE: |
$y = -$y; |
case $this->S_POLE: |
$tp = - $rho / $this->akm1; |
$phi_l = Proj4php::$common->HALF_PI - 2. * atan( $tp ); |
$pi2 = -Proj4php::$common->HALF_PI; |
$halfe = -.5 * $this->e; |
break; |
} |
for( $i = $this->NITER; $i--; $phi_l = $lat ) { //check this |
$sinphi = $this->e * sin( $phi_l ); |
$lat = 2. * atan( $tp * pow( (1. + $sinphi) / (1. - $sinphi), $halfe ) ) - $pi2; |
if( abs( phi_l - lat ) < $this->CONV ) { |
if( $this->mode == $this->S_POLE ) |
$lat = -$lat; |
$lon = ($x == 0. && $y == 0.) ? 0. : atan2( $x, $y ); |
$p->x = Proj4php::$common->adjust_lon( $lon + $this->long0 ); |
$p->y = $lat; |
return $p; |
} |
} |
} |
} |
} |
Proj4php::$proj['stere'] = new Proj4phpProjStere(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/gauss.php |
---|
New file |
0,0 → 1,74 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
class Proj4phpProjGauss { |
/** |
* |
*/ |
public function init() { |
$sphi = sin( $this->lat0 ); |
$cphi = cos( $this->lat0 ); |
$cphi *= $cphi; |
$this->rc = sqrt( 1.0 - $this->es ) / (1.0 - $this->es * $sphi * $sphi); |
$this->C = sqrt( 1.0 + $this->es * $cphi * $cphi / (1.0 - $this->es) ); |
$this->phic0 = asin( $sphi / $this->C ); |
$this->ratexp = 0.5 * $this->C * $this->e; |
$this->K = tan( 0.5 * $this->phic0 + Proj4php::$common->FORTPI ) / (pow( tan( 0.5 * $this->lat0 + Proj4php::$common->FORTPI ), $this->C ) * Proj4php::$common->srat( $this->e * $sphi, $this->ratexp )); |
} |
/** |
* |
* @param type $p |
* @return type |
*/ |
public function forward( $p ) { |
$lon = $p->x; |
$lat = $p->y; |
$p->y = 2.0 * atan( $this->K * pow( tan( 0.5 * $lat + Proj4php::$common->FORTPI ), $this->C ) * Proj4php::$common->srat( $this->e * sin( $lat ), $this->ratexp ) ) - Proj4php::$common->HALF_PI; |
$p->x = $this->C * $lon; |
return $p; |
} |
/** |
* |
* @param type $p |
* @return null |
*/ |
public function inverse( $p ) { |
$DEL_TOL = 1e-14; |
$lon = $p->x / $this->C; |
$lat = $p->y; |
$num = pow( tan( 0.5 * $lat + Proj4php::$common . FORTPI ) / $this->K, 1. / $this->C ); |
for( $i = Proj4php::$common . MAX_ITER; $i > 0; --$i ) { |
$lat = 2.0 * atan( $num * Proj4php::$common->srat( $this->e * sin( $p->y ), -0.5 * $this->e ) ) - Proj4php::$common->HALF_PI; |
if( abs( $lat - $p->y ) < $DEL_TOL ) |
break; |
$p->y = $lat; |
} |
/* convergence failed */ |
if( !$i ) { |
Proj4php::reportError( "gauss:inverse:convergence failed" ); |
return null; |
} |
$p->x = $lon; |
$p->y = $lat; |
return $p; |
} |
} |
Proj4php::$proj['gauss'] = new Proj4phpProjGauss(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/sterea.php |
---|
New file |
0,0 → 1,91 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
class Proj4phpProjSterea { |
protected $dependsOn = 'gauss'; |
/** |
* |
* @return void |
*/ |
public function init() { |
if( !$this->rc ) { |
Proj4php::reportError( "sterea:init:E_ERROR_0" ); |
return; |
} |
$this->sinc0 = sin( $this->phic0 ); |
$this->cosc0 = cos( $this->phic0 ); |
$this->R2 = 2.0 * $this->rc; |
if( !$this->title ) |
$this->title = "Oblique Stereographic Alternative"; |
} |
/** |
* |
* @param type $p |
* @return type |
*/ |
public function forward( $p ) { |
$p->x = Proj4php::$common->adjust_lon( $p->x - $this->long0 ); /* adjust del longitude */ |
$p = Proj4php::$proj['gauss']->forward( $p ); |
$sinc = sin( $p->y ); |
$cosc = cos( $p->y ); |
$cosl = cos( $p->x ); |
$k = $this->k0 * $this->R2 / (1.0 + $this->sinc0 * $sinc + $this->cosc0 * $cosc * $cosl); |
$p->x = $k * $cosc * sin( $p->x ); |
$p->y = $k * ($this->cosc0 * sinc - $this->sinc0 * $cosc * $cosl); |
$p->x = $this->a * $p->x + $this->x0; |
$p->y = $this->a * $p->y + $this->y0; |
return $p; |
} |
/** |
* |
* @param type $p |
* @return type |
*/ |
public function inverse( $p ) { |
#$lon; |
#$lat; |
$p->x = ($p->x - $this->x0) / $this->a; /* descale and de-offset */ |
$p->y = ($p->y - $this->y0) / $this->a; |
$p->x /= $this->k0; |
$p->y /= $this->k0; |
if( ($rho = sqrt( $p->x * $p->x + $p->y * $p->y ) ) ) { |
$c = 2.0 * atan2( $rho, $this->R2 ); |
$sinc = sin( $c ); |
$cosc = cos( $c ); |
$lat = asin( $cosc * $this->sinc0 + $p->y * $sinc * $this->cosc0 / $rho ); |
$lon = atan2( $p->x * $sinc, $rho * $this->cosc0 * $cosc - $p->y * $this->sinc0 * $sinc ); |
} else { |
$lat = $this->phic0; |
$lon = 0.; |
} |
$p->x = $lon; |
$p->y = $lat; |
$p = Proj4php::$proj['gauss']->inverse( $p ); |
$p->x = Proj4php::$common->adjust_lon( $p->x + $this->long0 ); /* adjust longitude to CM */ |
return $p; |
} |
} |
Proj4php::$proj['sterea'] = new Proj4phpProjSterea(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/poly.php |
---|
New file |
0,0 → 1,192 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/* Function to compute, phi4, the latitude for the inverse of the |
Polyconic projection. |
------------------------------------------------------------ */ |
function phi4z( $eccent, $e0, $e1, $e2, $e3, $a, $b, &$c, $phi ) { |
/* |
$sinphi; |
$sin2ph; |
$tanph; |
$ml; |
$mlp; |
$con1; |
$con2; |
$con3; |
$dphi; |
$i; |
*/ |
$phi = $a; |
for( $i = 1; $i <= 15; $i++ ) { |
$sinphi = sin( $phi ); |
$tanphi = tan( $phi ); |
$c = $tanphi * sqrt( 1.0 - $eccent * $sinphi * $sinphi ); |
$sin2ph = sin( 2.0 * $phi ); |
/* |
ml = e0 * *phi - e1 * sin2ph + e2 * sin (4.0 * *phi); |
mlp = e0 - 2.0 * e1 * cos (2.0 * *phi) + 4.0 * e2 * cos (4.0 * *phi); |
*/ |
$ml = $e0 * $phi - $e1 * $sin2ph + $e2 * sin( 4.0 * $phi ) - $e3 * sin( 6.0 * $phi ); |
$mlp = $e0 - 2.0 * $e1 * cos( 2.0 * $phi ) + 4.0 * $e2 * cos( 4.0 * $phi ) - 6.0 * $e3 * cos( 6.0 * $phi ); |
$con1 = 2.0 * $ml + $c * ($ml * $ml + $b) - 2.0 * $a * ($c * $ml + 1.0); |
$con2 = $eccent * $sin2ph * ($ml * $ml + $b - 2.0 * $a * $ml) / (2.0 * $c); |
$con3 = 2.0 * ($a - $ml) * ($c * $mlp - 2.0 / $sin2ph) - 2.0 * $mlp; |
$dphi = $con1 / ($con2 + $con3); |
$phi += $dphi; |
if( abs( $dphi ) <= .0000000001 ) |
return($phi); |
} |
Proj4php::reportError( "phi4z: No convergence" ); |
return null; |
} |
/* Function to compute the constant e4 from the input of the eccentricity |
of the spheroid, x. This constant is used in the Polar Stereographic |
projection. |
-------------------------------------------------------------------- */ |
function e4fn( $x ) { |
#$con; |
#$com; |
$con = 1.0 + $x; |
$com = 1.0 - $x; |
return (sqrt( (pow( $con, $con )) * (pow( $com, $com )) )); |
} |
/* * ***************************************************************************** |
NAME POLYCONIC |
PURPOSE: Transforms input longitude and latitude to Easting and |
Northing for the Polyconic projection. The |
longitude and latitude must be in radians. The Easting |
and Northing values will be returned in meters. |
PROGRAMMER DATE |
---------- ---- |
T. Mittan Mar, 1993 |
ALGORITHM REFERENCES |
1. 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. |
2. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", |
U.S. Geological Survey Professional Paper 1453 , United State Government |
Printing Office, Washington D.C., 1989. |
* ***************************************************************************** */ |
class Proj4phpProjPoly { |
/* Initialize the POLYCONIC projection |
---------------------------------- */ |
public function init() { |
#$temp; /* temporary variable */ |
if( $this->lat0 == 0 ) |
$this->lat0 = 90; //$this->lat0 ca |
/* Place parameters in static storage for common use |
------------------------------------------------- */ |
$this->temp = $this->b / $this->a; |
$this->es = 1.0 - pow( $this->temp, 2 ); // devait etre dans tmerc.js mais n y est pas donc je commente sinon retour de valeurs nulles |
$this->e = sqrt( $this->es ); |
$this->e0 = Proj4php::$common->e0fn( $this->es ); |
$this->e1 = Proj4php::$common->e1fn( $this->es ); |
$this->e2 = Proj4php::$common->e2fn( $this->es ); |
$this->e3 = Proj4php::$common->e3fn( $this->es ); |
$this->ml0 = Proj4php::$common->mlfn( $this->e0, $this->e1, $this->e2, $this->e3, $this->lat0 ); //si que des zeros le calcul ne se fait pas |
//if (!$this->ml0) {$this->ml0=0;} |
} |
/* Polyconic forward equations--mapping lat,long to x,y |
--------------------------------------------------- */ |
public function forward( $p ) { |
/* |
$sinphi; |
$cosphi; // sin and cos value |
$al; // temporary values |
$c; // temporary values |
$con; |
$ml; // cone constant, small m |
$ms; // small m |
$x; |
$y; |
*/ |
$lon = $p->x; |
$lat = $p->y; |
$con = Proj4php::$common->adjust_lon( $lon - $this->long0 ); |
if( abs( $lat ) <= .0000001 ) { |
$x = $this->x0 + $this->a * $con; |
$y = $this->y0 - $this->a * $this->ml0; |
} else { |
$sinphi = sin( $lat ); |
$cosphi = cos( $lat ); |
$ml = Proj4php::$common->mlfn( $this->e0, $this->e1, $this->e2, $this->e3, $lat ); |
$ms = Proj4php::$common->msfnz( $this->e, $sinphi, $cosphi ); |
$x = $this->x0 + $this->a * $ms * sin( $sinphi ) / $sinphi; |
$y = $this->y0 + $this->a * ($ml - $this->ml0 + $ms * (1.0 - cos( $sinphi )) / $sinphi); |
} |
$p->x = $x; |
$p->y = $y; |
return $p; |
} |
/* Inverse equations |
----------------- */ |
public function inverse( $p ) { |
/* |
$sin_phi; |
$cos_phi; // sin and cos values |
$al; // temporary values |
$b; // temporary values |
$c; // temporary values |
$con; |
$ml; // cone constant, small m |
$iflg; // error flag |
$lon; |
$lat; |
*/ |
$p->x -= $this->x0; |
$p->y -= $this->y0; |
$al = $this->ml0 + $p->y / $this->a; |
$iflg = 0; |
if( abs( $al ) <= .0000001 ) { |
$lon = $p->x / $this->a + $this->long0; |
$lat = 0.0; |
} else { |
$b = $al * $al + ($p->x / $this->a) * ($p->x / $this->a); |
$iflg = phi4z( $this->es, $this->e0, $this->e1, $this->e2, $this->e3, $this->al, $b, $c, $lat ); |
if( $iflg != 1 ) |
return($iflg); |
$lon = Proj4php::$common->adjust_lon( (Proj4php::$common->asinz( $p->x * $c / $this->a ) / sin( $lat )) + $this->long0 ); |
} |
$p->x = $lon; |
$p->y = $lat; |
return $p; |
} |
} |
Proj4php::$proj['poly'] = new Proj4phpProjPoly(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/aea.php |
---|
New file |
0,0 → 1,184 |
<?php |
/******************************************************************************* |
NAME ALBERS CONICAL EQUAL AREA |
PURPOSE: Transforms input longitude and latitude to Easting and Northing |
for the Albers Conical Equal Area projection. The longitude |
and latitude must be in radians. The Easting and Northing |
values will be returned in meters. |
PROGRAMMER DATE |
---------- ---- |
T. Mittan, Feb, 1992 |
ALGORITHM REFERENCES |
1. 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. |
2. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", |
U.S. Geological Survey Professional Paper 1453 , United State Government |
Printing Office, Washington D.C., 1989. |
*******************************************************************************/ |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
class Proj4phpProjAea { |
/** |
* |
* @return void |
*/ |
public function init() { |
if( abs( $this->lat1 + $this->lat2 ) < Proj4php::$common->EPSLN ) { |
Proj4php::reportError( "aeaInitEqualLatitudes" ); |
return; |
} |
$this->temp = $this->b / $this->a; |
$this->es = 1.0 - pow( $this->temp, 2 ); |
$this->e3 = sqrt( $this->es ); |
$this->sin_po = sin( $this->lat1 ); |
$this->cos_po = cos( $this->lat1 ); |
$this->t1 = $this->sin_po; |
$this->con = $this->sin_po; |
$this->ms1 = Proj4php::$common->msfnz( $this->e3, $this->sin_po, $this->cos_po ); |
$this->qs1 = Proj4php::$common->qsfnz( $this->e3, $this->sin_po, $this->cos_po ); |
$this->sin_po = sin( $this->lat2 ); |
$this->cos_po = cos( $this->lat2 ); |
$this->t2 = $this->sin_po; |
$this->ms2 = Proj4php::$common->msfnz( $this->e3, $this->sin_po, $this->cos_po ); |
$this->qs2 = Proj4php::$common->qsfnz( $this->e3, $this->sin_po, $this->cos_po ); |
$this->sin_po = sin( $this->lat0 ); |
$this->cos_po = cos( $this->lat0 ); |
$this->t3 = $this->sin_po; |
$this->qs0 = Proj4php::$common->qsfnz( $this->e3, $this->sin_po, $this->cos_po ); |
if( abs( $this->lat1 - $this->lat2 ) > Proj4php::$common->EPSLN ) { |
$this->ns0 = ($this->ms1 * $this->ms1 - $this->ms2 * $this->ms2) / ($this->qs2 - $this->qs1); |
} else { |
$this->ns0 = $this->con; |
} |
$this->c = $this->ms1 * $this->ms1 + $this->ns0 * $this->qs1; |
$this->rh = $this->a * sqrt( $this->c - $this->ns0 * $this->qs0 ) / $this->ns0; |
} |
/** |
* Albers Conical Equal Area forward equations--mapping lat,long to x,y |
* |
* @param Point $p |
* @return Point $p |
*/ |
public function forward( $p ) { |
$lon = $p->x; |
$lat = $p->y; |
$this->sin_phi = sin( $lat ); |
$this->cos_phi = cos( $lat ); |
$qs = Proj4php::$common->qsfnz( $this->e3, $this->sin_phi, $this->cos_phi ); |
$rh1 = $this->a * sqrt( $this->c - $this->ns0 * $qs ) / $this->ns0; |
$theta = $this->ns0 * Proj4php::$common->adjust_lon( $lon - $this->long0 ); |
$x = rh1 * sin( $theta ) + $this->x0; |
$y = $this->rh - $rh1 * cos( $theta ) + $this->y0; |
$p->x = $x; |
$p->y = $y; |
return $p; |
} |
/** |
* |
* @param Point $p |
* @return Point $p |
*/ |
public function inverse( $p ) { |
$p->x -= $this->x0; |
$p->y = $this->rh - $p->y + $this->y0; |
if( $this->ns0 >= 0 ) { |
$rh1 = sqrt( $p->x * $p->x + $p->y * $p->y ); |
$con = 1.0; |
} else { |
$rh1 = -sqrt( $p->x * $p->x + $p->y * $p->y ); |
$con = -1.0; |
} |
$theta = 0.0; |
if( $rh1 != 0.0 ) { |
$theta = 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) * log( (1.0 - $this->e3) / (1.0 + $this->e3) ) / $this->e3; |
if( abs( abs( $con ) - abs( $qs ) ) > .0000000001 ) { |
$lat = $this->phi1z( $this->e3, $qs ); |
} else { |
if( $qs >= 0 ) { |
$lat = .5 * Proj4php::$Common->PI; |
} else { |
$lat = -.5 * Proj4php::$Common->PI; |
} |
} |
} else { |
$lat = $this->phi1z( $this->e3, $qs ); |
} |
$lon = Proj4php::$common->adjust_lon( $theta / $this->ns0 + $this->long0 ); |
$p->x = $lon; |
$p->y = $lat; |
return $p; |
} |
/** |
* Function to compute phi1, the latitude for the inverse of the Albers Conical Equal-Area projection. |
* |
* @param type $eccent |
* @param type $qs |
* @return $phi or null on Convergence error |
*/ |
public function phi1z( $eccent, $qs ) { |
$phi = Proj4php::$common->asinz( .5 * $qs ); |
if( $eccent < Proj4php::$common->EPSLN ) |
return $phi; |
$eccnts = $eccent * $eccent; |
for( $i = 1; $i <= 25; ++$i ) { |
$sinphi = sin( $phi ); |
$cosphi = cos( $phi ); |
$con = $eccent * $sinphi; |
$com = 1.0 - $con * $con; |
$dphi = .5 * $com * $com / $cosphi * ($qs / (1.0 - $eccnts) - $sinphi / $com + .5 / $eccent * log( (1.0 - $con) / (1.0 + $con) )); |
$phi = $phi + $dphi; |
if( abs( $dphi ) <= 1e-7 ) |
return $phi; |
} |
Proj4php::reportError( "aea:phi1z:Convergence error" ); |
return null; |
} |
} |
Proj4php::$proj['aea'] = new Proj4phpProjAea(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/merc.php |
---|
New file |
0,0 → 1,129 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/* * ***************************************************************************** |
NAME MERCATOR |
PURPOSE: Transforms input longitude and latitude to Easting and |
Northing for the Mercator projection. The |
longitude and latitude must be in radians. The Easting |
and Northing values will be returned in meters. |
PROGRAMMER DATE |
---------- ---- |
D. Steinwand, EROS Nov, 1991 |
T. Mittan Mar, 1993 |
ALGORITHM REFERENCES |
1. 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. |
2. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", |
U.S. Geological Survey Professional Paper 1453 , United State Government |
Printing Office, Washington D.C., 1989. |
* ***************************************************************************** */ |
//static double r_major = a; /* major axis */ |
//static double r_minor = b; /* minor axis */ |
//static double lon_center = long0; /* Center longitude (projection center) */ |
//static double lat_origin = lat0; /* center latitude */ |
//static double e,es; /* eccentricity constants */ |
//static double m1; /* small value m */ |
//static double false_northing = y0; /* y offset in meters */ |
//static double false_easting = x0; /* x offset in meters */ |
//scale_fact = k0 |
class Proj4phpProjMerc { |
public function init() { |
//?$this->temp = $this->r_minor / $this->r_major; |
//$this->temp = $this->b / $this->a; |
//$this->es = 1.0 - sqrt($this->temp); |
//$this->e = sqrt( $this->es ); |
//?$this->m1 = cos($this->lat_origin) / (sqrt( 1.0 - $this->es * sin($this->lat_origin) * sin($this->lat_origin))); |
//$this->m1 = cos(0.0) / (sqrt( 1.0 - $this->es * sin(0.0) * sin(0.0))); |
if( $this->lat_ts ) { |
if( $this->sphere ) { |
$this->k0 = cos( $this->lat_ts ); |
} else { |
$this->k0 = Proj4php::$common->msfnz( $this->es, sin( $this->lat_ts ), cos( $this->lat_ts ) ); |
} |
} |
} |
/* Mercator forward equations--mapping lat,long to x,y |
-------------------------------------------------- */ |
public function forward( $p ) { |
//alert("ll2m coords : ".coords); |
$lon = $p->x; |
$lat = $p->y; |
// convert to radians |
if( $lat * Proj4php::$common->R2D > 90.0 && |
$lat * Proj4php::$common->R2D < -90.0 && |
$lon * Proj4php::$common->R2D > 180.0 && |
$lon * Proj4php::$common->R2D < -180.0 ) { |
Proj4php::reportError( "merc:forward: llInputOutOfRange: " . $lon . " : " . $lat ); |
return null; |
} |
if( abs( abs( $lat ) - Proj4php::$common->HALF_PI ) <= Proj4php::$common->EPSLN ) { |
Proj4php::reportError( "merc:forward: ll2mAtPoles" ); |
return null; |
} else { |
if( $this->sphere ) { |
$x = $this->x0 + $this->a * $this->k0 * Proj4php::$common->adjust_lon( $lon - $this->long0 ); |
$y = $this->y0 + $this->a * $this->k0 * log( tan( Proj4php::$common->FORTPI + 0.5 * $lat ) ); |
} else { |
$sinphi = sin( lat ); |
$ts = Proj4php::$common . tsfnz( $this->e, $lat, $sinphi ); |
$x = $this->x0 + $this->a * $this->k0 * Proj4php::$common->adjust_lon( $lon - $this->long0 ); |
$y = $this->y0 - $this->a * $this->k0 * log( $ts ); |
} |
$p->x = $x; |
$p->y = $y; |
return $p; |
} |
} |
/* Mercator inverse equations--mapping x,y to lat/long |
-------------------------------------------------- */ |
public function inverse( $p ) { |
$x = $p->x - $this->x0; |
$y = $p->y - $this->y0; |
if( $this->sphere ) { |
$lat = Proj4php::$common->HALF_PI - 2.0 * atan( exp( -$y / $this->a * $this->k0 ) ); |
} else { |
$ts = exp( -$y / ($this->a * $this->k0) ); |
$lat = Proj4php::$common->phi2z( $this->e, $ts ); |
if( $lat == -9999 ) { |
Proj4php::reportError( "merc:inverse: lat = -9999" ); |
return null; |
} |
} |
$lon = Proj4php::$common->adjust_lon( $this->long0 + $x / ($this->a * $this->k0) ); |
$p->x = $lon; |
$p->y = $lat; |
return $p; |
} |
} |
Proj4php::$proj['merc'] = new Proj4phpProjMerc(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/somerc.php |
---|
New file |
0,0 → 1,135 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/******************************************************************************* |
NAME SWISS OBLIQUE MERCATOR |
PURPOSE: Swiss projection. |
WARNING: X and Y are inverted (weird) in the swiss coordinate system. Not |
here, since we want X to be horizontal and Y vertical. |
ALGORITHM REFERENCES |
1. "Formules et constantes pour le Calcul pour la |
projection cylindrique conforme à axe oblique et pour la transformation entre |
des systèmes de référence". |
http://www.swisstopo.admin.ch/internet/swisstopo/fr/home/topics/survey/sys/refsys/switzerland.parsysrelated1.31216.downloadList.77004.DownloadFile.tmp/swissprojectionfr.pdf |
*******************************************************************************/ |
class Proj4phpProjSomerc { |
/** |
* |
*/ |
public function init() { |
$phy0 = $this->lat0; |
$this->lambda0 = $this->long0; |
$sinPhy0 = sin( $phy0 ); |
$semiMajorAxis = $this->a; |
$invF = $this->rf; |
$flattening = 1 / $invF; |
$e2 = 2 * $flattening - pow( $flattening, 2 ); |
$e = $this->e = sqrt( $e2 ); |
$this->R = $this->k0 * $semiMajorAxis * sqrt( 1 - $e2 ) / (1 - $e2 * pow( $sinPhy0, 2.0 )); |
$this->alpha = sqrt( 1 + $e2 / (1 - $e2) * pow( cos( $phy0 ), 4.0 ) ); |
$this->b0 = asin( $sinPhy0 / $this->alpha ); |
$this->K = log( tan( $PI / 4.0 + $this->b0 / 2.0 ) ) |
- $this->alpha |
* log( tan( $PI / 4.0 + $phy0 / 2.0 ) ) |
+ $this->alpha |
* $e / 2 |
* log( (1 + $e * $sinPhy0) |
/ (1 - $e * $sinPhy0) ); |
} |
/** |
* |
* @param type $p |
* @return type |
*/ |
public function forward( $p ) { |
$Sa1 = log( tan( $PI / 4.0 - $p->y / 2.0 ) ); |
$Sa2 = $this->e / 2.0 |
* log( (1 + $this->e * sin( $p->y )) |
/ (1 - $this->e * sin( $p->y )) ); |
$S = -$this->alpha * ($Sa1 + $Sa2) + $this->K; |
// spheric latitude |
$b = 2.0 * (atan( exp( $S ) ) - proj4phpCommon::PI / 4.0); |
// spheric longitude |
$I = $this->alpha * ($p->x - $this->lambda0); |
// psoeudo equatorial rotation |
$rotI = atan( sin( $I ) |
/ (sin( $this->b0 ) * tan( $b ) + |
cos( $this->b0 ) * cos( $I )) ); |
$rotB = asin( cos( $this->b0 ) * sin( $b ) - |
sin( $this->b0 ) * cos( $b ) * cos( $I ) ); |
$p->y = $this->R / 2.0 |
* log( (1 + sin( $rotB )) / (1 - sin( $rotB )) ) |
+ $this->y0; |
$p->x = $this->R * $rotI + $this->x0; |
return $p; |
} |
/** |
* |
* @param type $p |
* @return type |
*/ |
public function inverse( $p ) { |
$Y = $p->x - $this->x0; |
$X = $p->y - $this->y0; |
$rotI = $Y / $this->R; |
$rotB = 2 * (atan( exp( $X / $this->R ) ) - $PI / 4.0); |
$b = asin( cos( $this->b0 ) * sin( $rotB ) |
+ sin( $this->b0 ) * cos( $rotB ) * cos( $rotI ) ); |
$I = atan( sin( $rotI ) |
/ (cos( $this->b0 ) * cos( $rotI ) - sin( $this->b0 ) |
* tan( $rotB )) ); |
$lambda = $this->lambda0 + $I / $this->alpha; |
$S = 0.0; |
$phy = $b; |
$prevPhy = -1000.0; |
$iteration = 0; |
while( abs( $phy - $prevPhy ) > 0.0000001 ) { |
if( ++$iteration > 20 ) { |
Proj4php::reportError( "omercFwdInfinity" ); |
return; |
} |
//S = log(tan(PI / 4.0 + phy / 2.0)); |
$S = 1.0 |
/ $this->alpha |
* (log( tan( $PI / 4.0 + $b / 2.0 ) ) - $this->K) |
+ $this->e |
* log( tan( $PI / 4.0 |
+ asin( $this->e * sin( $phy ) ) |
/ 2.0 ) ); |
$prevPhy = $phy; |
$phy = 2.0 * atan( exp( $S ) ) - $PI / 2.0; |
} |
$p->x = $lambda; |
$p->y = $phy; |
return $p; |
} |
} |
Proj4php::$proj['somerc'] = new Proj4phpProjSomerc(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/cea.php |
---|
New file |
0,0 → 1,97 |
<?php |
/******************************************************************************* |
NAME LAMBERT CYLINDRICAL EQUAL AREA |
PURPOSE: Transforms input longitude and latitude to Easting and |
Northing for the Lambert Cylindrical Equal Area projection. |
This class of projection includes the Behrmann and |
Gall-Peters Projections. The |
longitude and latitude must be in radians. The Easting |
and Northing values will be returned in meters. |
PROGRAMMER DATE |
---------- ---- |
R. Marsden August 2009 |
Winwaed Software Tech LLC, http://www.winwaed.com |
This function was adapted from the Miller Cylindrical Projection in the Proj4php |
library. |
Note: This implementation assumes a Spherical Earth. The (commented) code |
has been included for the ellipsoidal forward transform, but derivation of |
the ellispoidal inverse transform is beyond me. Note that most of the |
Proj4php implementations do NOT currently support ellipsoidal figures. |
Therefore this is not seen as a problem - especially this lack of support |
is explicitly stated here. |
ALGORITHM REFERENCES |
1. "Cartographic Projection Procedures for the UNIX Environment - |
A User's Manual" by Gerald I. Evenden, USGS Open File Report 90-284 |
and Release 4 Interim Reports (2003) |
2. Snyder, John P., "Flattening the Earth - Two Thousand Years of Map |
Projections", Univ. Chicago Press, 1993 |
****************************************************************************** */ |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
class Proj4phpProjCea { |
/* Initialize the Cylindrical Equal Area projection |
------------------------------------------- */ |
public function init() { |
//no-op |
} |
/* Cylindrical Equal Area forward equations--mapping lat,long to x,y |
------------------------------------------------------------ */ |
public function forward( $p ) { |
$lon = $p->x; |
$lat = $p->y; |
/* Forward equations |
----------------- */ |
$dlon = Proj4php::$common->adjust_lon( $lon - $this->long0 ); |
$x = $this->x0 + $this->a * $dlon * cos( $this->lat_ts ); |
$y = $this->y0 + $this->a * sin( $lat ) / cos( $this->lat_ts ); |
/* Elliptical Forward Transform |
Not implemented due to a lack of a matchign inverse function |
{ |
$Sin_Lat = sin(lat); |
$Rn = $this->a * (sqrt(1.0e0 - $this->es * Sin_Lat * Sin_Lat )); |
x = $this->x0 + $this->a * dlon * cos($this->lat_ts); |
y = $this->y0 + Rn * sin(lat) / cos($this->lat_ts); |
} |
*/ |
$p->x = $x; |
$p->y = $y; |
return $p; |
} |
/** |
* Cylindrical Equal Area inverse equations--mapping x,y to lat/long |
* |
* @param type $p |
* @return type |
*/ |
public function inverse( $p ) { |
$p->x -= $this->x0; |
$p->y -= $this->y0; |
$p->x = Proj4php::$common->adjust_lon( $this->long0 + ($p->x / $this->a) / cos( $this->lat_ts ) ); |
$p->y = asin( ($p->y / $this->a) * cos( $this->lat_ts ) ); |
return $p; |
} |
} |
Proj4php::$proj['cea'] = new Proj4phpProjCea(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/cass.php |
---|
New file |
0,0 → 1,118 |
<?php |
/******************************************************************************* |
NAME CASSINI |
PURPOSE: Transforms input longitude and latitude to Easting and |
Northing for the Cassini projection. The |
longitude and latitude must be in radians. The Easting |
and Northing values will be returned in meters. |
Ported from PROJ.4. |
ALGORITHM REFERENCES |
1. 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. |
2. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", |
U.S. Geological Survey Professional Paper 1453 , United State Government |
****************************************************************************** */ |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
//Proj4php.defs["EPSG:28191"] = "+proj=cass +lat_0=31.73409694444445 +lon_0=35.21208055555556 +x_0=170251.555 +y_0=126867.909 +a=6378300.789 +b=6356566.435 +towgs84=-275.722,94.7824,340.894,-8.001,-4.42,-11.821,1 +units=m +no_defs"; |
// Initialize the Cassini projection |
// ----------------------------------------------------------------- |
class Proj4phpProjCass { |
public function init() { |
if( !$this->sphere ) { |
$this->en = Proj4php::$common->pj_enfn( $this->es ); |
$this->m0 = Proj4php::$common->pj_mlfn( $this->lat0, sin( $this->lat0 ), cos( $this->lat0 ), $this->en ); |
} |
} |
protected $C1 = .16666666666666666666; |
protected $C2 = .00833333333333333333; |
protected $C3 = .04166666666666666666; |
protected $C4 = .33333333333333333333; |
protected $C5 = .06666666666666666666; |
/* Cassini forward equations--mapping lat,long to x,y |
----------------------------------------------------------------------- */ |
public function forward( $p ) { |
/* Forward equations |
----------------- */ |
#$x; |
#$y; |
$lam = $p->x; |
$phi = $p->y; |
$lam = Proj4php::$common->adjust_lon( $lam - $this->long0 ); |
if( $this->sphere ) { |
$x = asin( cos( $phi ) * sin( $lam ) ); |
$y = atan2( tan( $phi ), cos( $lam ) ) - $this->phi0; |
} else { |
//ellipsoid |
$this->n = sin( $phi ); |
$this->c = cos( $phi ); |
$y = $this->pj_mlfn( $phi, $this->n, $this->c, $this->en ); |
$this->n = 1. / sqrt( 1. - $this->es * $this->n * $this->n ); |
$this->tn = 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); |
} |
$p->x = $this->a * $x + $this->x0; |
$p->y = $this->a * $y + $this->y0; |
return $p; |
} |
/* Inverse equations |
----------------- */ |
public function inverse( $p ) { |
$p->x -= $this->x0; |
$p->y -= $this->y0; |
$x = $p->x / $this->a; |
$y = $p->y / $this->a; |
if( $this->sphere ) { |
$this->dd = $y + $this->lat0; |
$phi = asin( sin( $this->dd ) * cos( $x ) ); |
$lam = atan2( tan( $x ), cos( $this->dd ) ); |
} else { |
/* ellipsoid */ |
$ph1 = Proj4php::$common->pj_inv_mlfn( $this->m0 + $y, $this->es, $this->en ); |
$this->tn = tan( $ph1 ); |
$this->t = $this->tn * $this->tn; |
$this->n = sin( $ph1 ); |
$this->r = 1. / (1. - $this->es * $this->n * $this->n); |
$this->n = 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)) / cos( $ph1 ); |
} |
$p->x = Proj4php::$common->adjust_lon( $this->long0 + $lam ); |
$p->y = $phi; |
return $p; |
} |
} |
Proj4php::$proj['cass'] = new Proj4phpProjCass(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/ortho.php |
---|
New file |
0,0 → 1,139 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/* * ***************************************************************************** |
NAME ORTHOGRAPHIC |
PURPOSE: Transforms input longitude and latitude to Easting and |
Northing for the Orthographic projection. The |
longitude and latitude must be in radians. The Easting |
and Northing values will be returned in meters. |
PROGRAMMER DATE |
---------- ---- |
T. Mittan Mar, 1993 |
ALGORITHM REFERENCES |
1. 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. |
2. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", |
U.S. Geological Survey Professional Paper 1453 , United State Government |
Printing Office, Washington D.C., 1989. |
* ***************************************************************************** */ |
class Proj4phpProjOrtho { |
/* Initialize the Orthographic projection |
------------------------------------- */ |
public function init( $def ) { |
//double temp; /* temporary variable */ |
/* Place parameters in static storage for common use |
------------------------------------------------- */; |
$this->sin_p14 = sin( $this->lat0 ); |
$this->cos_p14 = cos( $this->lat0 ); |
} |
/* Orthographic forward equations--mapping lat,long to x,y |
--------------------------------------------------- */ |
public function forward( $p ) { |
/* |
$sinphi; |
$cosphi; // sin and cos value |
$dlon; // delta longitude value |
$coslon; // cos of longitude |
$ksp; // scale factor |
$g; |
*/ |
$lon = $p->x; |
$lat = $p->y; |
/* Forward equations |
----------------- */ |
$dlon = Proj4php::$common->adjust_lon( $lon - $this->long0 ); |
$sinphi = sin( $lat ); |
$cosphi = cos( $lat ); |
$coslon = cos( $dlon ); |
$g = $this->sin_p14 * sinphi + $this->cos_p14 * $cosphi * $coslon; |
$ksp = 1.0; |
if( ($g > 0) || (abs( $g ) <= Proj4php::$common->EPSLN) ) { |
$x = $this->a * $ksp * $cosphi * sin( $dlon ); |
$y = $this->y0 + $this->a * $ksp * ($this->cos_p14 * $sinphi - $this->sin_p14 * $cosphi * $coslon); |
} else { |
Proj4php::reportError( "orthoFwdPointError" ); |
} |
$p->x = $x; |
$p->y = $y; |
return $p; |
} |
/** |
* |
* @param type $p |
* @return type |
*/ |
public function inverse( $p ) { |
/* |
$rh; // height above ellipsoid |
$z; // angle |
$sinz; |
$cosz; // sin of z and cos of z |
$temp; |
$con; |
$lon; |
$lat; |
*/ |
/* Inverse equations |
----------------- */ |
$p->x -= $this->x0; |
$p->y -= $this->y0; |
$rh = sqrt( $p->x * $p->x + $p->y * $p->y ); |
if( $rh > $this->a + .0000001 ) { |
Proj4php::reportError( "orthoInvDataError" ); |
} |
$z = Proj4php::$common . asinz( $rh / $this->a ); |
$sinz = sin( $z ); |
$cosz = cos( $z ); |
$lon = $this->long0; |
if( abs( $rh ) <= Proj4php::$common->EPSLN ) { |
$lat = $this->lat0; |
} |
$lat = Proj4php::$common . asinz( $cosz * $this->sin_p14 + ($p->y * $sinz * $this->cos_p14) / $rh ); |
$con = abs( $this->lat0 ) - Proj4php::$common->HALF_PI; |
if( abs( con ) <= Proj4php::$common->EPSLN ) { |
if( $this->lat0 >= 0 ) { |
$lon = Proj4php::$common->adjust_lon( $this->long0 + atan2( $p->x, -$p->y ) ); |
} else { |
$lon = Proj4php::$common->adjust_lon( $this->long0 - atan2( -$p->x, $p->y ) ); |
} |
} |
$con = $cosz - $this->sin_p14 * sin( $lat ); |
$p->x = $lon; |
$p->y = $lat; |
return $p; |
} |
} |
Proj4php::$proj['ortho'] = new Proj4phpProjOrtho(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/mill.php |
---|
New file |
0,0 → 1,82 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/******************************************************************************* |
NAME MILLER CYLINDRICAL |
PURPOSE: Transforms input longitude and latitude to Easting and |
Northing for the Miller Cylindrical projection. The |
longitude and latitude must be in radians. The Easting |
and Northing values will be returned in meters. |
PROGRAMMER DATE |
---------- ---- |
T. Mittan March, 1993 |
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. |
* ***************************************************************************** */ |
class Proj4phpProjMill { |
/* Initialize the Miller Cylindrical projection |
------------------------------------------- */ |
public function init() { |
//no-op |
} |
/* Miller Cylindrical forward equations--mapping lat,long to x,y |
------------------------------------------------------------ */ |
public function forward( $p ) { |
$lon = $p->x; |
$lat = $p->y; |
/* Forward equations |
----------------- */ |
$dlon = Proj4php::$common->adjust_lon( $lon - $this->long0 ); |
$x = $this->x0 + $this->a * $dlon; |
$y = $this->y0 + $this->a * log( tan( (Proj4php::$common->PI / 4.0) + ($lat / 2.5) ) ) * 1.25; |
$p->x = $x; |
$p->y = $y; |
return $p; |
} |
/* Miller Cylindrical inverse equations--mapping x,y to lat/long |
------------------------------------------------------------ */ |
public function inverse( $p ) { |
$p->x -= $this->x0; |
$p->y -= $this->y0; |
$lon = Proj4php::$common->adjust_lon( $this->long0 + $p->x / $this->a ); |
$lat = 2.5 * (atan( exp( 0.8 * $p->y / $this->a ) ) - Proj4php::$common->PI / 4.0); |
$p->x = $lon; |
$p->y = $lat; |
return $p; |
} |
} |
Proj4php::$proj['mill'] = new Proj4phpProjMill(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/krovak.php |
---|
New file |
0,0 → 1,156 |
<?php |
/** |
NOTES: According to EPSG the full Krovak projection method should have |
the following parameters. Within PROJ.4 the azimuth, and pseudo |
standard parallel are hardcoded in the algorithm and can't be |
altered from outside. The others all have defaults to match the |
common usage with Krovak projection. |
lat_0 = latitude of centre of the projection |
lon_0 = longitude of centre of the projection |
* * = azimuth (true) of the centre line passing through the centre of the projection |
* * = latitude of pseudo standard parallel |
k = scale factor on the pseudo standard parallel |
x_0 = False Easting of the centre of the projection at the apex of the cone |
y_0 = False Northing of the centre of the projection at the apex of the cone |
**/ |
class Proj4phpProjKrovak { |
/** |
* |
*/ |
public function init() { |
/* we want Bessel as fixed ellipsoid */ |
$this->a = 6377397.155; |
$this->es = 0.006674372230614; |
$this->e = sqrt( $this->es ); |
/* if latitude of projection center is not set, use 49d30'N */ |
if( !$this->lat0 ) { |
$this->lat0 = 0.863937979737193; |
} |
if( !$this->long0 ) { |
$this->long0 = 0.7417649320975901 - 0.308341501185665; |
} |
/* if scale not set default to 0.9999 */ |
if( !$this->k0 ) { |
$this->k0 = 0.9999; |
} |
$this->s45 = 0.785398163397448; /* 45° */ |
$this->s90 = 2 * $this->s45; |
$this->fi0 = $this->lat0; /* Latitude of projection centre 49° 30' */ |
/* Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128, |
e2=0.006674372230614; |
*/ |
$this->e2 = $this->es; /* 0.006674372230614; */ |
$this->e = sqrt( $this->e2 ); |
$this->alfa = sqrt( 1. + ($this->e2 * pow( cos( $this->fi0 ), 4 )) / (1. - $this->e2) ); |
$this->uq = 1.04216856380474; /* DU(2, 59, 42, 42.69689) */ |
$this->u0 = asin( sin( $this->fi0 ) / $this->alfa ); |
$this->g = pow( (1. + $this->e * sin( $this->fi0 )) / (1. - $this->e * sin( $this->fi0 )), $this->alfa * $this->e / 2. ); |
$this->k = tan( $this->u0 / 2. + $this->s45 ) / pow( tan( $this->fi0 / 2. + $this->s45 ), $this->alfa ) * $this->g; |
$this->k1 = $this->k0; |
$this->n0 = $this->a * sqrt( 1. - $this->e2 ) / (1. - $this->e2 * pow( sin( $this->fi0 ), 2 )); |
$this->s0 = 1.37008346281555; /* Latitude of pseudo standard parallel 78° 30'00" N */ |
$this->n = sin( $this->s0 ); |
$this->ro0 = $this->k1 * $this->n0 / tan( $this->s0 ); |
$this->ad = $this->s90 - $this->uq; |
} |
/** |
* ellipsoid |
* calculate xy from lat/lon |
* Constants, identical to inverse transform function |
* |
* @param type $p |
* @return type |
*/ |
public function forward( $p ) { |
$lon = $p->x; |
$lat = $p->y; |
$delta_lon = Proj4php::$common->adjust_lon( $lon - $this->long0 ); // Delta longitude |
/* Transformation */ |
$gfi = pow( ((1. + $this->e * sin( $lat )) / (1. - $this->e * sin( $lat )) ), ($this->alfa * $this->e / 2. ) ); |
$u = 2. * (atan( $this->k * pow( tan( $lat / 2. + $this->s45 ), $this->alfa ) / $gfi ) - $this->s45); |
$deltav = - $delta_lon * $this->alfa; |
$s = asin( cos( $this->ad ) * sin( $u ) + sin( $this->ad ) * cos( $u ) * cos( $deltav ) ); |
$d = asin( cos( $u ) * sin( $deltav ) / cos( $s ) ); |
$eps = $this->n * $d; |
$ro = $this->ro0 * pow( tan( $this->s0 / 2. + $this->s45 ), $this->n ) / pow( tan( $s / 2. + $this->s45 ), $this->n ); |
/* x and y are reverted! */ |
//$p->y = ro * cos(eps) / a; |
//$p->x = ro * sin(eps) / a; |
$p->y = $ro * cos( $eps ) / 1.0; |
$p->x = $ro * sin( $eps ) / 1.0; |
if( $this->czech ) { |
$p->y *= -1.0; |
$p->x *= -1.0; |
} |
return $p; |
} |
/** |
* calculate lat/lon from xy |
* |
* @param Point $p |
* @return Point $p |
*/ |
public function inverse( $p ) { |
/* Transformation */ |
/* revert y, x */ |
$tmp = $p->x; |
$p->x = $p->y; |
$p->y = $tmp; |
if( $this->czech ) { |
$p->y *= -1.0; |
$p->x *= -1.0; |
} |
$ro = sqrt( $p->x * $p->x + $p->y * $p->y ); |
$eps = atan2( $p->y, $p->x ); |
$d = $eps / sin( $this->s0 ); |
$s = 2. * (atan( pow( $this->ro0 / $ro, 1. / $this->n ) * tan( $this->s0 / 2. + $this->s45 ) ) - $this->s45); |
$u = asin( cos( $this->ad ) * sin( s ) - sin( $this->ad ) * cos( s ) * cos( d ) ); |
$deltav = asin( cos( $s ) * sin( $d ) / cos( $u ) ); |
$p->x = $this->long0 - $deltav / $this->alfa; |
/* ITERATION FOR $lat */ |
$fi1 = $u; |
$ok = 0; |
$iter = 0; |
do { |
$p->y = 2. * ( atan( pow( $this->k, -1. / $this->alfa ) * |
pow( tan( $u / 2. + $this->s45 ), 1. / $this->alfa ) * |
pow( (1. + $this->e * sin( $fi1 )) / (1. - $this->e * sin( $fi1 )), $this->e / 2. ) |
) - $this->s45); |
if( abs( $fi1 - $p->y ) < 0.0000000001 ) |
$ok = 1; |
$fi1 = $p->y; |
$iter += 1; |
} while( $ok == 0 && $iter < 15 ); |
if( $iter >= 15 ) { |
Proj4php::reportError( "PHI3Z-CONV:Latitude failed to converge after 15 iterations" ); |
//console.log('iter:', iter); |
return null; |
} |
return $p; |
} |
} |
Proj4php::$proj['krovak'] = new Proj4phpProjKrovak(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/vandg.php |
---|
New file |
0,0 → 1,167 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/******************************************************************************* |
NAME VAN DER GRINTEN |
PURPOSE: Transforms input Easting and Northing to longitude and |
latitude for the Van der Grinten projection. The |
Easting and Northing must be in meters. The longitude |
and latitude values will be returned in radians. |
PROGRAMMER DATE |
---------- ---- |
T. Mittan March, 1993 |
This function was adapted from the Van Der Grinten 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. |
* ***************************************************************************** */ |
class Proj4phpProjVandg { |
/* Initialize the Van Der Grinten projection |
---------------------------------------- */ |
public function init() { |
$this->R = 6370997.0; //Radius of earth |
} |
/** |
* |
* @param type $p |
* @return type |
*/ |
public function forward( $p ) { |
$lon = $p->x; |
$lat = $p->y; |
/* Forward equations |
----------------- */ |
$dlon = Proj4php::$common->adjust_lon( $lon - $this->long0 ); |
$x; |
$y; |
if( abs( $lat ) <= Proj4php::$common->EPSLN ) { |
$x = $this->x0 + $this->R * $dlon; |
$y = $this->y0; |
} |
$theta = Proj4php::$common . asinz( 2.0 * abs( $lat / Proj4php::$common->PI ) ); |
if( (abs( $dlon ) <= Proj4php::$common->EPSLN) || (abs( abs( $lat ) - Proj4php::$common->HALF_PI ) <= Proj4php::$common->EPSLN) ) { |
$x = $this->x0; |
if( $lat >= 0 ) { |
$y = $this->y0 + Proj4php::$common->PI * $this->R * tan( .5 * $theta ); |
} else { |
$y = $this->y0 + Proj4php::$common->PI * $this->R * - tan( .5 * $theta ); |
} |
// return(OK); |
} |
$al = .5 * abs( (Proj4php::$common->PI / $dlon) - ($dlon / Proj4php::$common->PI) ); |
$asq = $al * $al; |
$sinth = sin( $theta ); |
$costh = cos( $theta ); |
$g = $costh / ($sinth + $costh - 1.0); |
$gsq = $g * $g; |
$m = $g * (2.0 / $sinth - 1.0); |
$msq = $m * $m; |
$con = Proj4php::$common->PI * $this->R * ($al * ($g - $msq) + sqrt( $asq * ($g - $sq) * ($g - $msq) - ($msq + $asq) * ($gsq - $msq) )) / ($msq + $asq); |
if( $dlon < 0 ) { |
$con = -$con; |
} |
$x = $this->x0 + $con; |
$con = abs( $con / (Proj4php::$common->PI * $this->R) ); |
if( $lat >= 0 ) { |
$y = $this->y0 + Proj4php::$common->PI * $this->R * sqrt( 1.0 - $con * $con - 2.0 * $al * $con ); |
} else { |
$y = $this->y0 - Proj4php::$common->PI * $this->R * sqrt( 1.0 - $con * $con - 2.0 * $al * $con ); |
} |
$p->x = $x; |
$p->y = $y; |
return $p; |
} |
/* Van Der Grinten inverse equations--mapping x,y to lat/long |
--------------------------------------------------------- */ |
public function inverse( $p ) { |
/* |
$dlon; |
$xx; |
$yy; |
$xys; |
$c1; |
$c2; |
$c3; |
$al; |
$asq; |
$a1; |
$m1; |
$con; |
$th1; |
$d; |
*/ |
/* inverse equations |
----------------- */ |
$p->x -= $this->x0; |
$p->y -= $this->y0; |
$con = Proj4php::$common->PI * $this->R; |
$xx = $p->x / $con; |
$yy = $p->y / $con; |
$xys = $xx * $xx + $yy * $yy; |
$c1 = -abs( $yy ) * (1.0 + $xys); |
$c2 = $c1 - 2.0 * $yy * $yy + $xx * $xx; |
$c3 = -2.0 * $c1 + 1.0 + 2.0 * $yy * $yy + $xys * $xys; |
$d = $yy * $yy / $c3 + (2.0 * $c2 * $c2 * $c2 / $c3 / $c3 / $c3 - 9.0 * $c1 * $c2 / $c3 / $c3) / 27.0; |
$a1 = ($c1 - $c2 * $c2 / 3.0 / $c3) / $c3; |
$m1 = 2.0 * sqrt( -$a1 / 3.0 ); |
$con = ((3.0 * $d) / $a1) / $m1; |
if( abs( $con ) > 1.0 ) { |
if( $con >= 0.0 ) { |
$con = 1.0; |
} else { |
$con = -1.0; |
} |
} |
$th1 = acos( $con ) / 3.0; |
if( $p->$y >= 0 ) { |
$lat = (-$m1 * cos( $th1 + Proj4php::$common->PI / 3.0 ) - $c2 / 3.0 / $c3) * Proj4php::$common->PI; |
} else { |
$lat = -(-m1 * cos( $th1 + Proj4php::$common->PI / 3.0 ) - $c2 / 3.0 / $c3) * Proj4php::$common->PI; |
} |
if( abs( $xx ) < Proj4php::$common->EPSLN ) { |
$lon = $this->$long0; |
} |
$lon = Proj4php::$common->adjust_lon( $this->long0 + Proj4php::$common->PI * ($xys - 1.0 + sqrt( 1.0 + 2.0 * ($xx * $xx - $yy * $yy) + $xys * $xys )) / 2.0 / $xx ); |
$p->x = $lon; |
$p->y = $lat; |
return $p; |
} |
} |
Proj4php::$proj['vandg'] = new Proj4phpProjVandg(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/gnom.php |
---|
New file |
0,0 → 1,145 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/***************************************************************************** |
NAME GNOMONIC |
PURPOSE: Transforms input longitude and latitude to Easting and |
Northing for the Gnomonic Projection. |
Implementation based on the existing sterea and ortho |
implementations. |
PROGRAMMER DATE |
---------- ---- |
Richard Marsden November 2009 |
ALGORITHM REFERENCES |
1. Snyder, John P., "Flattening the Earth - Two Thousand Years of Map |
Projections", University of Chicago Press 1993 |
2. Wolfram Mathworld "Gnomonic Projection" |
http://mathworld.wolfram.com/GnomonicProjection.html |
Accessed: 12th November 2009 |
******************************************************************************/ |
class Proj4phpProjGnom { |
/** |
* Initialize the Gnomonic projection |
* |
* @todo $def not used in context...? |
* @param type $def |
*/ |
public function init( $def ) { |
/* Place parameters in static storage for common use |
------------------------------------------------- */ |
$this->sin_p14 = sin( $this->lat0 ); |
$this->cos_p14 = cos( $this->lat0 ); |
// Approximation for projecting points to the horizon (infinity) |
$this->infinity_dist = 1000 * $this->a; |
$this->rc = 1; |
} |
/* Gnomonic forward equations--mapping lat,long to x,y |
--------------------------------------------------- */ |
public function forward( $p ) { |
/* |
$sinphi; |
$cosphi; // sin and cos value |
$dlon; // delta longitude value |
$coslon; // cos of longitude |
$ksp; // scale factor |
$g; |
*/ |
$lon = $p->x; |
$lat = $p->y; |
/* Forward equations |
----------------- */ |
$dlon = Proj4php::$common->adjust_lon( $lon - $this->long0 ); |
$sinphi = sin( $lat ); |
$cosphi = cos( $lat ); |
$coslon = cos( $dlon ); |
$g = $this->sin_p14 * $sinphi + $this->cos_p14 * $cosphi * $coslon; |
$ksp = 1.0; |
if( (g > 0) || (abs( g ) <= Proj4php::$common->EPSLN) ) { |
$x = $this->x0 + $this->a * $ksp * $cosphi * sin( $dlon ) / $g; |
$y = $this->y0 + $this->a * $ksp * ($this->cos_p14 * $sinphi - $this->sin_p14 * $cosphi * $coslon) / $g; |
} else { |
Proj4php::reportError( "orthoFwdPointError" ); |
// Point is in the opposing hemisphere and is unprojectable |
// We still need to return a reasonable point, so we project |
// to infinity, on a bearing |
// equivalent to the northern hemisphere equivalent |
// This is a reasonable approximation for short shapes and lines that |
// straddle the horizon. |
$x = $this->x0 + $this->infinity_dist * $cosphi * sin( $dlon ); |
$y = $this->y0 + $this->infinity_dist * ($this->cos_p14 * $sinphi - $this->sin_p14 * $cosphi * $coslon); |
} |
$p->x = $x; |
$p->y = $y; |
return $p; |
} |
/** |
* |
* @param type $p |
* @return type |
*/ |
public function inverse( $p ) { |
/* |
$rh; // Rho |
$z; // angle |
$sinc; |
$cosc; |
$c; |
$lon; |
$lat; |
*/ |
/* Inverse equations |
----------------- */ |
$p->x = ($p->x - $this->x0) / $this->a; |
$p->y = ($p->y - $this->y0) / $this->a; |
$p->x /= $this->k0; |
$p->y /= $this->k0; |
if( ($rh = sqrt( $p->x * $p->x + $p->y * $p->y ) ) ) { |
$c = atan2( $rh, $this->rc ); |
$sinc = sin( $c ); |
$cosc = cos( $c ); |
$lat = Proj4php::$common->asinz( $cosc * $this->sin_p14 + ($p->y * $sinc * $this->cos_p14) / $rh ); |
$lon = atan2( $p->x * sinc, rh * $this->cos_p14 * $cosc - $p->y * $this->sin_p14 * $sinc ); |
$lon = Proj4php::$common->adjust_lon( $this->long0 + $lon ); |
} else { |
$lat = $this->phic0; |
$lon = 0.0; |
} |
$p->x = $lon; |
$p->y = $lat; |
return $p; |
} |
} |
Proj4php::$proj['gnom'] = new Proj4phpProjGnom(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/lcc.php |
---|
New file |
0,0 → 1,166 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/******************************************************************************* |
NAME LAMBERT CONFORMAL CONIC |
PURPOSE: Transforms input longitude and latitude to Easting and |
Northing for the Lambert Conformal Conic projection. The |
longitude and latitude must be in radians. The Easting |
and Northing values will be returned in meters. |
ALGORITHM REFERENCES |
1. 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. |
2. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", |
U.S. Geological Survey Professional Paper 1453 , United State Government |
*******************************************************************************/ |
//<2104> +proj=lcc +lat_1=10.16666666666667 +lat_0=10.16666666666667 +lon_0=-71.60561777777777 +k_0=1 +x0=-17044 +x0=-23139.97 +ellps=intl +units=m +no_defs no_defs |
// Initialize the Lambert Conformal conic projection |
// ----------------------------------------------------------------- |
//class Proj4phpProjlcc = Class.create(); |
class Proj4phpProjLcc { |
public function init() { |
// array of: r_maj,r_min,lat1,lat2,c_lon,c_lat,false_east,false_north |
//double c_lat; /* center latitude */ |
//double c_lon; /* center longitude */ |
//double lat1; /* first standard parallel */ |
//double lat2; /* second standard parallel */ |
//double r_maj; /* major axis */ |
//double r_min; /* minor axis */ |
//double false_east; /* x offset in meters */ |
//double false_north; /* y offset in meters */ |
//if lat2 is not defined |
if( !isset($this->lat2) ) { |
$this->lat2 = $this->lat0; |
} |
//if k0 is not defined |
if( !isset($this->k0) ) |
$this->k0 = 1.0; |
// Standard Parallels cannot be equal and on opposite sides of the equator |
if( abs( $this->lat1 + $this->lat2 ) < Proj4php::$common->EPSLN ) { |
Proj4php::reportError( "lcc:init: Equal Latitudes" ); |
return; |
} |
$temp = $this->b / $this->a; |
$this->e = sqrt( 1.0 - $temp * $temp ); |
$sin1 = sin( $this->lat1 ); |
$cos1 = cos( $this->lat1 ); |
$ms1 = Proj4php::$common->msfnz( $this->e, $sin1, $cos1 ); |
$ts1 = Proj4php::$common->tsfnz( $this->e, $this->lat1, $sin1 ); |
$sin2 = sin( $this->lat2 ); |
$cos2 = cos( $this->lat2 ); |
$ms2 = Proj4php::$common->msfnz( $this->e, $sin2, $cos2 ); |
$ts2 = Proj4php::$common->tsfnz( $this->e, $this->lat2, $sin2 ); |
$ts0 = Proj4php::$common->tsfnz( $this->e, $this->lat0, sin( $this->lat0 ) ); |
if( abs( $this->lat1 - $this->lat2 ) > Proj4php::$common->EPSLN ) { |
$this->ns = log( $ms1 / $ms2 ) / log( $ts1 / $ts2 ); |
} else { |
$this->ns = $sin1; |
} |
$this->f0 = $ms1 / ($this->ns * pow( $ts1, $this->ns )); |
$this->rh = $this->a * $this->f0 * pow( $ts0, $this->ns ); |
if( !isset($this->title) ) |
$this->title = "Lambert Conformal Conic"; |
} |
// Lambert Conformal conic forward equations--mapping lat,long to x,y |
// ----------------------------------------------------------------- |
public function forward( $p ) { |
$lon = $p->x; |
$lat = $p->y; |
// convert to radians |
if( $lat <= 90.0 && $lat >= -90.0 && $lon <= 180.0 && $lon >= -180.0 ) { |
//lon = lon * Proj4php::$common.D2R; |
//lat = lat * Proj4php::$common.D2R; |
} else { |
Proj4php::reportError( "lcc:forward: llInputOutOfRange: " . $lon . " : " . $lat ); |
return null; |
} |
$con = abs( abs( $lat ) - Proj4php::$common->HALF_PI ); |
if( $con > Proj4php::$common->EPSLN ) { |
$ts = Proj4php::$common->tsfnz( $this->e, $lat, sin( $lat ) ); |
$rh1 = $this->a * $this->f0 * pow( $ts, $this->ns ); |
} else { |
$con = $lat * $this->ns; |
if( $con <= 0 ) { |
Proj4php::reportError( "lcc:forward: No Projection" ); |
return null; |
} |
$rh1 = 0; |
} |
$theta = $this->ns * Proj4php::$common->adjust_lon( $lon - $this->long0 ); |
$p->x = $this->k0 * ($rh1 * sin( $theta )) + $this->x0; |
$p->y = $this->k0 * ($this->rh - $rh1 * cos( $theta )) + $this->y0; |
return $p; |
} |
/** |
* Lambert Conformal Conic inverse equations--mapping x,y to lat/long |
* |
* @param type $p |
* @return null |
*/ |
public function inverse( $p ) { |
$x = ($p->x - $this->x0) / $this->k0; |
$y = ($this->rh - ($p->y - $this->y0) / $this->k0); |
if( $this->ns > 0 ) { |
$rh1 = sqrt( $x * $x + $y * $y ); |
$con = 1.0; |
} else { |
$rh1 = -sqrt( $x * $x + $y * $y ); |
$con = -1.0; |
} |
$theta = 0.0; |
if( $rh1 != 0 ) { |
$theta = atan2( ($con * $x ), ($con * $y ) ); |
} |
if( ($rh1 != 0) || ($this->ns > 0.0) ) { |
$con = 1.0 / $this->ns; |
$ts = pow( ($rh1 / ($this->a * $this->f0) ), $con ); |
$lat = Proj4php::$common->phi2z( $this->e, $ts ); |
if( $lat == -9999 ) |
return null; |
} else { |
$lat = -Proj4php::$common->HALF_PI; |
} |
$lon = Proj4php::$common->adjust_lon( $theta / $this->ns + $this->long0 ); |
$p->x = $lon; |
$p->y = $lat; |
return $p; |
} |
} |
Proj4php::$proj['lcc'] = new Proj4phpProjLcc(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/laea.php |
---|
New file |
0,0 → 1,398 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/* * ***************************************************************************** |
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. |
* ***************************************************************************** */ |
class Proj4phpProjLaea { |
protected $S_POLE = 1; |
protected $N_POLE = 2; |
protected $EQUIT = 3; |
protected $OBLIQ = 4; |
protected $P00 = .33333333333333333333; |
protected $P01 = .17222222222222222222; |
protected $P02 = .10257936507936507936; |
protected $P10 = .06388888888888888888; |
protected $P11 = .06640211640211640211; |
protected $P20 = .01641501294219154443; |
/* Initialize the Lambert Azimuthal Equal Area projection |
------------------------------------------------------ */ |
public function init() { |
$t = abs( $this->lat0 ); |
if( abs( $t - Proj4php::$common->HALF_PI ) < Proj4php::$common->EPSLN ) { |
$this->mode = $this->lat0 < 0. ? $this->S_POLE : $this->N_POLE; |
} else if( abs( $t ) < Proj4php::$common->EPSLN ) { |
$this->mode = $this->EQUIT; |
} else { |
$this->mode = $this->OBLIQ; |
} |
if( $this->es > 0 ) { |
#$sinphi; |
$this->qp = Proj4php::$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 = sqrt( .5 * $this->qp ); |
$this->dd = 1. / $this->rq; |
$this->xmf = 1.; |
$this->ymf = .5 * $this->qp; |
break; |
case $this->OBLIQ: |
$this->rq = sqrt( .5 * $this->qp ); |
$sinphi = sin( $this->lat0 ); |
$this->sinb1 = Proj4php::$common->qsfnz( $this->e, $sinphi ) / $this->qp; |
$this->cosb1 = sqrt( 1. - $this->sinb1 * $this->sinb1 ); |
$this->dd = cos( $this->lat0 ) / (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 = sin( $this->lat0 ); |
$this->cosph0 = cos( $this->lat0 ); |
} |
} |
} |
/* Lambert Azimuthal Equal Area forward equations--mapping lat,long to x,y |
----------------------------------------------------------------------- */ |
public function forward( $p ) { |
/* Forward equations |
----------------- */ |
#$x; |
#$y; |
$lam = $p->x; |
$phi = $p->y; |
$lam = Proj4php::$common->adjust_lon( $lam - $this->long0 ); |
if( $this->sphere ) { |
/* |
$coslam; |
$cosphi; |
$sinphi; |
*/ |
$sinphi = sin( $phi ); |
$cosphi = cos( $phi ); |
$coslam = 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 <= Proj4php::$common->EPSLN ) { |
Proj4php::reportError( "laea:fwd:y less than eps" ); |
return null; |
} |
$y = sqrt( 2. / $y ); |
$x = $y * cosphi * 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( abs( $phi + $this->phi0 ) < Proj4php::$common->EPSLN ) { |
Proj4php::reportError( "laea:fwd:phi < eps" ); |
return null; |
} |
$y = Proj4php::$common->FORTPI - $phi * .5; |
$y = 2. * (($this->mode == $this->S_POLE) ? cos( $y ) : sin( $y )); |
$x = $y * sin( $lam ); |
$y *= $coslam; |
break; |
} |
} else { |
/* |
$coslam; |
$sinlam; |
$sinphi; |
$q; |
*/ |
$sinb = 0.0; |
$cosb = 0.0; |
$b = 0.0; |
$coslam = cos( $lam ); |
$sinlam = sin( $lam ); |
$sinphi = sin( $phi ); |
$q = Proj4php::$common->qsfnz( $this->e, $sinphi ); |
if( $this->mode == $this->OBLIQ || $this->mode == $this->EQUIT ) { |
$sinb = $q / $this->qp; |
$cosb = 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 = Proj4php::$common->HALF_PI + $phi; |
$q = $this->qp - $q; |
break; |
case $this->S_POLE: |
$b = $phi - Proj4php::$common->HALF_PI; |
$q = $this->qp + $q; |
break; |
} |
if( abs( $b ) < Proj4php::$common->EPSLN ) { |
Proj4php::reportError( "laea:fwd:b < eps" ); |
return null; |
} |
switch( $this->mode ) { |
case $this->OBLIQ: |
case $this->EQUIT: |
$b = sqrt( 2. / $b ); |
if( $this->mode == $this->OBLIQ ) { |
$y = $this->ymf * $b * ($this->cosb1 * $sinb - $this->sinb1 * $cosb * $coslam); |
} else { |
$y = ($b = 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 = sqrt( $q )) * $sinlam; |
$y = $coslam * (($this->mode == $this->S_POLE) ? $b : -$b); |
} else { |
$x = $y = 0.; |
} |
break; |
} |
} |
//v 1.0 |
/* |
$sin_lat=sin(lat); |
$cos_lat=cos(lat); |
$sin_delta_lon=sin(delta_lon); |
$cos_delta_lon=cos(delta_lon); |
$g =$this->sin_lat_o * sin_lat +$this->cos_lat_o * cos_lat * cos_delta_lon; |
if (g == -1.0) { |
Proj4php::reportError("laea:fwd:Point projects to a circle of radius "+ 2.0 * R); |
return null; |
} |
$ksp = $this->a * sqrt(2.0 / (1.0 + g)); |
$x = ksp * cos_lat * sin_delta_lon + $this->x0; |
$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; |
} |
/* Inverse equations |
----------------- */ |
public function inverse( $p ) { |
$p->x -= $this->x0; |
$p->y -= $this->y0; |
$x = $p->x / $this->a; |
$y = $p->y / $this->a; |
if( $this->sphere ) { |
$cosz = 0.0; |
#$rh; |
$sinz = 0.0; |
$rh = sqrt( $x * $x + $y * $y ); |
$phi = $rh * .5; |
if( $phi > 1. ) { |
Proj4php::reportError( "laea:Inv:DataError" ); |
return null; |
} |
$phi = 2. * asin( $phi ); |
if( $this->mode == $this->OBLIQ || $this->mode == $this->EQUIT ) { |
$sinz = sin( $phi ); |
$cosz = cos( $phi ); |
} |
switch( $this->mode ) { |
case $this->EQUIT: |
$phi = (abs( $rh ) <= Proj4php::$common->EPSLN) ? 0. : asin( $y * $sinz / $rh ); |
$x *= $sinz; |
$y = $cosz * $rh; |
break; |
case $this->OBLIQ: |
$phi = (abs( $rh ) <= Proj4php::$common->EPSLN) ? $this->phi0 : asin( $cosz * $this->sinph0 + $y * $sinz * $this->cosph0 / $rh ); |
$x *= $sinz * $this->cosph0; |
$y = ($cosz - sin( $phi ) * $this->sinph0) * $rh; |
break; |
case $this->N_POLE: |
$y = -$y; |
$phi = Proj4php::$common->HALF_PI - $phi; |
break; |
case $this->S_POLE: |
$phi -= Proj4php::$common->HALF_PI; |
break; |
} |
$lam = ($y == 0. && ($this->mode == $this->EQUIT || $this->mode == $this->OBLIQ)) ? 0. : atan2( $x, $y ); |
} else { |
/* |
$cCe; |
$sCe; |
$q; |
$rho; |
*/ |
$ab = 0.0; |
switch( $this->mode ) { |
case $this->EQUIT: |
case $this->OBLIQ: |
$x /= $this->dd; |
$y *= $this->dd; |
$rho = sqrt( $x * $x + $y * $y ); |
if( $rho < Proj4php::$common->EPSLN ) { |
$p->x = 0.; |
$p->y = $this->phi0; |
return $p; |
} |
$sCe = 2. * asin( .5 * $rho / $this->rq ); |
$cCe = cos( $sCe ); |
$x *= ($sCe = 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 = atan2( $x, $y ); |
$phi = $this->authlat( asin( $ab ), $this->apa ); |
} |
/* |
$Rh = sqrt($p->x *$p->x +$p->y * $p->y); |
$temp = Rh / (2.0 * $this->a); |
if (temp > 1) { |
Proj4php::reportError("laea:Inv:DataError"); |
return null; |
} |
$z = 2.0 * Proj4php::$common.asinz(temp); |
$sin_z=sin(z); |
$cos_z=cos(z); |
$lon =$this->long0; |
if (abs(Rh) > Proj4php::$common->EPSLN) { |
$lat = Proj4php::$common.asinz($this->sin_lat_o * cos_z +$this-> cos_lat_o * sin_z *$p->y / Rh); |
$temp =abs($this->lat0) - Proj4php::$common->HALF_PI; |
if (abs(temp) > Proj4php::$common->EPSLN) { |
temp = cos_z -$this->sin_lat_o * sin(lat); |
if(temp!=0.0) lon=Proj4php::$common->adjust_lon($this->long0+atan2($p->x*sin_z*$this->cos_lat_o,temp*Rh)); |
} else if ($this->lat0 < 0.0) { |
lon = Proj4php::$common->adjust_lon($this->long0 - atan2(-$p->x,$p->y)); |
} else { |
lon = Proj4php::$common->adjust_lon($this->long0 + atan2($p->x, -$p->y)); |
} |
} else { |
lat = $this->lat0; |
} |
*/ |
//return(OK); |
$p->x = Proj4php::$common->adjust_lon( $this->long0 + $lam ); |
$p->y = $phi; |
return $p; |
} |
/** |
* determine latitude from authalic latitude |
* |
* @param type $es |
* @return type |
*/ |
public function authset( $es ) { |
#$t; |
$APA = 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; |
} |
/** |
* |
* @param type $beta |
* @param type $APA |
* @return type |
*/ |
public function authlat( $beta, $APA ) { |
$t = $beta + $beta; |
return($beta + $APA[0] * sin( $t ) + $APA[1] * sin( $t + $t ) + $APA[2] * sin( $t + $t + $t )); |
} |
} |
Proj4php::$proj['laea'] = new Proj4phpProjLaea(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/moll.php |
---|
New file |
0,0 → 1,121 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/******************************************************************************* |
NAME MOLLWEIDE |
PURPOSE: Transforms input longitude and latitude to Easting and |
Northing for the MOllweide projection. The |
longitude and latitude must be in radians. The Easting |
and Northing values will be returned in meters. |
PROGRAMMER DATE |
---------- ---- |
D. Steinwand, EROS May, 1991; Updated Sept, 1992; Updated Feb, 1993 |
S. Nelson, EDC Jun, 2993; Made corrections in precision and |
number of iterations. |
ALGORITHM REFERENCES |
1. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", |
U.S. Geological Survey Professional Paper 1453 , United State Government |
Printing Office, Washington D.C., 1989. |
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. |
****************************************************************************** */ |
class Proj4phpProjMoll { |
/* Initialize the Mollweide projection |
------------------------------------ */ |
public function init() { |
//no-op |
} |
/* Mollweide forward equations--mapping lat,long to x,y |
---------------------------------------------------- */ |
public function forward( $p ) { |
/* Forward equations |
----------------- */ |
$lon = $p->x; |
$lat = $p->y; |
$delta_lon = Proj4php::$common->adjust_lon( $lon - $this->long0 ); |
$theta = $lat; |
$con = Proj4php::$common->PI * sin( $lat ); |
/* Iterate using the Newton-Raphson method to find theta |
----------------------------------------------------- */ |
for( $i = 0; true; ++$i ) { |
$delta_theta = -($theta + sin( $theta ) - $con) / (1.0 + cos( $theta )); |
$theta += $delta_theta; |
if( abs( $delta_theta ) < Proj4php::$common->EPSLN ) |
break; |
if( $i >= 50 ) { |
Proj4php::reportError( "moll:Fwd:IterationError" ); |
//return(241); |
} |
} |
$theta /= 2.0; |
/* If the latitude is 90 deg, force the x coordinate to be "0 . false easting" |
this is done here because of precision problems with "cos(theta)" |
-------------------------------------------------------------------------- */ |
if( Proj4php::$common->PI / 2 - abs( $lat ) < Proj4php::$common->EPSLN ) |
$delta_lon = 0; |
$x = 0.900316316158 * $this->a * $delta_lon * cos( $theta ) + $this->x0; |
$y = 1.4142135623731 * $this->a * sin( $theta ) + $this->y0; |
$p->x = $x; |
$p->y = $y; |
return $p; |
} |
/** |
* |
* @param type $p |
* @return type |
*/ |
public function inverse( $p ) { |
#$theta; |
#$arg; |
/* Inverse equations |
----------------- */ |
$p->x-= $this->x0; |
//~ $p->y -= $this->y0; |
$arg = $p->y / (1.4142135623731 * $this->a); |
/* Because of division by zero problems, 'arg' can not be 1.0. Therefore |
a number very close to one is used instead. |
------------------------------------------------------------------- */ |
if( abs( $arg ) > 0.999999999999 ) |
$arg = 0.999999999999; |
$theta = asin( $arg ); |
$lon = Proj4php::$common->adjust_lon( $this->long0 + ($p->x / (0.900316316158 * $this->a * cos( $theta ))) ); |
if( $lon < (-Proj4php::$common->PI) ) |
$lon = -Proj4php::$common->PI; |
if( $lon > Proj4php::$common->PI ) |
$lon = Proj4php::$common->PI; |
$arg = (2.0 * $theta + sin( 2.0 * $theta )) / Proj4php::$common->PI; |
if( abs( $arg ) > 1.0 ) |
$arg = 1.0; |
$lat = asin( $arg ); |
//return(OK); |
$p->x = $lon; |
$p->y = $lat; |
return $p; |
} |
} |
Proj4php::$proj['moll'] = new Proj4phpProjMoll(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/equi.php |
---|
New file |
0,0 → 1,80 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/******************************************************************************* |
NAME EQUIRECTANGULAR |
PURPOSE: Transforms input longitude and latitude to Easting and |
Northing for the Equirectangular projection. The |
longitude and latitude must be in radians. The Easting |
and Northing values will be returned in meters. |
PROGRAMMER DATE |
---------- ---- |
T. Mittan Mar, 1993 |
ALGORITHM REFERENCES |
1. 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. |
2. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", |
U.S. Geological Survey Professional Paper 1453 , United State Government |
Printing Office, Washington D.C., 1989. |
*******************************************************************************/ |
class Proj4phpProjEqui { |
public function init() { |
if( !$this->x0 ) |
$this->x0 = 0; |
if( !$this->y0 ) |
$this->y0 = 0; |
if( !$this->lat0 ) |
$this->lat0 = 0; |
if( !$this->long0 ) |
$this->long0 = 0; |
///$this->t2; |
} |
/* Equirectangular forward equations--mapping lat,long to x,y |
--------------------------------------------------------- */ |
public function forward( $p ) { |
$lon = $p->x; |
$lat = $p->y; |
$dlon = Proj4php::$common->adjust_lon( $lon - $this->long0 ); |
$x = $this->x0 + $this->a * $dlon * cos( $this->lat0 ); |
$y = $this->y0 + $this->a * $lat; |
$this->t1 = $x; |
$this->t2 = cos( $this->lat0 ); |
$p->x = $x; |
$p->y = $y; |
return $p; |
} |
/* Equirectangular inverse equations--mapping x,y to lat/long |
--------------------------------------------------------- */ |
public function inverse( $p ) { |
$p->x -= $this->x0; |
$p->y -= $this->y0; |
$lat = $p->y / $this->a; |
if( abs( $lat ) > Proj4php::$common->HALF_PI ) { |
Proj4php::reportError( "equi:Inv:DataError" ); |
} |
$lon = Proj4php::$common->adjust_lon( $this->long0 + $p->x / ($this->a * cos( $this->lat0 )) ); |
$p->x = $lon; |
$p->y = $lat; |
} |
} |
Proj4php::$proj['equi'] = new Proj4phpProjEqui(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/gstmerc.php |
---|
New file |
0,0 → 1,63 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
class Proj4phpProjGstmerc { |
public function init() { |
// array of: a, b, lon0, lat0, k0, x0, y0 |
$temp = $this->b / $this->a; |
$this->e = sqrt( 1.0 - $temp * $temp ); |
$this->lc = $this->long0; |
$this->rs = sqrt( 1.0 + $this->e * $this->e * pow( cos( $this->lat0 ), 4.0 ) / (1.0 - $this->e * $this->e) ); |
$sinz = sin( $this->lat0 ); |
$pc = asin( $sinz / $this->rs ); |
$sinzpc = sin( $pc ); |
$this->cp = Proj4php::$common->latiso( 0.0, $pc, $sinzpc ) - $this->rs * Proj4php::$common->latiso( $this->e, $this->lat0, $sinz ); |
$this->n2 = $this->k0 * $this->a * sqrt( 1.0 - $this->e * $this->e ) / (1.0 - $this->e * $this->e * $sinz * $sinz); |
$this->xs = $this->x0; |
$this->ys = $this->y0 - $this->n2 * $pc; |
if( !$this->title ) |
$this->title = "Gauss Schreiber transverse mercator"; |
} |
// forward equations--mapping lat,long to x,y |
// ----------------------------------------------------------------- |
public function forward( $p ) { |
$lon = $p->x; |
$lat = $p->y; |
$L = $this->rs * ($lon - $this->lc); |
$Ls = $this->cp + ($this->rs * Proj4php::$common->latiso( $this->e, $lat, sin( $lat ) )); |
$lat1 = asin( sin( $L ) / Proj4php::$common . cosh( $Ls ) ); |
$Ls1 = Proj4php::$common . latiso( 0.0, $lat1, sin( $lat1 ) ); |
$p->x = $this->xs + ($this->n2 * $Ls1); |
$p->y = $this->ys + ($this->n2 * atan( Proj4php::$common->sinh( $Ls ) / cos( $L ) )); |
return $p; |
} |
// inverse equations--mapping x,y to lat/long |
// ----------------------------------------------------------------- |
public function inverse( $p ) { |
$x = $p->x; |
$y = $p->y; |
$L = atan( Proj4php::$common . sinh( ($x - $this->xs) / $this->n2 ) / cos( ($y - $this->ys) / $this->n2 ) ); |
$lat1 = asin( sin( ($y - $this->ys) / $this->n2 ) / Proj4php::$common . cosh( ($x - $this->xs) / $this->n2 ) ); |
$LC = Proj4php::$common . latiso( 0.0, $lat1, sin( $lat1 ) ); |
$p->x = $this->lc + $L / $this->rs; |
$p->y = Proj4php::$common . invlatiso( $this->e, ($LC - $this->cp) / $this->rs ); |
return $p; |
} |
} |
Proj4php::$proj['gstmerc'] = new Proj4phpProjGestmerc(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/omerc.php |
---|
New file |
0,0 → 1,301 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/* * ***************************************************************************** |
NAME OBLIQUE MERCATOR (HOTINE) |
PURPOSE: Transforms input longitude and latitude to Easting and |
Northing for the Oblique Mercator projection. The |
longitude and latitude must be in radians. The Easting |
and Northing values will be returned in meters. |
PROGRAMMER DATE |
---------- ---- |
T. Mittan Mar, 1993 |
ALGORITHM REFERENCES |
1. 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. |
2. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", |
U.S. Geological Survey Professional Paper 1453 , United State Government |
Printing Office, Washington D.C., 1989. |
* ***************************************************************************** */ |
class Proj4phpProjOmerc { |
/* Initialize the Oblique Mercator projection |
------------------------------------------ */ |
public function init() { |
if( !$this->mode ) |
$this->mode = 0; |
if( !$this->lon1 ) { |
$this->lon1 = 0; |
$this->mode = 1; |
} |
if( !$this->lon2 ) |
$this->lon2 = 0; |
if( !$this->lat2 ) |
$this->lat2 = 0; |
/* Place parameters in static storage for common use |
------------------------------------------------- */ |
$temp = $this->b / $this->a; |
$es = 1.0 - pow( $temp, 2 ); |
$e = sqrt( $es ); |
$this->sin_p20 = sin( $this->lat0 ); |
$this->cos_p20 = cos( $this->lat0 ); |
$this->con = 1.0 - $this->es * $this->sin_p20 * $this->sin_p20; |
$this->com = sqrt( 1.0 - $es ); |
$this->bl = sqrt( 1.0 + $this->es * pow( $this->cos_p20, 4.0 ) / (1.0 - $es) ); |
$this->al = $this->a * $this->bl * $this->k0 * $this->com / $this->con; |
if( abs( $this->lat0 ) < Proj4php::$common->EPSLN ) { |
$this->ts = 1.0; |
$this->d = 1.0; |
$this->el = 1.0; |
} else { |
$this->ts = Proj4php::$common->tsfnz( $this->e, $this->lat0, $this->sin_p20 ); |
$this->con = sqrt( $this->con ); |
$this->d = $this->bl * $this->com / ($this->cos_p20 * $this->con); |
if( ($this->d * $this->d - 1.0) > 0.0 ) { |
if( $this->lat0 >= 0.0 ) { |
$this->f = $this->d + sqrt( $this->d * $this->d - 1.0 ); |
} else { |
$this->f = $this->d - sqrt( $this->d * $this->d - 1.0 ); |
} |
} else { |
$this->f = $this->d; |
} |
$this->el = $this->f * pow( $this->ts, $this->bl ); |
} |
//$this->longc=52.60353916666667; |
if( $this->mode != 0 ) { |
$this->g = .5 * ($this->f - 1.0 / $this->f); |
$this->gama = Proj4php::$common->asinz( sin( $this->alpha ) / $this->d ); |
$this->longc = $this->longc - Proj4php::$common->asinz( $this->g * tan( $this->gama ) ) / $this->bl; |
/* Report parameters common to format B |
------------------------------------- */ |
//genrpt(azimuth * R2D,"Azimuth of Central Line: "); |
//cenlon(lon_origin); |
// cenlat(lat_origin); |
$this->con = abs( $this->lat0 ); |
if( ($this->con > Proj4php::$common->EPSLN) && (abs( $this->con - Proj4php::$common->HALF_PI ) > Proj4php::$common->EPSLN) ) { |
$this->singam = sin( $this->gama ); |
$this->cosgam = cos( $this->gama ); |
$this->sinaz = sin( $this->alpha ); |
$this->cosaz = cos( $this->alpha ); |
if( $this->lat0 >= 0 ) { |
$this->u = ($this->al / $this->bl) * atan( sqrt( $this->d * $this->d - 1.0 ) / $this->cosaz ); |
} else { |
$this->u = -($this->al / $this->bl) * atan( sqrt( $this->d * $this->d - 1.0 ) / $this->cosaz ); |
} |
} else { |
Proj4php::reportError( "omerc:Init:DataError" ); |
} |
} else { |
$this->sinphi = sin( $this->at1 ); |
$this->ts1 = Proj4php::$common->tsfnz( $this->e, $this->lat1, $this->sinphi ); |
$this->sinphi = sin( $this->lat2 ); |
$this->ts2 = Proj4php::$common->tsfnz( $this->e, $this->lat2, $this->sinphi ); |
$this->h = pow( $this->ts1, $this->bl ); |
$this->l = pow( $this->ts2, $this->bl ); |
$this->f = $this->el / $this->h; |
$this->g = .5 * ($this->f - 1.0 / $this->f); |
$this->j = ($this->el * $this->el - $this->l * $this->h) / ($this->el * $this->el + $this->l * $this->h); |
$this->p = ($this->l - $this->h) / ($this->l + $this->h); |
$this->dlon = $this->lon1 - $this->lon2; |
if( $this->dlon < -Proj4php::$common->PI ) |
$this->lon2 = $this->lon2 - 2.0 * Proj4php::$common->PI; |
if( $this->dlon > Proj4php::$common->PI ) |
$this->lon2 = $this->lon2 + 2.0 * Proj4php::$common->PI; |
$this->dlon = $this->lon1 - $this->lon2; |
$this->longc = .5 * ($this->lon1 + $this->lon2) - atan( $this->j * tan( .5 * $this->bl * $this->dlon ) / $this->p ) / $this->bl; |
$this->dlon = Proj4php::$common->adjust_lon( $this->lon1 - $this->longc ); |
$this->gama = atan( sin( $this->bl * $this->dlon ) / $this->g ); |
$this->alpha = Proj4php::$common->asinz( $this->d * sin( $this->gama ) ); |
/* Report parameters common to format A |
------------------------------------- */ |
if( abs( $this->lat1 - $this->lat2 ) <= Proj4php::$common->EPSLN ) { |
Proj4php::reportError( "omercInitDataError" ); |
//return(202); |
} else { |
$this->con = abs( $this->lat1 ); |
} |
if( ($this->con <= Proj4php::$common->EPSLN) || (abs( $this->con - Proj4php::$common->HALF_PI ) <= Proj4php::$common->EPSLN) ) { |
Proj4php::reportError( "omercInitDataError" ); |
//return(202); |
} else { |
if( abs( abs( $this->lat0 ) - Proj4php::$common->HALF_PI ) <= Proj4php::$common->EPSLN ) { |
Proj4php::reportError( "omercInitDataError" ); |
//return(202); |
} |
} |
$this->singam = sin( $this->gam ); |
$this->cosgam = cos( $this->gam ); |
$this->sinaz = sin( $this->alpha ); |
$this->cosaz = cos( $this->alpha ); |
if( $this->lat0 >= 0 ) { |
$this->u = ($this->al / $this->bl) * atan( sqrt( $this->d * $this->d - 1.0 ) / $this->cosaz ); |
} else { |
$this->u = -($this->al / $this->bl) * atan( sqrt( $this->d * $this->d - 1.0 ) / $this->cosaz ); |
} |
} |
} |
/* Oblique Mercator forward equations--mapping lat,long to x,y |
---------------------------------------------------------- */ |
public function forward( $p ) { |
/* |
$theta; // angle |
$sin_phi; |
$cos_phi; // sin and cos value |
$b; // temporary values |
$c; |
$t; |
$tq; // temporary values |
$con; |
$n; |
$ml; // cone constant, small m |
$q; |
$us; |
$vl; |
$ul; |
$vs; |
$s; |
$dlon; |
$ts1; |
*/ |
$lon = $p->x; |
$lat = $p->y; |
/* Forward equations |
----------------- */ |
$sin_phi = sin( $lat ); |
$dlon = Proj4php::$common->adjust_lon( $lon - $this->longc ); |
$vl = sin( $this->bl * $dlon ); |
if( abs( abs( $lat ) - Proj4php::$common->HALF_PI ) > Proj4php::$common->EPSLN ) { |
$ts1 = Proj4php::$common->tsfnz( $this->e, $lat, $sin_phi ); |
$q = $this->el / (pow( $ts1, $this->bl )); |
$s = .5 * ($q - 1.0 / $q); |
$t = .5 * ($q + 1.0 / $q); |
$ul = ($s * $this->singam - $vl * $this->cosgam) / $t; |
$con = cos( $this->bl * $dlon ); |
if( abs( con ) < .0000001 ) { |
$us = $this->al * $this->bl * $dlon; |
} else { |
$us = $this->al * atan( ($s * $this->cosgam + $vl * $this->singam) / $con ) / $this->bl; |
if( $con < 0 ) |
$us = $us + Proj4php::$common->PI * $this->al / $this->bl; |
} |
} else { |
if( $lat >= 0 ) { |
$ul = $this->singam; |
} else { |
$ul = -$this->singam; |
} |
$us = $this->al * $lat / $this->bl; |
} |
if( abs( abs( $ul ) - 1.0 ) <= Proj4php::$common->EPSLN ) { |
//alert("Point projects into infinity","omer-for"); |
Proj4php::reportError( "omercFwdInfinity" ); |
//return(205); |
} |
$vs = .5 * $this->al * log( (1.0 - $ul) / (1.0 + $ul) ) / $this->bl; |
$us = $us - $this->u; |
$p->x = $this->x0 + $vs * $this->cosaz + $us * $this->sinaz; |
$p->y = $this->y0 + $us * $this->cosaz - $vs * $this->sinaz; |
return $p; |
} |
/** |
* |
* @param type $p |
* @return type |
*/ |
public function inverse( $p ) { |
/* |
$delta_lon; /* Delta longitude (Given longitude - center |
$theta; /* angle |
$delta_theta; /* adjusted longitude |
$sin_phi; |
$cos_phi; /* sin and cos value |
$b; /* temporary values |
$c; |
$t; |
$tq; /* temporary values |
$con; |
$n; |
$ml; /* cone constant, small m |
$vs; |
$us; |
$q; |
$s; |
$ts1; |
$vl; |
$ul; |
$bs; |
$dlon; |
$flag; |
*/ |
/* Inverse equations |
----------------- */ |
$p->x -= $this->x0; |
$p->y -= $this->y0; |
#$flag = 0; |
$vs = $p->x * $this->cosaz - $p->y * $this->sinaz; |
$us = $p->y * $this->cosaz + $p->x * $this->sinaz; |
$us = $us + $this->u; |
$q = exp( -$this->bl * $vs / $this->al ); |
$s = .5 * ($q - 1.0 / $q); |
$t = .5 * ($q + 1.0 / $q); |
$vl = sin( $this->bl * $us / $this->al ); |
$ul = ($vl * $this->cosgam + $s * $this->singam) / $t; |
if( abs( abs( $ul ) - 1.0 ) <= Proj4php::$common->EPSLN ) { |
$lon = $this->longc; |
if( ul >= 0.0 ) { |
$lat = Proj4php::$common->HALF_PI; |
} else { |
$lat = -Proj4php::$common->HALF_PI; |
} |
} else { |
$con = 1.0 / $this->bl; |
$ts1 = pow( ($this->el / sqrt( (1.0 + $ul) / (1.0 - $ul) ) ), $con ); |
$lat = Proj4php::$common->phi2z( $this->e, $ts1 ); |
//if ($flag != 0) |
//return($flag); |
//~ con = cos($this->bl * us /al); |
$theta = $this->longc - atan2( ($s * $this->cosgam - $vl * $this->singam ), $con ) / $this->bl; |
$lon = Proj4php::$common->adjust_lon( $theta ); |
} |
$p->x = $lon; |
$p->y = $lat; |
return $p; |
} |
} |
Proj4php::$proj['omerc'] = new Proj4phpProjOmerc(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/utm.php |
---|
New file |
0,0 → 1,74 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/******************************************************************************* |
NAME TRANSVERSE MERCATOR |
PURPOSE: Transforms input longitude and latitude to Easting and |
Northing for the Transverse Mercator projection. The |
longitude and latitude must be in radians. The Easting |
and Northing values will be returned in meters. |
ALGORITHM REFERENCES |
1. 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. |
2. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", |
U.S. Geological Survey Professional Paper 1453 , United State Government |
Printing Office, Washington D.C., 1989. |
*******************************************************************************/ |
/** |
Initialize Transverse Mercator projection |
*/ |
class Proj4phpProjUtm { |
public $dependsOn = 'tmerc'; |
public $utmSouth = false; // UTM north/south |
/** |
* |
* @return void |
*/ |
public function init() { |
if( !isset($this->zone) ) { |
Proj4php::reportError( "utm:init: zone must be specified for UTM" ); |
return; |
} |
$this->lat0 = 0.0; |
$this->long0 = ((6 * abs( $this->zone )) - 183) * Proj4php::$common->D2R; |
$this->x0 = 500000.0; |
$this->y0 = $this->utmSouth ? 10000000.0 : 0.0; |
$this->k0 = 0.9996; |
} |
/** |
* |
* @param type $p |
* @return type |
*/ |
public function forward( $p ) { |
return Proj4php::$proj['tmerc']->forward( $p ); |
} |
/** |
* |
* @param type $p |
* @return type |
*/ |
public function inverse( $p ) { |
return Proj4php::$proj['tmerc']->inverse( $p ); |
} |
} |
Proj4php::$proj['utm'] = new Proj4phpProjUtm(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/eqc.php |
---|
New file |
0,0 → 1,58 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/* similar to equi.js FIXME proj4 uses eqc */ |
class Proj4phpProjEqc { |
public function init() { |
if( !$this->x0 ) |
$this->x0 = 0; |
if( !$this->y0 ) |
$this->y0 = 0; |
if( !$this->lat0 ) |
$this->lat0 = 0; |
if( !$this->long0 ) |
$this->long0 = 0; |
if( !$this->lat_ts ) |
$this->lat_ts = 0; |
if( !$this->title ) |
$this->title = "Equidistant Cylindrical (Plate Carre)"; |
$this->rc = cos( $this->lat_ts ); |
} |
// forward equations--mapping lat,long to x,y |
// ----------------------------------------------------------------- |
public function forward( $p ) { |
$lon = $p->x; |
$lat = $p->y; |
$dlon = Proj4php::$common->adjust_lon( $lon - $this->long0 ); |
$dlat = Proj4php::$common . adjust_lat( $lat - $this->lat0 ); |
$p->x = $this->x0 + ($this->a * $dlon * $this->rc); |
$p->y = $this->y0 + ($this->a * $dlat ); |
return $p; |
} |
// inverse equations--mapping x,y to lat/long |
// ----------------------------------------------------------------- |
public function inverse( $p ) { |
$x = $p->x; |
$y = $p->y; |
$p->x = Proj4php::$common->adjust_lon( $this->long0 + (($x - $this->x0) / ($this->a * $this->rc)) ); |
$p->y = Proj4php::$common->adjust_lat( $this->lat0 + (($y - $this->y0) / ($this->a )) ); |
return $p; |
} |
} |
Proj4php::$proj['eqc'] = new Proj4phpProjEqc(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/tmerc.php |
---|
New file |
0,0 → 1,166 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/******************************************************************************* |
NAME TRANSVERSE MERCATOR |
PURPOSE: Transforms input longitude and latitude to Easting and |
Northing for the Transverse Mercator projection. The |
longitude and latitude must be in radians. The Easting |
and Northing values will be returned in meters. |
ALGORITHM REFERENCES |
1. 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. |
2. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", |
U.S. Geological Survey Professional Paper 1453 , United State Government |
Printing Office, Washington D.C., 1989. |
*******************************************************************************/ |
/** |
Initialize Transverse Mercator projection |
*/ |
class Proj4phpProjTmerc { |
private $e0, $e1, $e2, $e3, $ml0; |
/** |
* |
*/ |
public function init() { |
$this->e0 = Proj4php::$common->e0fn( $this->es ); |
$this->e1 = Proj4php::$common->e1fn( $this->es ); |
$this->e2 = Proj4php::$common->e2fn( $this->es ); |
$this->e3 = Proj4php::$common->e3fn( $this->es ); |
$this->ml0 = $this->a * Proj4php::$common->mlfn( $this->e0, $this->e1, $this->e2, $this->e3, $this->lat0 ); |
} |
/** |
Transverse Mercator Forward - long/lat to x/y |
long/lat in radians |
*/ |
public function forward( $p ) { |
$lon = $p->x; |
$lat = $p->y; |
$delta_lon = Proj4php::$common->adjust_lon( $lon - $this->long0 ); // Delta longitude |
#$con = 0; // cone constant |
#$x = 0; |
#$y = 0; |
$sin_phi = sin( $lat ); |
$cos_phi = cos( $lat ); |
if( isset($this->sphere) && $this->sphere === true ) { /* spherical form */ |
$b = $cos_phi * sin( $delta_lon ); |
if( (abs( abs( $b ) - 1.0 )) < .0000000001 ) { |
Proj4php::reportError( "tmerc:forward: Point projects into infinity" ); |
return(93); |
} else { |
$x = .5 * $this->a * $this->k0 * log( (1.0 + $b) / (1.0 - $b) ); |
$con = acos( $cos_phi * cos( $delta_lon ) / sqrt( 1.0 - $b * $b ) ); |
if( $lat < 0 ) |
$con = - $con; |
$y = $this->a * $this->k0 * ($con - $this->lat0); |
} |
} else { |
$al = $cos_phi * $delta_lon; |
$als = pow( $al, 2 ); |
$c = $this->ep2 * pow( $cos_phi, 2 ); |
$tq = tan( $lat ); |
$t = pow( $tq, 2 ); |
$con = 1.0 - $this->es * pow( $sin_phi, 2 ); |
$n = $this->a / sqrt( $con ); |
$ml = $this->a * Proj4php::$common->mlfn( $this->e0, $this->e1, $this->e2, $this->e3, $lat ); |
$x = $this->k0 * $n * $al * (1.0 + $als / 6.0 * (1.0 - $t + $c + $als / 20.0 * (5.0 - 18.0 * $t + pow( $t, 2 ) + 72.0 * $c - 58.0 * $this->ep2))) + $this->x0; |
$y = $this->k0 * ($ml - $this->ml0 + $n * $tq * ($als * (0.5 + $als / 24.0 * (5.0 - $t + 9.0 * $c + 4.0 * pow( $c, 2 ) + $als / 30.0 * (61.0 - 58.0 * $t + pow( $t, 2 ) + 600.0 * $c - 330.0 * $this->ep2))))) + $this->y0; |
} |
$p->x = $x; |
$p->y = $y; |
return $p; |
} |
/** |
Transverse Mercator Inverse - x/y to long/lat |
*/ |
public function inverse( $p ) { |
#$phi; /* temporary angles */ |
#$delta_phi; /* difference between longitudes */ |
$max_iter = 6; /* maximun number of iterations */ |
if( isset($this->sphere) && $this->sphere === true ) { /* spherical form */ |
$f = exp( $p->x / ($this->a * $this->k0) ); |
$g = .5 * ($f - 1 / $f); |
$temp = $this->lat0 + $p->y / ($this->a * $this->k0); |
$h = cos( $temp ); |
$con = sqrt( (1.0 - $h * $h) / (1.0 + $g * $g) ); |
$lat = Proj4php::$common->asinz( $con ); |
if( $temp < 0 ) |
$lat = -$lat; |
if( ($g == 0) && ($h == 0) ) { |
$lon = $this->long0; |
} else { |
$lon = Proj4php::$common->adjust_lon( atan2( $g, $h ) + $this->long0 ); |
} |
} else { // ellipsoidal form |
$x = $p->x - $this->x0; |
$y = $p->y - $this->y0; |
$con = ($this->ml0 + $y / $this->k0) / $this->a; |
$phi = $con; |
for( $i = 0; true; $i++ ) { |
$delta_phi = (($con + $this->e1 * sin( 2.0 * $phi ) - $this->e2 * sin( 4.0 * $phi ) + $this->e3 * sin( 6.0 * $phi )) / $this->e0) - $phi; |
$phi += $delta_phi; |
if( abs( $delta_phi ) <= Proj4php::$common->EPSLN ) |
break; |
if( $i >= $max_iter ) { |
Proj4php::reportError( "tmerc:inverse: Latitude failed to converge" ); |
return(95); |
} |
} // for() |
if( abs( $phi ) < Proj4php::$common->HALF_PI ) { |
// sincos(phi, &sin_phi, &cos_phi); |
$sin_phi = sin( $phi ); |
$cos_phi = cos( $phi ); |
$tan_phi = tan( $phi ); |
$c = $this->ep2 * pow( $cos_phi, 2 ); |
$cs = pow( $c, 2 ); |
$t = pow( $tan_phi, 2 ); |
$ts = pow( $t, 2 ); |
$con = 1.0 - $this->es * pow( $sin_phi, 2 ); |
$n = $this->a / sqrt( $con ); |
$r = $n * (1.0 - $this->es) / $con; |
$d = $x / ($n * $this->k0); |
$ds = pow( $d, 2 ); |
$lat = $phi - ($n * $tan_phi * $ds / $r) * (0.5 - $ds / 24.0 * (5.0 + 3.0 * $t + 10.0 * $c - 4.0 * $cs - 9.0 * $this->ep2 - $ds / 30.0 * (61.0 + 90.0 * $t + 298.0 * $c + 45.0 * $ts - 252.0 * $this->ep2 - 3.0 * $cs))); |
$lon = Proj4php::$common->adjust_lon( $this->long0 + ($d * (1.0 - $ds / 6.0 * (1.0 + 2.0 * $t + $c - $ds / 20.0 * (5.0 - 2.0 * $c + 28.0 * $t - 3.0 * $cs + 8.0 * $this->ep2 + 24.0 * $ts))) / $cos_phi) ); |
} else { |
$lat = Proj4php::$common->HALF_PI * Proj4php::$common->sign( $y ); |
$lon = $this->long0; |
} |
} |
$p->x = $lon; |
$p->y = $lat; |
return $p; |
} |
} |
Proj4php::$proj['tmerc'] = new Proj4phpProjTmerc(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/aeqd.php |
---|
New file |
0,0 → 1,101 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
class Proj4phpProjAeqd { |
public function init() { |
$this->sin_p12 = sin( $this->lat0 ); |
$this->cos_p12 = cos( $this->lat0 ); |
} |
/** |
* |
* @param type $p |
* @return type |
*/ |
public function forward( $p ) { |
#$lon = $p->x; |
#$lat = $p->y; |
#$ksp; |
$sinphi = sin( $p->y ); |
$cosphi = cos( $p->y ); |
$dlon = Proj4php::$common->adjust_lon( lon - $this->long0 ); |
$coslon = cos( $dlon ); |
$g = $this->sin_p12 * $sinphi + $this->cos_p12 * $cosphi * $coslon; |
if( abs( abs( $g ) - 1.0 ) < Proj4php::$common->EPSLN ) { |
$ksp = 1.0; |
if( $g < 0.0 ) { |
Proj4php::reportError( "aeqd:Fwd:PointError" ); |
return; |
} |
} else { |
$z = acos( $g ); |
$ksp = $z / sin( $z ); |
} |
$p->x = $this->x0 + $this->a * $ksp * $cosphi * sin( $dlon ); |
$p->y = $this->y0 + $this->a * $ksp * ($this->cos_p12 * $sinphi - $this->sin_p12 * $cosphi * $coslon); |
return $p; |
} |
/** |
* |
* @param type $p |
* @return type |
*/ |
public function inverse( $p ) { |
$p->x -= $this->x0; |
$p->y -= $this->y0; |
$rh = sqrt( $p->x * $p->x + $p->y * $p->y ); |
if( $rh > (2.0 * Proj4php::$common->HALF_PI * $this->a) ) { |
Proj4php::reportError( "aeqdInvDataError" ); |
return; |
} |
$z = $rh / $this->a; |
$sinz = sin( $z ); |
$cosz = cos( $z ); |
$lon = $this->long0; |
#$lat; |
if( abs( $rh ) <= Proj4php::$common->EPSLN ) { |
$lat = $this->lat0; |
} else { |
$lat = Proj4php::$common->asinz( $cosz * $this->sin_p12 + ($p->y * $sinz * $this->cos_p12) / $rh ); |
$con = abs( $this->lat0 ) - Proj4php::$common->HALF_PI; |
if( abs( $con ) <= Proj4php::$common->EPSLN ) { |
if( $this->lat0 >= 0.0 ) { |
$lon = Proj4php::$common->adjust_lon( $this->long0 + atan2( $p->x, -$p->y ) ); |
} else { |
$lon = Proj4php::$common->adjust_lon( $this->long0 - atan2( -$p->x, $p->y ) ); |
} |
} else { |
$con = $cosz - $this->sin_p12 * sin( $lat ); |
if( (abs( $con ) < Proj4php::$common->EPSLN) && (abs( $p->x ) < Proj4php::$common->EPSLN) ) { |
//no-op, just keep the lon value as is |
} else { |
#$temp = atan2( ($p->x * $sinz * $this->cos_p12 ), ($con * $rh ) ); // $temp is unused !?! |
$lon = Proj4php::$common->adjust_lon( $this->long0 + atan2( ($p->x * $sinz * $this->cos_p12 ), ($con * $rh ) ) ); |
} |
} |
} |
$p->x = $lon; |
$p->y = $lat; |
return $p; |
} |
} |
Proj4php::$proj['aeqd'] = new Proj4phpProjAeqd(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/nzmg.php |
---|
New file |
0,0 → 1,338 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/******************************************************************************* |
NAME NEW ZEALAND MAP GRID |
PURPOSE: Transforms input longitude and latitude to Easting and |
Northing for the New Zealand Map Grid projection. The |
longitude and latitude must be in radians. The Easting |
and Northing values will be returned in meters. |
ALGORITHM REFERENCES |
1. Department of Land and Survey Technical Circular 1973/32 |
http://www.linz.govt.nz/docs/miscellaneous/nz-map-definition.pdf |
2. OSG Technical Report 4.1 |
http://www.linz.govt.nz/docs/miscellaneous/nzmg.pdf |
IMPLEMENTATION NOTES |
The two references use different symbols for the calculated values. This |
implementation uses the variable names similar to the symbols in reference [1]. |
The alogrithm uses different units for delta latitude and delta longitude. |
The delta latitude is assumed to be in units of seconds of arc x 10^-5. |
The delta longitude is the usual radians. Look out for these conversions. |
The algorithm is described using complex arithmetic. There were three |
options: |
* find and use a Javascript library for complex arithmetic |
* write my own complex library |
* expand the complex arithmetic by hand to simple arithmetic |
This implementation has expanded the complex multiplication operations |
into parallel simple arithmetic operations for the real and imaginary parts. |
The imaginary part is way over to the right of the display; this probably |
violates every coding standard in the world, but, to me, it makes it much |
more obvious what is going on. |
The following complex operations are used: |
- addition |
- multiplication |
- division |
- complex number raised to integer power |
- summation |
A summary of complex arithmetic operations: |
(from http://en.wikipedia.org/wiki/Complex_arithmetic) |
addition: (a + bi) + (c + di) = (a + c) + (b + d)i |
subtraction: (a + bi) - (c + di) = (a - c) + (b - d)i |
multiplication: (a + bi) x (c + di) = (ac - bd) + (bc + ad)i |
division: (a + bi) / (c + di) = [(ac + bd)/(cc + dd)] + [(bc - ad)/(cc + dd)]i |
The algorithm needs to calculate summations of simple and complex numbers. This is |
implemented using a for-loop, pre-loading the summed value to zero. |
The algorithm needs to calculate theta^2, theta^3, etc while doing a summation. |
There are three possible implementations: |
- use pow in the summation loop - except for complex numbers |
- precalculate the values before running the loop |
- calculate theta^n = theta^(n-1) * theta during the loop |
This implementation uses the third option for both real and complex arithmetic. |
For example |
psi_n = 1; |
sum = 0; |
for (n = 1; n <=6; n++) { |
psi_n1 = psi_n * psi; // calculate psi^(n+1) |
psi_n = psi_n1; |
sum = sum + A[n] * psi_n; |
} |
TEST VECTORS |
NZMG E, N: 2487100.638 6751049.719 metres |
NZGD49 long, lat: 172.739194 -34.444066 degrees |
NZMG E, N: 2486533.395 6077263.661 metres |
NZGD49 long, lat: 172.723106 -40.512409 degrees |
NZMG E, N: 2216746.425 5388508.765 metres |
NZGD49 long, lat: 169.172062 -46.651295 degrees |
Note that these test vectors convert from NZMG metres to lat/long referenced |
to NZGD49, not the more usual WGS84. The difference is about 70m N/S and about |
10m E/W. |
These test vectors are provided in reference [1]. Many more test |
vectors are available in |
http://www.linz.govt.nz/docs/topography/topographicdata/placenamesdatabase/nznamesmar08.zip |
which is a catalog of names on the 260-series maps. |
EPSG CODES |
NZMG EPSG:27200 |
NZGD49 EPSG:4272 |
http://spatialreference.org/ defines these as |
Proj4php.defs["EPSG:4272"] = "+proj=longlat +ellps=intl +datum=nzgd49 +no_defs "; |
Proj4php.defs["EPSG:27200"] = "+proj=nzmg +lat_0=-41 +lon_0=173 +x_0=2510000 +y_0=6023150 +ellps=intl +datum=nzgd49 +units=m +no_defs "; |
LICENSE |
Copyright: Stephen Irons 2008 |
Released under terms of the LGPL as per: http://www.gnu.org/copyleft/lesser.html |
* ***************************************************************************** */ |
/** |
Initialize New Zealand Map Grip projection |
*/ |
class Proj4phpProjNzmg { |
/** |
* iterations: Number of iterations to refine inverse transform. |
* 0 -> km accuracy |
* 1 -> m accuracy -- suitable for most mapping applications |
* 2 -> mm accuracy |
*/ |
protected $iterations = 1; |
/** |
* |
*/ |
public function init() { |
$this->A = array( ); |
$this->A[1] = +0.6399175073; |
$this->A[2] = -0.1358797613; |
$this->A[3] = +0.063294409; |
$this->A[4] = -0.02526853; |
$this->A[5] = +0.0117879; |
$this->A[6] = -0.0055161; |
$this->A[7] = +0.0026906; |
$this->A[8] = -0.001333; |
$this->A[9] = +0.00067; |
$this->A[10] = -0.00034; |
$this->B_re = array( ); |
$this->B_im = array( ); |
$this->B_re[1] = +0.7557853228; |
$this->B_im[1] = 0.0; |
$this->B_re[2] = +0.249204646; |
$this->B_im[2] = +0.003371507; |
$this->B_re[3] = -0.001541739; |
$this->B_im[3] = +0.041058560; |
$this->B_re[4] = -0.10162907; |
$this->B_im[4] = +0.01727609; |
$this->B_re[5] = -0.26623489; |
$this->B_im[5] = -0.36249218; |
$this->B_re[6] = -0.6870983; |
$this->B_im[6] = -1.1651967; |
$this->C_re = array( ); |
$this->C_im = array( ); |
$this->C_re[1] = +1.3231270439; |
$this->C_im[1] = 0.0; |
$this->C_re[2] = -0.577245789; |
$this->C_im[2] = -0.007809598; |
$this->C_re[3] = +0.508307513; |
$this->C_im[3] = -0.112208952; |
$this->C_re[4] = -0.15094762; |
$this->C_im[4] = +0.18200602; |
$this->C_re[5] = +1.01418179; |
$this->C_im[5] = +1.64497696; |
$this->C_re[6] = +1.9660549; |
$this->C_im[6] = +2.5127645; |
$this->D = array( ); |
$this->D[1] = +1.5627014243; |
$this->D[2] = +0.5185406398; |
$this->D[3] = -0.03333098; |
$this->D[4] = -0.1052906; |
$this->D[5] = -0.0368594; |
$this->D[6] = +0.007317; |
$this->D[7] = +0.01220; |
$this->D[8] = +0.00394; |
$this->D[9] = -0.0013; |
} |
/** |
New Zealand Map Grid Forward - long/lat to x/y |
long/lat in radians |
*/ |
public function forward( $p ) { |
$lon = $p->x; |
$lat = $p->y; |
$delta_lat = $lat - $this->lat0; |
$delta_lon = $lon - $this->long0; |
// 1. Calculate d_phi and d_psi ... // and d_lambda |
// For this algorithm, delta_latitude is in seconds of arc x 10-5, so we need to scale to those units. Longitude is radians. |
$d_phi = $delta_lat / Proj4php::$common->SEC_TO_RAD * 1E-5; |
$d_lambda = $delta_lon; |
$d_phi_n = 1; // d_phi^0 |
$d_psi = 0; |
for( $n = 1; $n <= 10; $n++ ) { |
$d_phi_n = $d_phi_n * $d_phi; |
$d_psi = $d_psi + $this->A[$n] * $d_phi_n; |
} |
// 2. Calculate theta |
$th_re = $d_psi; |
$th_im = $d_lambda; |
// 3. Calculate z |
$th_n_re = 1; |
$th_n_im = 0; // theta^0 |
#$th_n_re1; |
#$th_n_im1; |
$z_re = 0; |
$z_im = 0; |
for( $n = 1; $n <= 6; $n++ ) { |
$th_n_re1 = $th_n_re * $th_re - $th_n_im * $th_im; |
$th_n_im1 = $th_n_im * $th_re + $th_n_re * $th_im; |
$th_n_re = $th_n_re1; |
$th_n_im = $th_n_im1; |
$z_re = $z_re + $this->B_re[$n] * $th_n_re - $this->B_im[$n] * $th_n_im; |
$z_im = $z_im + $this->B_im[$n] * $th_n_re + $this->B_re[$n] * $th_n_im; |
} |
// 4. Calculate easting and northing |
$p->x = ($z_im * $this->a) + $this->x0; |
$p->y = ($z_re * $this->a) + $this->y0; |
return $p; |
} |
/** |
New Zealand Map Grid Inverse - x/y to long/lat |
*/ |
public function inverse( $p ) { |
$x = $p->x; |
$y = $p->y; |
$delta_x = $x - $this->x0; |
$delta_y = $y - $this->y0; |
// 1. Calculate z |
$z_re = $delta_y / $this->a; |
$z_im = $delta_x / $this->a; |
// 2a. Calculate theta - first approximation gives km accuracy |
$z_n_re = 1; |
$z_n_im = 0; // z^0 |
$z_n_re1; |
$z_n_im1; |
$th_re = 0; |
$th_im = 0; |
for( $n = 1; $n <= 6; $n++ ) { |
$z_n_re1 = $z_n_re * $z_re - $z_n_im * $z_im; |
$z_n_im1 = $z_n_im * $z_re + $z_n_re * $z_im; |
$z_n_re = $z_n_re1; |
$z_n_im = $z_n_im1; |
$th_re = $th_re + $this->C_re[$n] * $z_n_re - $this->C_im[$n] * $z_n_im; |
$th_im = $th_im + $this->C_im[$n] * $z_n_re + $this->C_re[$n] * $z_n_im; |
} |
// 2b. Iterate to refine the accuracy of the calculation |
// 0 iterations gives km accuracy |
// 1 iteration gives m accuracy -- good enough for most mapping applications |
// 2 iterations bives mm accuracy |
for( $i = 0; $i < $this->iterations; $i++ ) { |
$th_n_re = $th_re; |
$th_n_im = $th_im; |
$th_n_re1; |
$th_n_im1; |
$num_re = $z_re; |
$num_im = $z_im; |
for( $n = 2; $n <= 6; $n++ ) { |
$th_n_re1 = $th_n_re * th_re - $th_n_im * $th_im; |
$th_n_im1 = $th_n_im * $th_re + $th_n_re * $th_im; |
$th_n_re = $th_n_re1; |
$th_n_im = $th_n_im1; |
$num_re = $num_re + ($n - 1) * ($this->B_re[$n] * $th_n_re - $this->B_im[$n] * $th_n_im); |
$num_im = $num_im + (n - 1) * ($this->B_im[$n] * $th_n_re + $this->B_re[$n] * $th_n_im); |
} |
$th_n_re = 1; |
$th_n_im = 0; |
$den_re = $this->B_re[1]; |
$den_im = $this->B_im[1]; |
for( $n = 2; $n <= 6; $n++ ) { |
$th_n_re1 = $th_n_re * $th_re - $th_n_im * $th_im; |
$th_n_im1 = $th_n_im * $th_re + $th_n_re * $th_im; |
$th_n_re = $th_n_re1; |
$th_n_im = $th_n_im1; |
$den_re = $den_re + $n * ($this->B_re[$n] * $th_n_re - $this->B_im[$n] * $th_n_im); |
$den_im = $den_im + $n * ($this->B_im[n] * $th_n_re + $this->B_re[$n] * $th_n_im); |
} |
// Complex division |
$den2 = $den_re * $den_re + $den_im * $den_im; |
$th_re = ($num_re * $den_re + $num_im * $den_im) / $den2; |
$th_im = ($num_im * $den_re - $num_re * $den_im) / $den2; |
} |
// 3. Calculate d_phi ... // and d_lambda |
$d_psi = $th_re; |
$d_lambda = $th_im; |
$d_psi_n = 1; // d_psi^0 |
$d_phi = 0; |
for( $n = 1; $n <= 9; $n++ ) { |
$d_psi_n = $d_psi_n * $d_psi; |
$d_phi = $d_phi + $this->D[$n] * $d_psi_n; |
} |
// 4. Calculate latitude and longitude |
// d_phi is calcuated in second of arc * 10^-5, so we need to scale back to radians. d_lambda is in radians. |
$lat = $this->lat0 + ($d_phi * Proj4php::$common->SEC_TO_RAD * 1E5); |
$lon = $this->long0 + $d_lambda; |
$p->x = $lon; |
$p->y = $lat; |
return $p; |
} |
} |
Proj4php::$proj['nzmg'] = new Proj4phpProjNzmg(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/eqdc.php |
---|
New file |
0,0 → 1,152 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/******************************************************************************* |
NAME EQUIDISTANT CONIC |
PURPOSE: Transforms input longitude and latitude to Easting and Northing |
for the Equidistant Conic projection. The longitude and |
latitude must be in radians. The Easting and Northing values |
will be returned in meters. |
PROGRAMMER DATE |
---------- ---- |
T. Mittan Mar, 1993 |
ALGORITHM REFERENCES |
1. 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. |
2. Snyder, John $p-> and Voxland, Philip M., "An Album of Map Projections", |
U.S. Geological Survey Professional Paper 1453 , United State Government |
Printing Office, Washington D.C., 1989. |
*******************************************************************************/ |
/* Variables common to all subroutines in this code file |
-----------------------------------------------------*/ |
class Proj4phpProjEqdc { |
/* Initialize the Equidistant Conic projection |
------------------------------------------ */ |
public function init() { |
/* Place parameters in static storage for common use |
------------------------------------------------- */ |
if( !$this->mode ) |
$this->mode = 0; //chosen default mode |
$this->temp = $this->b / $this->a; |
$this->es = 1.0 - pow( $this->temp, 2 ); |
$this->e = sqrt( $this->es ); |
$this->e0 = Proj4php::$common->e0fn( $this->es ); |
$this->e1 = Proj4php::$common->e1fn( $this->es ); |
$this->e2 = Proj4php::$common->e2fn( $this->es ); |
$this->e3 = Proj4php::$common->e3fn( $this->es ); |
$this->sinphi = sin( $this->lat1 ); |
$this->cosphi = cos( $this->lat1 ); |
$this->ms1 = Proj4php::$common->msfnz( $this->e, $this->sinphi, $this->cosphi ); |
$this->ml1 = Proj4php::$common->mlfn( $this->e0, $this->e1, $this->e2, $this->e3, $this->lat1 ); |
/* format B |
--------- */ |
if( $this->mode != 0 ) { |
if( abs( $this->lat1 + $this->lat2 ) < Proj4php::$common->EPSLN ) { |
Proj4php::reportError( "eqdc:Init:EqualLatitudes" ); |
//return(81); |
} |
$this->sinphi = sin( $this->lat2 ); |
$this->cosphi = cos( $this->lat2 ); |
$this->ms2 = Proj4php::$common->msfnz( $this->e, $this->sinphi, $this->cosphi ); |
$this->ml2 = Proj4php::$common->mlfn( $this->e0, $this->e1, $this->e2, $this->e3, $this->lat2 ); |
if( abs( $this->lat1 - $this->lat2 ) >= Proj4php::$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 = Proj4php::$common->mlfn( $this->e0, $this->e1, $this->e2, $this->e3, $this->lat0 ); |
$this->rh = $this->a * ($this->g - $this->ml0); |
} |
/* Equidistant Conic forward equations--mapping lat,long to x,y |
----------------------------------------------------------- */ |
public function forward( $p ) { |
$lon = $p->x; |
$lat = $p->y; |
/* Forward equations |
----------------- */ |
$ml = Proj4php::$common->mlfn( $this->e0, $this->e1, $this->e2, $this->e3, $lat ); |
$rh1 = $this->a * ($this->g - $ml); |
$theta = $this->ns * Proj4php::$common->adjust_lon( $lon - $this->long0 ); |
$x = $this->x0 + $rh1 * sin( $theta ); |
$y = $this->y0 + $this->rh - $rh1 * cos( $theta ); |
$p->x = $x; |
$p->y = $y; |
return $p; |
} |
/* Inverse equations |
----------------- */ |
public function inverse( $p ) { |
$p->x -= $this->x0; |
$p->y = $this->rh - $p->y + $this->y0; |
if( $this->ns >= 0 ) { |
$rh1 = sqrt( $p->x * $p->x + $p->y * $p->y ); |
$con = 1.0; |
} else { |
$rh1 = -sqrt( $p->x * $p->x + $p->y * $p->y ); |
$con = -1.0; |
} |
$theta = 0.0; |
if( $rh1 != 0.0 ) |
$theta = atan2( $con * $p->x, $con * $p->y ); |
$ml = $this->g - $rh1 / $this->a; |
$lat = $this->phi3z( $ml, $this->e0, $this->e1, $this->e2, $this->e3 ); |
$lon = Proj4php::$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. |
----------------------------------------------------------------- */ |
public function phi3z( $ml, $e0, $e1, $e2, $e3 ) { |
$phi = $ml; |
for( $i = 0; $i < 15; $i++ ) { |
$dphi = ($ml + $e1 * sin( 2.0 * $phi ) - $e2 * sin( 4.0 * $phi ) + $e3 * sin( 6.0 * $phi )) / $e0 - $phi; |
$phi += $dphi; |
if( abs( $dphi ) <= .0000000001 ) { |
return $phi; |
} |
} |
Proj4php::reportError( "PHI3Z-CONV:Latitude failed to converge after 15 iterations" ); |
return null; |
} |
} |
Proj4php::$proj['eqdc'] = new Proj4phpProjEqdc(); |
/tags/v5.9-aulnaie/scripts/modules/ifn/bibliotheque/proj4php/projCode/sinu.php |
---|
New file |
0,0 → 1,139 |
<?php |
/** |
* Author : Julien Moquet |
* |
* Inspired by Proj4php from Mike Adair madairATdmsolutions.ca |
* and Richard Greenwood rich@greenwoodma$p->com |
* License: LGPL as per: http://www.gnu.org/copyleft/lesser.html |
*/ |
/******************************************************************************* |
NAME SINUSOIDAL |
PURPOSE: Transforms input longitude and latitude to Easting and |
Northing for the Sinusoidal projection. The |
longitude and latitude must be in radians. The Easting |
and Northing values will be returned in meters. |
PROGRAMMER DATE |
---------- ---- |
D. Steinwand, EROS May, 1991 |
This function was adapted from the Sinusoidal 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. 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. |
2. "Software Documentation for GCTP General Cartographic Transformation |
Package", U.S. Geological Survey National Mapping Division, May 1982. |
* ***************************************************************************** */ |
class Proj4phpProjSinu { |
/* Initialize the Sinusoidal projection |
------------------------------------ */ |
public function init() { |
/* Place parameters in static storage for common use |
------------------------------------------------- */ |
#$this->R = 6370997.0; //Radius of earth |
if( !$this->sphere ) { |
$this->en = Proj4php::$common->pj_enfn( $this->es ); |
} else { |
$this->n = 1.; |
$this->m = 0.; |
$this->es = 0; |
$this->C_y = sqrt( ($this->m + 1.) / $this->n ); |
$this->C_x = $this->C_y / ($this->m + 1.); |
} |
} |
/* Sinusoidal forward equations--mapping lat,long to x,y |
----------------------------------------------------- */ |
public function forward( $p ) { |
#$x,y,delta_lon; |
$lon = $p->x; |
$lat = $p->y; |
/* Forward equations |
----------------- */ |
$lon = Proj4php::$common->adjust_lon( $lon - $this->long0 ); |
if( isset($this->sphere) ) { |
if( !$this->m ) { |
$lat = $this->n != 1. ? asin( $this->n * sin( $lat ) ) : $lat; |
} else { |
$k = $this->n * sin( $lat ); |
for( $i = Proj4php::$common->MAX_ITER; $i; --$i ) { |
$V = ($this->m * $lat + sin( $lat ) - $k) / ($this->m + cos( $lat )); |
$lat -= $V; |
if( abs( $V ) < Proj4php::$common->EPSLN ) |
break; |
} |
} |
$x = $this->a * $this->C_x * $lon * ($this->m + cos( $lat )); |
$y = $this->a * $this->C_y * $lat; |
} else { |
$s = sin( $lat ); |
$c = cos( $lat ); |
$y = $this->a * Proj4php::$common->pj_mlfn( $lat, $s, $c, $this->en ); |
$x = $this->a * $lon * $c / sqrt( 1. - $this->es * $s * $s ); |
} |
$p->x = $x; |
$p->y = $y; |
return $p; |
} |
/** |
* |
* @param type $p |
* @return type |
*/ |
public function inverse( $p ) { |
#$lat; |
#$temp; |
#$lon; |
/* Inverse equations |
----------------- */ |
$p->x -= $this->x0; |
$p->y -= $this->y0; |
$lat = $p->y / $this->a; |
if( isset($this->sphere) ) { |
$p->y /= $this->C_y; |
$lat = $this->m ? asin( ($this->m * $p->y + sin( $p->y )) / $this->n ) : ( $this->n != 1. ? asin( sin( $p->y ) / $this->n ) : $p->y ); |
$lon = $p->x / ($this->C_x * ($this->m + cos( $p->y ))); |
} |
else { |
$lat = Proj4php::$common->pj_inv_mlfn( $p->y / $this->a, $this->es, $this->en ); |
$s = abs( $lat ); |
if( $s < Proj4php::$common->HALF_PI ) { |
$s = sin( $lat ); |
$temp = $this->long0 + $p->x * sqrt( 1. - $this->es * $s * $s ) / ($this->a * cos( $lat )); |
//temp = $this->long0 + $p->x / ($this->a * cos($lat)); |
$lon = Proj4php::$common->adjust_lon( $temp ); |
} else if( ($s - Proj4php::$common->EPSLN) < Proj4php::$common->HALF_PI ) { |
$lon = $this->long0; |
} |
} |
$p->x = $lon; |
$p->y = $lat; |
return $p; |
} |
} |
Proj4php::$proj['sinu'] = new Proj4phpProjSinu(); |