391 lines
13 KiB
PHP
391 lines
13 KiB
PHP
|
<?php
|
||
|
|
function supprDecimales($dec) {
|
||
|
|
if ($dec>0 )
|
||
|
|
return floor($dec);
|
||
|
|
else
|
||
|
|
return ceil($dec);
|
||
|
|
}
|
||
|
|
|
||
|
|
function dec2dms($dec) {
|
||
|
|
$d = supprDecimales($dec);
|
||
|
|
$m = supprDecimales(($dec - $d) * 60);
|
||
|
|
$s = abs(round(((($dec - $d) * 60) - $m) * 60));
|
||
|
|
$m = abs($m);
|
||
|
|
return $d.'°'.$m."'".$s.'"';
|
||
|
|
}
|
||
|
|
|
||
|
|
/**/
|
||
|
|
function ALG0001($phi,$e) {
|
||
|
|
$temp = ( 1 - ( $e * sin( $phi ) ) ) / ( 1 + ( $e * sin( $phi ) ) );
|
||
|
|
$L = log ( tan ( (pi()/4) + ($phi/2) ) * pow ($temp, ($e/2) ));
|
||
|
|
return $L;
|
||
|
|
}
|
||
|
|
|
||
|
|
/** Calcul de la latitude à partir de la latitude isométrique
|
||
|
|
**/
|
||
|
|
function ALG0002($L,$e,$epsilon) {
|
||
|
|
$phi[0] = 2 * atan(exp($L)) - (pi()/2);
|
||
|
|
|
||
|
|
$i=0;
|
||
|
|
do {
|
||
|
|
$i++;
|
||
|
|
$temp = ( 1 + ( $e * sin( $phi[$i-1] ) ) ) / ( 1 - ( $e * sin( $phi[$i-1] ) ) );
|
||
|
|
$phi[$i] = 2 * atan ( pow ($temp, ($e/2)) * exp ($L) ) - pi()/2;
|
||
|
|
}
|
||
|
|
while (abs($phi[$i] - $phi[$i - 1]) >= $epsilon);
|
||
|
|
|
||
|
|
return $phi[$i];
|
||
|
|
}
|
||
|
|
|
||
|
|
/** Transformation de coordonnées en projection conique conforme de Lambert, en coordonnées géographiques
|
||
|
|
** @param double $X Coordonnée X en projection conique conforme de Lambert du point
|
||
|
|
** @param double $Y Coordonnée Y en projection conique conforme de Lambert du point
|
||
|
|
** @param double $n Exposant de la projection
|
||
|
|
** @param double $c Constante de la projection
|
||
|
|
** @param double $Xs Coordonnée Xs en projection du pôle
|
||
|
|
** @param double $Ys Coordonnée Ys en projection du pôle
|
||
|
|
** @param double $lambdac Longitude de l'origine par rapport au méridien origine
|
||
|
|
** @param double $e Première excentricité de l'ellipsoïde
|
||
|
|
** @param double $epsilon Tolérance de convergence
|
||
|
|
** @return array lambda Longitude par rapport au méridien origine
|
||
|
|
** phi Latitude
|
||
|
|
**/
|
||
|
|
function ALG0004($X,$Y,$n,$c,$Xs,$Ys,$lambdac,$e,$epsilon) {
|
||
|
|
$R = sqrt( pow(($X - $Xs),2) + pow(($Y - $Ys),2) );
|
||
|
|
$gamma = atan(($X - $Xs)/($Ys - $Y));
|
||
|
|
|
||
|
|
$lambda = $lambdac + ($gamma / $n);
|
||
|
|
|
||
|
|
$L = (-1 / $n) * log(abs($R/$c));
|
||
|
|
|
||
|
|
$phi = ALG0002($L,$e,$epsilon);
|
||
|
|
|
||
|
|
$coords['lambda']=$lambda;
|
||
|
|
$coords['phi']=$phi;
|
||
|
|
|
||
|
|
return $coords;
|
||
|
|
}
|
||
|
|
|
||
|
|
/** Transformation des coordonnées géographiques ellipsoïdales en coordonnées cartésiennes
|
||
|
|
**
|
||
|
|
**/
|
||
|
|
function ALG0009($lambda,$phi,$he,$a,$e) {
|
||
|
|
$N = ALG0021($phi,$a,$e);
|
||
|
|
|
||
|
|
$X = ($N + $he) * cos($phi) * cos($lambda);
|
||
|
|
|
||
|
|
$Y = ($N + $he) * cos($phi) * sin($lambda);
|
||
|
|
|
||
|
|
$Z = ($N * (1 - $e*$e) + $he) * sin ($phi);
|
||
|
|
|
||
|
|
$coords['X']=$X;
|
||
|
|
$coords['Y']=$Y;
|
||
|
|
$coords['Z']=$Z;
|
||
|
|
|
||
|
|
return $coords;
|
||
|
|
}
|
||
|
|
|
||
|
|
/**/
|
||
|
|
function ALG0012($X,$Y,$Z,$a,$e,$epsilon) {
|
||
|
|
$lambda = atan ($Y/$X);
|
||
|
|
|
||
|
|
$P = sqrt($X*$X + $Y*$Y);
|
||
|
|
$phi[0] = atan ($Z/ ($P * (1 - ( ($a*$e*$e)/sqrt($X*$X + $Y*$Y + $Z*$Z) ) ) ) );
|
||
|
|
|
||
|
|
$i = 0;
|
||
|
|
do
|
||
|
|
{
|
||
|
|
$i++;
|
||
|
|
$temp = pow((1 - ( $a * $e*$e * cos($phi[$i - 1] )/( $P * sqrt(1 - $e*$e * sin($phi[$i - 1])*sin($phi[$i - 1]))))),-1);
|
||
|
|
$phi[$i] = atan( $temp * $Z / $P );
|
||
|
|
}
|
||
|
|
while (abs($phi[$i] - $phi[$i - 1]) >= $epsilon);
|
||
|
|
|
||
|
|
$phix = $phi[$i];
|
||
|
|
|
||
|
|
$he = ($P/cos($phix)) - ($a/sqrt(1 - $e*$e * sin($phix)*sin($phix)));
|
||
|
|
|
||
|
|
$coords['lambda']=$lambda;
|
||
|
|
$coords['phi']=$phix;
|
||
|
|
$coords['he']=$he;
|
||
|
|
|
||
|
|
return $coords;
|
||
|
|
}
|
||
|
|
|
||
|
|
/** Transformation d'un jeu de 7 paramètres entre 2 systèmes géodésiques
|
||
|
|
** @param double $Tx Translation suivant l'axe des X (de 1 vers 2)
|
||
|
|
** @param double $Ty Translation suivant l'axe des Y (de 1 vers 2)
|
||
|
|
** @param double $Tz Translation suivant l'axe des Z (de 1 vers 2)
|
||
|
|
** @param double $D Facteur d'échelle de 1 vers 2
|
||
|
|
** @param double $Rx Angle de rotation autour de l'axe des X en radian (de 1 vers 2)
|
||
|
|
** @param double $Ry Angle de rotation autour de l'axe des Y en radian (de 1 vers 2)
|
||
|
|
** @param double $Rz Angle de rotation autour de l'axe des Z en radian (de 1 vers 2)
|
||
|
|
** @param double $U Vecteur de coordonnées cartésiennes tridimensionnelles dans le système 1 : U=(Ux,Uy,Uz)
|
||
|
|
** @return array Vecteur de coordonnées cartésiennes tridimensionnelles dans le système 2, V=(Vx,Vy,Vz)
|
||
|
|
*/
|
||
|
|
function ALG0013($Tx,$Ty,$Tz,$D,$Rx,$Ry,$Rz,$U) {
|
||
|
|
$V=array();
|
||
|
|
$V['X'] = $Tx + $U['X'] * (1 + $D) + $U['Z'] * $Ry - $U['Y'] * $Rz;
|
||
|
|
$V['Y'] = $Ty + $U['Y'] * (1 + $D) + $U['X'] * $Rz - $U['Z'] * $Rx;
|
||
|
|
$V['Z'] = $Tz + $U['Z'] * (1 + $D) + $U['Y'] * $Rx - $U['X'] * $Ry;
|
||
|
|
return $V;
|
||
|
|
}
|
||
|
|
|
||
|
|
/** Détermination des paramètres de calcul d'une projection Lambert conique
|
||
|
|
** dans le cas tangent, avec ou sans facteur d'échelle en fonction des paramètres de définition usuels
|
||
|
|
** @return array n, C, lambdac, Xs, Ys
|
||
|
|
**/
|
||
|
|
function ALG0019($lambda0,$phi0,$k0,$X0,$Y0,$a,$e) {
|
||
|
|
$lambdac = $lambda0;
|
||
|
|
$n = sin($phi0);
|
||
|
|
$C = $k0 * ALG0021($phi0,$a,$e) * tan (pi()/2 - $phi0) * exp ( $n * ALG0001($phi0,$e) );
|
||
|
|
$Xs = $X0;
|
||
|
|
$Ys = $Y0 + $k0 * ALG0021($phi0,$a,$e) * tan (pi()/2 - $phi0) ;
|
||
|
|
|
||
|
|
$tab ['e'] = $e;
|
||
|
|
$tab ['n'] = $n;
|
||
|
|
$tab ['C'] = $C;
|
||
|
|
$tab ['lambdac'] = $lambdac;
|
||
|
|
$tab ['Xs'] = $Xs;
|
||
|
|
$tab ['Ys'] = $Ys;
|
||
|
|
|
||
|
|
return $tab;
|
||
|
|
|
||
|
|
}
|
||
|
|
|
||
|
|
/** Calcul de la grande normale de l'ellipsoïde
|
||
|
|
** @param double $phi Latitude
|
||
|
|
** @param double $a Demi-grand axe de l'ellipsoïde
|
||
|
|
** @param double $e Première excentricité de l'ellipsoïde
|
||
|
|
** @return double
|
||
|
|
**/
|
||
|
|
function ALG0021($phi,$a,$e) {
|
||
|
|
$N = $a/sqrt( 1 - $e * $e * sin($phi) * sin($phi) );
|
||
|
|
return $N;
|
||
|
|
}
|
||
|
|
|
||
|
|
/** Calcul des constantes d'une projection Lambert conique conforme dans le cas sécant
|
||
|
|
** @param double $lambda0 Longitude origine en radian par rapport au méridien origine
|
||
|
|
** @param double $phi0 Latitude origine
|
||
|
|
** @param double $X0 Coordonnée X en projection du point origine
|
||
|
|
** @param double $Y0 Coordonnée Y en projection du point origine
|
||
|
|
** @param double $phi1 Latitude en radian du 1er parallèle automécoïque
|
||
|
|
** @param double $phi2 Latitude en radian du 2ème parallèle automécoïque
|
||
|
|
** @param double $a Demi-grand axe de l'ellipsoïde
|
||
|
|
** @param double $e Première excentricité de l'ellipsoïde
|
||
|
|
** @return array n, C, lambdac, Xs, Ys
|
||
|
|
**/
|
||
|
|
function ALG0054($lambda0,$phi0,$X0,$Y0,$phi1,$phi2,$a,$e) {
|
||
|
|
$lambdac = $lambda0;
|
||
|
|
$n = ( (log( (ALG0021($phi2,$a,$e)*cos($phi2))/(ALG0021($phi1,$a,$e)*cos($phi1)) ))/(ALG0001($phi1,$e) - ALG0001($phi2,$e) ));
|
||
|
|
$C = ((ALG0021($phi1,$a,$e)* cos($phi1))/$n) * exp($n * ALG0001($phi1,$e));
|
||
|
|
|
||
|
|
if ($phi0 == (pi()/2)) {
|
||
|
|
$Xs = $X0;
|
||
|
|
$Ys = $Y0;
|
||
|
|
} else {
|
||
|
|
$Xs = $X0;
|
||
|
|
$Ys = $Y0 + $C * exp(-1 * $n * ALG0001($phi0,$e));
|
||
|
|
}
|
||
|
|
$tab=array( 'e' => $e,
|
||
|
|
'n' => $n,
|
||
|
|
'C' => $C,
|
||
|
|
'lambdac' => $lambdac,
|
||
|
|
'Xs'=> $Xs,
|
||
|
|
'Ys'=> $Ys);
|
||
|
|
return $tab;
|
||
|
|
}
|
||
|
|
|
||
|
|
function Lambert2WGS84($X,$Y,$orig='L93') {
|
||
|
|
$epsilon = 0.00000000001;
|
||
|
|
|
||
|
|
$lambdac = 0.04079234433; // pour greenwich
|
||
|
|
$e = 0.08248325676; // première excentricité de l ellipsoïde Clarke 1880 français
|
||
|
|
$he = 100;
|
||
|
|
$a = 6378249.2; // demi-grand axe de l'ellipsoide
|
||
|
|
|
||
|
|
$Tx = -168;
|
||
|
|
$Ty = -60;
|
||
|
|
$Tz = +320;
|
||
|
|
$D = 0;
|
||
|
|
$Rx = $Ry = $Rz = 0;
|
||
|
|
|
||
|
|
$orig=strtoupper($orig);
|
||
|
|
switch ($orig) {
|
||
|
|
case 'LI':
|
||
|
|
case 'L1': $n = 0.7604059656;
|
||
|
|
$c = 11603796.98;
|
||
|
|
$Xs = 600000;
|
||
|
|
$Ys = 5657616.674;
|
||
|
|
break;
|
||
|
|
case 'LII':
|
||
|
|
case 'LIIE':
|
||
|
|
case 'L2E':
|
||
|
|
case 'L2': $n = 0.7289686274;
|
||
|
|
$c = 11745793.39;
|
||
|
|
$Xs = 600000;
|
||
|
|
if ($orig=='L2E' || $orig=='LIIE')
|
||
|
|
$Ys = 8199695.768;
|
||
|
|
else $Ys = 6199695.768;
|
||
|
|
break;
|
||
|
|
case 'LIII':
|
||
|
|
case 'L3': $n = 0.6959127966;
|
||
|
|
$c = 11947992.52;
|
||
|
|
$Xs = 600000;
|
||
|
|
$Ys = 6791905.085;
|
||
|
|
break;
|
||
|
|
case 'LIV':
|
||
|
|
case 'L4': $n = 0.6712679322;
|
||
|
|
$c = 12136281.99;
|
||
|
|
$Xs = 234.358;
|
||
|
|
$Ys = 7239161.542;
|
||
|
|
break;
|
||
|
|
case 'L93':
|
||
|
|
default: $n = 0.7256077650;
|
||
|
|
$c = 11745255.426;
|
||
|
|
$Xs = 700000;
|
||
|
|
$Ys = 12655612.050;
|
||
|
|
break;
|
||
|
|
}
|
||
|
|
|
||
|
|
$coords = ALG0004($X,$Y,$n,$c,$Xs,$Ys,$lambdac,$e,$epsilon);
|
||
|
|
|
||
|
|
$coords = ALG0009($coords['lambda'],$coords['phi'],$he,$a,$e);
|
||
|
|
|
||
|
|
$coords = ALG0013($Tx,$Ty,$Tz,$D,$Rx,$Ry,$Rz,$coords);
|
||
|
|
|
||
|
|
$a = 6378137.0; // ellipsoïdes WGS84
|
||
|
|
$f = 1/298.257223563;
|
||
|
|
$b = $a*(1-$f);
|
||
|
|
$e = sqrt(($a*$a - $b*$b)/($a*$a));
|
||
|
|
|
||
|
|
$X = $coords['X'];
|
||
|
|
$Y = $coords['Y'];
|
||
|
|
$Z = $coords['Z'];
|
||
|
|
$coords = ALG0012($X,$Y,$Z,$a,$e,$epsilon);
|
||
|
|
|
||
|
|
$xy['long'] = rad2deg($coords['lambda']);
|
||
|
|
$xy['lat'] = rad2deg($coords['phi']);
|
||
|
|
return $xy;
|
||
|
|
}
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Function to convert geographic coordinates to "lambert 2 etendue" coordinates
|
||
|
|
* from: http://www.forumsig.org/showthread.php?p=64050#post64050
|
||
|
|
**/
|
||
|
|
function geos2lambert($latitude,$longitude) {
|
||
|
|
//0)degres-minutes-secondes + orientation (d,m,s,o) en radian
|
||
|
|
$lambda_w = $longitude * pi()/180 ;
|
||
|
|
$phi_w = $latitude * pi()/180 ;
|
||
|
|
|
||
|
|
//1) coordonnées géographiques WGS84 (phi_w,lambda_w) en coordonnées cartésiennes WGS84 (X_w,Y_w,Z_w)
|
||
|
|
$a_w = 6378137.0;
|
||
|
|
$b_w = 6356752.314;
|
||
|
|
|
||
|
|
$e2_w = ($a_w*$a_w-$b_w*$b_w)/($a_w*$a_w);
|
||
|
|
|
||
|
|
//et on a la grande normale de l'ellipsoide WGS84
|
||
|
|
$N = $a_w/sqrt(1-$e2_w*pow(sin($phi_w),2));
|
||
|
|
|
||
|
|
//ainsi on obtient
|
||
|
|
$X_w = $N * cos($phi_w) * cos($lambda_w);
|
||
|
|
$Y_w = $N * cos($phi_w) * sin($lambda_w);
|
||
|
|
$Z_w = $N * (1-$e2_w) * sin($phi_w);
|
||
|
|
//Ref.: http://www.ign.fr/telechargement/MPr...RCE/NTG_80.pdf and http://de.wikipedia.org/wiki/WGS84
|
||
|
|
|
||
|
|
//2) coordonnées cartésiennes WGS84 (X_w,Y_w,Z_w) en coordonnées cartésiennes NTF (X_n,Y_n,Z_n)
|
||
|
|
|
||
|
|
$dX = 168.0;
|
||
|
|
$dY = 60.0;
|
||
|
|
$dZ = -320.0;
|
||
|
|
|
||
|
|
$X_n = $X_w + $dX;
|
||
|
|
$Y_n = $Y_w + $dY;
|
||
|
|
$Z_n = $Z_w + $dZ;
|
||
|
|
//ref.: http://support.esrifrance.fr/Documents/Generalites/Projections/Generalites/Generalites.htm#2
|
||
|
|
|
||
|
|
//3) coordonnées cartésiennes NTF (X_n,Y_n,Z_n) en coordonnées géographiques NTF (phi_n,lambda_n)
|
||
|
|
$a_n = 6378249.2;
|
||
|
|
$b_n = 6356515.0;
|
||
|
|
|
||
|
|
$e2_n = ($a_n*$a_n-$b_n*$b_n)/($a_n*$a_n);
|
||
|
|
//on définit une tolérance de convergence
|
||
|
|
$epsilon = pow(10,-10);
|
||
|
|
|
||
|
|
//puis on amorce une boucle de calcul
|
||
|
|
$p0=atan($Z_n/sqrt($X_n*$X_n+$Y_n*$Y_n)*(1-($a_n*$e2_n)/(sqrt($X_n*$X_n+$Y_n*$Y_n+$Z_n*$Z_n))));
|
||
|
|
$p1=atan(($Z_n/sqrt($X_n*$X_n+$Y_n*$Y_n))/(1-($a_n*$e2_n*cos($p0))/(sqrt(($X_n*$X_n+$Y_n*$Y_n)*(1-$e2_n*pow(sin($p0),2))))));
|
||
|
|
while(!(abs($p1-$p0)<$epsilon)){
|
||
|
|
$p0 = $p1;
|
||
|
|
$p1 = atan(($Z_n/sqrt($X_n*$X_n+$Y_n*$Y_n))/(1-($a_n*$e2_n*cos($p0))/(sqrt(($X_n*$X_n+$Y_n*$Y_n)*(1-$e2_n*pow(sin($p0),2))))));
|
||
|
|
}
|
||
|
|
$phi_n = $p1;
|
||
|
|
$lambda_n = atan($Y_n/$X_n);
|
||
|
|
|
||
|
|
//4) coordonnées géographiques NTF (phi_n,lambda_n) en coordonnées projetées en Lambert II étendu (X_l2e, Y_l2e)
|
||
|
|
$n = 0.7289686274;
|
||
|
|
$c = 11745793.39;
|
||
|
|
$Xs = 600000.0;
|
||
|
|
$Ys = 8199695.768;
|
||
|
|
|
||
|
|
$e_n = sqrt($e2_n);
|
||
|
|
$lambda0 = 0.04079234433198; //correspond à la longitude en radian de Paris (2°20'14.025" E) par rapport à Greenwich
|
||
|
|
//puis on calcule la latitude isométrique
|
||
|
|
$L = log(tan(pi()/4 + $phi_n/2) * pow(((1-$e_n*sin($phi_n))/(1+$e_n*sin($phi_n))),($e_n/2)));
|
||
|
|
|
||
|
|
//enfin on projette
|
||
|
|
$X_l2e = $Xs + $c*exp((-$n*$L))*sin($n*($lambda_n-$lambda0));
|
||
|
|
$Y_l2e = $Ys - $c*exp((-$n*$L))*cos($n*($lambda_n-$lambda0));
|
||
|
|
|
||
|
|
return array("x_l2e"=>$X_l2e,"y_l2e"=>$Y_l2e);
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
/** Conversion WGS84 en Lambert93
|
||
|
|
** NTG_71.pdf
|
||
|
|
** http://www.ign.fr/affiche_rubrique.asp?rbr_id=1700&lng_id=FR
|
||
|
|
**/
|
||
|
|
function geos2lambert93($latitude,$longitude) {
|
||
|
|
|
||
|
|
// Système WGS84
|
||
|
|
$a=6378137; // demi grand axe de l'ellipsoide (m)
|
||
|
|
$e=0.08181919106; // première excentricitè de l'ellipsoide
|
||
|
|
|
||
|
|
// Paramètres de projections
|
||
|
|
$l0=$lc=deg2rad(3); // longitude de référence
|
||
|
|
$phi0=deg2rad(46.5); // latitude d'origine en radian
|
||
|
|
$phi1=deg2rad(44); // 1er parallele automécoïque
|
||
|
|
$phi2=deg2rad(49); // 2eme parallele automécoïque
|
||
|
|
|
||
|
|
$x0=700000; //coordonnées à l'origine
|
||
|
|
$y0=6600000; //coordonnées à l'origine
|
||
|
|
|
||
|
|
// coordonnées du point à traduire
|
||
|
|
$phi=deg2rad($latitude);
|
||
|
|
$l=deg2rad($longitude);
|
||
|
|
|
||
|
|
// calcul des grandes normales
|
||
|
|
$gN1=$a/sqrt(1-$e*$e*sin($phi1)*sin($phi1));
|
||
|
|
$gN2=$a/sqrt(1-$e*$e*sin($phi2)*sin($phi2));
|
||
|
|
|
||
|
|
//calculs de slatitudes isométriques
|
||
|
|
$gl1=log(tan(pi()/4+$phi1/2)*pow((1-$e*sin($phi1))/(1+$e*sin($phi1)),$e/2));
|
||
|
|
$gl2=log(tan(pi()/4+$phi2/2)*pow((1-$e*sin($phi2))/(1+$e*sin($phi2)),$e/2));
|
||
|
|
$gl0=log(tan(pi()/4+$phi0/2)*pow((1-$e*sin($phi0))/(1+$e*sin($phi0)),$e/2));
|
||
|
|
$gl=log(tan(pi()/4+$phi/2)*pow((1-$e*sin($phi))/(1+$e*sin($phi)),$e/2));
|
||
|
|
|
||
|
|
//calcul de l'exposant de la projection
|
||
|
|
$n=(log(($gN2*cos($phi2))/($gN1*cos($phi1))))/($gl1-$gl2);//ok
|
||
|
|
|
||
|
|
//calcul de la constante de projection
|
||
|
|
$c=(($gN1*cos($phi1))/$n)*exp($n*$gl1);//ok
|
||
|
|
|
||
|
|
//calcul des coordonnées
|
||
|
|
$ys=$y0+$c*exp(-1*$n*$gl0);
|
||
|
|
|
||
|
|
//calcul des coordonnées lambert
|
||
|
|
$x93=$x0+$c*exp(-1*$n*$gl)*sin($n*($l-$lc));
|
||
|
|
$y93=$ys-$c*exp(-1*$n*$gl)*cos($n*($l-$lc));
|
||
|
|
|
||
|
|
return array("x_93"=>$x93,"y_93"=>$y93);
|
||
|
|
}
|