Skip to content

Instantly share code, notes, and snippets.

@BackFront
Created January 16, 2020 19:27
Show Gist options
  • Save BackFront/ec5cc716faba0dbbe4b8cf2e933e17f8 to your computer and use it in GitHub Desktop.
Save BackFront/ec5cc716faba0dbbe4b8cf2e933e17f8 to your computer and use it in GitHub Desktop.
A seguir estão as funções escritas em PHP que realizam conversão de coordenadas de projeção para coordenadas geográficas e vice versa.
<?php
/********************************************************************
A seguir estão as funções escritas em PHP que realizam conversão de
coordenadas de projeção para coordenadas geográficas e vice versa e
que são utilizadas na Calculadora Geografica do INPE
(www.dpi.inpe.br/calcula/)
Referencia:
SNYDER, John Parr. Map projections--A working manual.
US Government Printing Office, 1987.
********************************************************************/
/********************************************************************
Parametros dos datuns utilizados
********************************************************************/
$datum_parametros = array
(
//id datum semi_eixo achat deltax deltay deltaz
"0 Null 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00",
"1 SAD69 6.3781600e+06 3.35289187e-03 -6.735000e+01 3.880000e+00 -3.822000e+01",
"2 CorregoAlegre 6.3783880e+06 3.36700337e-03 -2.060500e+02 1.682800e+02 -3.820000e+00",
"3 AstroChua 6.3783880e+06 3.36700337e-03 -1.443500e+02 2.433700e+02 -3.322000e+01",
"4 WGS84 6.3781370e+06 3.35281066e-03 0.000000e+00 0.000000e+00 0.000000e+00",
"5 SIRGAS2000 6.3781370e+06 3.35281068e-03 0.000000e+00 0.000000e+00 0.000000e+00",
);
/********************************************************************
geo_to_utm
Autor : Julio Cesar Lima d'Alge jul-88
Conversao p/ Php : Luis Maurano 2008
Resumo : transforma coordenadas geodesicas em coordenadas UTM
Entradas :
lat coordenadas geodesicas (em radianos)
lon coordenadas geodesicas (em radianos)
lon_mc meridiano central (em radianos)
semi_eixo semi_eixo_maior do elipsoide
achat achatamento do elipsoide
hemis norte ou sul
Saidas : x y coordenadas UTM (em metros)
********************************************************************/
function geo_to_utm($lat,$lon,$semi_eixo,$achat,$hemis,$lon_mc,&$x,&$y)
{
// define PI
$pi = pi();
// verifica hemisferio
if ($hemis == "norte")
$offy = 0.;
else
$offy = 10000000;
// Converte Lat Long em UTM
$k0 = 1. - (1./2500.);
$equad = 2.*$achat - pow($achat,(double)2);
$elinquad = $equad/(1. - $equad);
$aux1 = $equad*$equad;
$aux2 = $aux1*$equad;
$aux3 = sin((double)2*$lat);
$aux4 = sin((double)4*$lat);
$aux5 = sin((double)6*$lat);
$aux6 = (1. - $equad/4. - 3.*$aux1/64. - 5.*$aux2/256.)*$lat;
$aux7 = (3.*$equad/8. + 3.*$aux1/32. + 45.*$aux2/1024.)*$aux3;
$aux8 = (15.*$aux1/256. + 45.*$aux2/1024.)*$aux4;
$aux9 = (35.*$aux2/3072.)*$aux5;
$n = $semi_eixo/sqrt((double)1-$equad*pow(sin($lat),(double)2));
$t = pow(tan($lat),(double)2);
$c = $elinquad*pow(cos($lat),(double)2);
$ag = ($lon-$lon_mc)*cos($lat);
$m = $semi_eixo*($aux6 - $aux7 + $aux8 - $aux9);
$aux10 = (1.-$t+$c)*$ag*$ag*$ag/6.;
$aux11 = (5.-18.*$t+$t*$t+72.*$c-58.*$elinquad)*(pow($ag,(double)5))/120.;
$aux12 = (5.-$t+9.*$c+4.*$c*$c)*$ag*$ag*$ag*$ag/24.;
$aux13 = (61.-58.*$t+$t*$t+600.*$c-330.*$elinquad)*(pow($ag,(double)6))/720.;
$x = 500000. + $k0*$n*($ag + $aux10 + $aux11);
$y = $offy + $k0*($m + $n*tan($lat)*($ag*$ag/2. + $aux12 + $aux13));
}
/********************************************************************
utm_to_geo
Autor : Julio Cesar Lima d'Alge jul-88
Conversao p/ Php : Luis Maurano 2008
Resumo : transforma coordenadas UTM em coordenadas geodesicas
Entradas :
x coordenadas UTM (em metros)
y coordenadas UTM (em metros)
lon_mc meridiano central (em radianos)
semi_eixo semi_eixo maior do elipsoide
achat achatamento do elipsoide
hemis norte ou sul
Saidas : lat lon coordenadas geodesicas (em radianos)
********************************************************************/
function utm_to_geo($x,$y,$semi_eixo,$achat,$hemis,$lon_mc,&$lat,&$lon)
{
// define PI
$pi = pi();
// verifica hemisferio
if ($hemis == "norte")
$y = $y + 10000000;
// Converte UTM para Lat/Long;
$k0 = 1. - (1./2500.);
$equad = 2.*$achat - pow($achat,(double)2);
$elinquad = $equad/(1. - $equad);
$e1 = (1.-sqrt((double)1-$equad))/(1.+sqrt((double)1-$equad));
$aux1 = $equad*$equad;
$aux2 = $aux1*$equad;
$aux3 = $e1*$e1;
$aux4 = $e1*$aux3;
$aux5 = $aux4*$e1;
$m = ($y - 10000000.)/$k0;
$mi = $m/($semi_eixo*(1.-$equad/4.-3.*$aux1/64.-5.*$aux2/256.));
$aux6 = (3.*$e1/2. - 27.*$aux4/32.)*sin((double)2*$mi);
$aux7 = (21.*$aux3/16. - 55.*$aux5/32.)*sin((double)4*$mi);
$aux8 = (151.*$aux4/96.)*sin((double)6*$mi);
$lat1 = $mi + $aux6 + $aux7 + $aux8;
$c1 = $elinquad*pow(cos($lat1),(double)2);
$t1 = pow(tan($lat1),(double)2);
$n1 = $semi_eixo/sqrt((double)1-$equad*pow(sin($lat1),(double)2));
$quoc = pow(((double)1-$equad*sin($lat1)*sin($lat1)),(double)3);
$r1 = $semi_eixo*(1.-$equad)/sqrt($quoc);
$d = ($x - 500000.)/($n1*$k0);
$aux9 = (5.+3.*$t1+10.*$c1-4.*$c1*$c1-9.*$elinquad)*$d*$d*$d*$d/24.;
$aux10 = (61.+90.*$t1+298.*$c1+45.*$t1*$t1-252.*$elinquad-3.*$c1*$c1)*pow($d,(double)6)/720.;
$aux11 = $d - (1.+ 2.*$t1 + $c1)*$d*$d*$d/6.;
$aux12 = (5.-2.*$c1+28.*$t1-3.*$c1*$c1+8.*$elinquad+24.*$t1*$t1)*pow($d,(double)5)/120.;
$lat = $lat1 - ($n1*tan($lat1)/$r1)*($d*$d/2. - $aux9 + $aux10);
$lon = $lon_mc + ($aux11 + $aux12)/cos($lat1);
}
/********************************************************************
poli_to_geo
Autor : Julio Cesar Lima d'Alge fev-90
Conversao p/ Php : Luis Maurano 2008
Resumo : transforma coordenadas policonicas em coordenadas geodesicas
Entradas :
x coordenadas policonicas (em metros)
y coordenadas policonicas (em metros)
lat0 paralelo origem (em radianos)
long0 meridiano origem (em radianos)
semi_eixo semi_eixo maior do elipsoide
achat achatamento do elipsoide
Saidas : lat lon coordenadas geodesicas (em radianos)
********************************************************************/
function poli_to_geo($x,$y,$semi_eixo,$achat,$latO,$longO,&$lat,&$lon)
{
$equad = 2.*$achat - pow($achat,(double)2);
/* Para calculo de mo lat vira lat0 */
$aux01 = (1.-$equad/4.-3.*$equad*$$equad/64.-5.*$equad*$equad*$equad/256.)*$latO;
$aux02 = (3.*$equad/8+3.*$equad*$equad/32.+45.*$equad*$equad*$equad/1024.)*sin((double)2*$latO);
$aux03 = (15.*$equad*$equad/256.+45.*$equad*$equad*$equad/1024.)*sin((double)4*$latO);
$aux04 = (35.*$equad*$equad*$equad/3072.)*sin((double)6*$latO);
$m0 = $semi_eixo*($aux01 - $aux02 + $aux03 - $aux04);
if ($y == (-$m0))
{
$lat = 0.;
$lon = $x/$semi_eixo + $longO;
}
else
{
$A = ($m0 + $y)/$semi_eixo;
$B = (($x*$x)/($semi_eixo*$semi_eixo)) +($A*$A);
$lat2 = $A; /* Inicializando a $latitude para a iteracao */
do
{
$C = (sqrt(1.- $equad*sin($lat2)*sin($lat2)))*tan($lat2);
/* Calculo de mn */
$aux21 = (1.-$equad/4.-3.*$equad*$equad/64.-5.*$equad*$equad*$equad/256.)*$lat2;
$aux22 = (3.*$equad/8.+3.*$equad*$equad/32.+45.*$equad*$equad*$equad/1024.)*sin((double)2*$lat2);
$aux23 = (15.*$equad*$equad/256.+45.*$equad*$equad*$equad/1024.)*sin((double)4*$lat2);
$aux24 = (35.*$equad*$equad*$equad/3072.)*sin((double)6*$lat2);
$mn = $semi_eixo*($aux21 - $aux22 + $aux23 - $aux24);
/* Calculo de mnl*/
$aux05 = 1.- $equad/4.-3.*$equad*$equad/64.-5.*$equad*$equad*$equad/256.;
$aux06 = 2.*(3.*$equad/8.+3.*$equad*$equad/32.+45.*$equad*$equad*$equad/1024.)*cos((double)2*$lat2);
$aux07 = 4.*(15.*$equad*$equad/256.+45.*$equad*$equad*$equad/1024.)*cos((double)4*$lat2);
$aux08 = 6.*(35.*$equad*$equad*$equad/3072.)*cos((double)6*$lat2);
$mnl = $aux05 - $aux06 + $aux07- $aux08;
/* Calculo de ma */
$ma = $mn/$semi_eixo;
$aux09 = ($A*($C*$ma+1)-$ma)-(0.5*($ma*$ma+$B)*$C);
$aux10 = $equad*sin((double)2*$lat2)*($ma*$ma+$B-2.*$A*$ma);
$aux11 = 4.*$C+($A-$ma);
$aux12 = $C*$mnl-(2./sin((double)2*$lat2));
/* Calculo da nova $latitude */
$lat1 = $lat2 - ($aux09/(($aux10/($aux11*$aux12)) - $mnl));
$cp = abs($lat1-$lat2) ;
$lat2 = $lat1;
}
while($cp > 0.000000001);
$lat = $lat1;
$lon = ((asin(($x*$C)/$semi_eixo))/(sin($lat1))) + $longO;
}
}
/********************************************************************
geo_to_poli
Autor : Julio Cesar Lima d'Alge fev-90
Conversao p/ Php : Luis Maurano 2008
Resumo : transforma coordenadas geodesicas em coordenadas policonicas
Entradas :
lat coordenadas geodesicas (em radianos)
lon coordenadas geodesicas (em radianos)
lat0 paralelo origem (em radianos)
long0 meridiano origem (em radianos)
semi_eixo semi_eixo maior do elipsoide
achat achatamento do elipsoide
Saidas : x y coordenadas policonicas (em metros)
********************************************************************/
function geo_to_poli($lat,$lon,$semi_eixo,$achat,$latO,$longO,&$x,&$y)
{
$equad = 2.*$achat - pow($achat,(double)2);
$aux01 = (1.-$equad/4.-3.*$equad*$equad/64.-5.*$equad*$equad*$equad/256.)*$latO;
$aux02 = (3.*$equad/8+3.*$equad*$equad/32.+45.*$equad*$equad*$equad/1024.)*sin((double)2*$latO);
$aux03 = (15.*$equad*$equad/256.+45.*$equad*$equad*$equad/1024.)*sin((double)4*$latO);
$aux04 = (35.*$equad*$equad*$equad/3072.)*sin((double)6*$latO);
$m0 = $semi_eixo*($aux01 - $aux02 + $aux03 - $aux04);
if ($lat == 0.)
{
$x = $semi_eixo*($lon - $longO);
$y = -$m0;
}
else
{
$aux1 = (1.-$equad/4.-3.*$equad*$equad/64.-5.*$equad*$equad*$equad/256.)*$lat;
$aux2 = (3.*$equad/8+3.*$equad*$equad/32.+45.*$equad*$equad*$equad/1024.)*sin((double)2*$lat);
$aux3 = (15.*$equad*$equad/256.+45.*$equad*$equad*$equad/1024.)*sin((double)4*$lat);
$aux4 = (35.*$equad*$equad*$equad/3072.)*sin((double)6*$lat);
$m = $semi_eixo*($aux1 - $aux2 + $aux3 - $aux4);
$n = $semi_eixo/sqrt((double)1 - $equad*pow(sin($lat),(double)2));
$e = ($lon - $longO)*sin($lat);
$x = $n*sin($e)/tan($lat);
$y = $m - $m0 + ($n*(1. - cos($e))/tan($lat));
}
}
/********************************************************************
lamb_to_geo
Autor : Julio Cesar Lima d'Alge set-88
Conversao p/ Php : Luis Maurano 2008
Resumo : transforma coordenadas Lambert em coordenadas geodesicas
Entradas :
x coordenadas Lambert (em metros)
y coordenadas Lambert (em metros)
lat0 paralelo origem (em radianos)
long0 meridiano origem (em radianos)
lat1 paralelo padrao 1 (em radianos)
lat2 paralelo padrao 2 (em radianos)
semi_eixo semi_eixo maior do elipsoide
achat achatamento do elipsoide
Saidas : lat lon coordenadas geodesicas (em radianos)
********************************************************************/
function lamb_to_geo($x,$y,$latO,$longO,$lat1,$lat2,$semi_eixo,$achat,&$lat,&$lon)
{
// define PI
$pi = pi();
//$pi = 4.*atan((double)1);
$equad = 2.*$achat - pow($achat,(double)2);
$e = sqrt($equad);
$m1 = cos($lat1)/sqrt((double)1-$equad*pow(sin($lat1),(double)2));
$m2 = cos($lat2)/sqrt((double)1-$equad*pow(sin($lat2),(double)2));
$aux1 = sqrt(((double)1-$e*sin($lat1))/((double)1+$e*sin($lat1)));
$aux2 = sqrt(((double)1-$e*sin($lat2))/((double)1+$e*sin($lat2)));
$aux0 = sqrt(((double)1-$e*sin($latO))/((double)1+$e*sin($latO)));
$t1 = ((1.-tan($lat1/(double)2))/(1.+tan($lat1/(double)2)))/pow($aux1,$e);
$t2 = ((1.-tan($lat2/(double)2))/(1.+tan($lat2/(double)2)))/pow($aux2,$e);
$t0 = ((1.-tan($latO/(double)2))/(1.+tan($latO/(double)2)))/pow($aux0,$e);
if ($lat1 == $lat2)
$n = sin($lat1);
else
$n = (log($m1)-log($m2))/(log($t1)-log($t2));
$efe = $m1/($n*pow($t1,$n));
$ro0 = $semi_eixo*$efe*pow($t0,$n);
$sinal = (int)($n/abs($n));
$ro = sqrt($x*$x + ($ro0-$y)*($ro0-$y));
$ro *= $sinal;
$teta = atan($x/($ro0-$y));
$t = pow(($ro/($semi_eixo*$efe)),(double)1/$n);
$xx = $pi/2. - 2.*atan($t);
$aux3 = $equad/2. + 5.*$equad*$equad/24. + $equad*$equad*$equad/12.;
$aux4 = 7.*$equad*$equad/48. + 29.*$equad*$equad*$equad/240.;
$aux5 = (7.*$equad*$equad*$equad/120.)*sin(12.*atan($t));
$lat = $xx + $aux3*sin(4.*atan($t)) - $aux4*sin(8.*atan($t)) + $aux5;
$lon = $teta/$n + $longO;
}
/********************************************************************
geo_to_lamb
Autor : Julio Cesar Lima d'Alge set-88
Conversao p/ Php : Luis Maurano 2008
Resumo : transforma coordenadas geodesicas em coordenadas Lambert
Entradas :
lat coordenadas geodesicas (em radianos)
lon coordenadas geodesicas (em radianos)
lat0 paralelo origem (em radianos)
long0 meridiano origem (em radianos)
lat1 paralelo padrao 1 (em radianos)
lat2 paralelo padrao 2 (em radianos)
semi_eixo semi_eixo maior do elipsoide
achat achatamento do elipsoide
Saidas : x y coordenadas Lambert (em metros)
********************************************************************/
function geo_to_lamb($lat,$lon,$latO,$longO,$lat1,$lat2,$semi_eixo,$achat,&$x,&$y)
{
$equad = 2.*$achat - pow($achat,(double)2);
$e = sqrt($equad);
$m1 = cos($lat1)/sqrt((double)1-$equad*pow(sin($lat1),(double)2));
$m2 = cos($lat2)/sqrt((double)1-$equad*pow(sin($lat2),(double)2));
$aux1 = sqrt(((double)1-$e*sin($lat1))/((double)1+$e*sin($lat1)));
$aux2 = sqrt(((double)1-$e*sin($lat2))/((double)1+$e*sin($lat2)));
$aux0 = sqrt(((double)1-$e*sin($latO))/((double)1+$e*sin($latO)));
$t1 = ((1.-tan($lat1/(double)2))/(1.+tan($lat1/(double)2)))/pow($aux1,$e);
$t2 = ((1.-tan($lat2/(double)2))/(1.+tan($lat2/(double)2)))/pow($aux2,$e);
$t0 = ((1.-tan($latO/(double)2))/(1.+tan($latO/(double)2)))/pow($aux0,$e);
if ($lat1 == $lat2)
$n = sin($lat1);
else
$n = (log($m1)-log($m2))/(log($t1)-log($t2));
$efe = $m1/($n*pow($t1,$n));
$ro0 = $semi_eixo*$efe*pow($t0,$n);
$aux = sqrt(((double)1-$e*sin($lat))/((double)1+$e*sin($lat)));
$t = ((1.-tan($lat/(double)2))/(1.+tan($lat/(double)2)))/pow($aux,$e);
$ro = $semi_eixo*$efe*pow($t,$n);
$teta = $n*($lon-$longO);
$x = $ro*sin($teta);
$y = $ro0 - $ro*cos($teta);
}
/********************************************************************
merc_to_geo
Autor : Julio Cesar Lima d'Alge jul-88
Conversao p/ Php : Luis Maurano 2008
Resumo : transforma coordenadas Mercator em coordenadas geodesicas
Entradas :
x coordenadas Mercator (em metros)
y coordenadas Mercator (em metros)
lat1 latitude padrao ou reduzida (em radianos)
long0 meridiano origem (em radianos)
semi_eixo semi_eixo maior do elipsoide
achat achatamento
Saidas : lat lon coordenadas geodesicas (em radianos)
********************************************************************/
function merc_to_geo($x,$y,$lat1,$longO,$semi_eixo,$achat,&$lat,&$lon)
{
// define PI
$pi = pi();
$equad = 2.*$achat - pow($achat,(double)2);
$aux1 = cos($lat1);
$aux2 = 1./sqrt((double)1-$equad*pow(sin($lat1),(double)2));
$x = $x/($aux1*$aux2);
$y = $y/($aux1*$aux2);
$t = exp(-$y/$semi_eixo);
$xx = $pi/2. - 2.*atan($t);
$aux3 = ($equad/2.+5.*$equad*$equad/24.+$equad*$equad*$equad/12.)*sin(4.*atan($t));
$aux4 = -(7.*$equad*$equad/48.+29.*$equad*$equad*$equad/240.)*sin(8.*atan($t));
$aux5 = (7.*$equad*$equad*$equad/120.)*sin(12.*atan($t));
$lat = $xx + $aux3 + $aux4 + $aux5;
$lon = $x/$semi_eixo + $longO;
}
/********************************************************************
geo_to_merc
Autor : Julio Cesar Lima d'Alge jul-88
Conversao p/ Php : Luis Maurano 2008
Resumo : transforma coordenadas geodesicas em coordenadas Mercator
Entradas :
lat coordenadas geodesicas (em radianos)
lon coordenadas geodesicas (em radianos)
lat1 latitude padrao ou reduzida (em radianos)
long0 meridiano origem (em radianos)
semi_eixo semi_eixo maior do elipsoide
achat achatamento do elipsoide
Saidas : x y coordenadas Mercator (em metros)
********************************************************************/
function geo_to_merc($lat,$lon,$lat1,$longO,$semi_eixo,$achat,&$x,&$y)
{
$equad = 2.*$achat - pow($achat,(double)2);
$aux1 = (1. + tan($lat/(double)2))/(1. - tan($lat/(double)2));
$aux2 = ($equad+$equad*$equad/4.+$equad*$equad*$equad/8.)*sin($lat);
$aux3 = ($equad*$equad/12.+$equad*$equad*$equad/16.)*sin((double)3*$lat);
$aux4 = ($equad*$equad*$equad/80.)*sin((double)5*$lat);
$aux5 = cos($lat1);
$aux6 = 1./sqrt((double)1-$equad*pow(sin($lat1),(double)2));
$x = $semi_eixo*($lon - $longO)*$aux5*$aux6;
$y = $semi_eixo*(log($aux1) - $aux2 + $aux3 - $aux4)*$aux5*$aux6;
}
/********************************************************************
alb_to_geo
Autor : Julio Cesar Lima d'Alge mar-90
Conversao p/ Php : Luis Maurano 2008
Resumo : transforma coordenadas Albers em coordenadas geodesicas
Entradas :
x coordenadas Albers (em metros)
y coordenadas Albers (em metros)
lat0 paralelo origem (em radianos)
long0 meridiano origem (em radianos)
lat1 paralelo padrao 1 (em radianos)
lat2 paralelo padrao 2 (em radianos)
semi_eixo semi_eixo_maior do elipsoide
achat achatamento do elipsoide
Saidas : lat lon coordenadas geodesicas (em radianos)
********************************************************************/
function alb_to_geo($x,$y,$latO,$longO,$lat1,$lat2,$semi_eixo,$achat,&$lat,&$lon)
{
$sinal = (int)($lat2/abs($lat2));
$equad = 2.*$achat - pow($achat,(double)2);
$e = sqrt($equad);
$m1 = cos($lat1)/sqrt((double)1-$equad*pow(sin($lat1),(double)2));
$m2 = cos($lat2)/sqrt((double)1-$equad*pow(sin($lat2),(double)2));
$aux10 = sin($latO)/((double)1-$equad*pow(sin($latO),(double)2));
$aux11 = sin($lat1)/((double)1-$equad*pow(sin($lat1),(double)2));
$aux12 = sin($lat2)/((double)1-$equad*pow(sin($lat2),(double)2));
$aux20 = log((1. - $e*sin($latO))/(1. + $e*sin($latO)));
$aux21 = log((1. - $e*sin($lat1))/(1. + $e*sin($lat1)));
$aux22 = log((1. - $e*sin($lat2))/(1. + $e*sin($lat2)));
$q0 = (1. - $equad)*($aux10 - (1./(2.*$e))*$aux20);
$q1 = (1. - $equad)*($aux11 - (1./(2.*$e))*$aux21);
$q2 = (1. - $equad)*($aux12 - (1./(2.*$e))*$aux22);
$n = ($m1*$m1 - $m2*$m2)/($q2 - $q1);
$c = $m1*$m1 + $n*$q1;
$ro0 = $semi_eixo*sqrt($c - $n*$q0)/$n;
$ro = sqrt($x*$x + ($ro0 - $y)*($ro0 - $y));
$q = ($c - ($ro*$ro*$n*$n/($semi_eixo*$semi_eixo)))/$n;
$aux = ((1. - $equad)/(2.*$e))*log((1. - $e)/(1. + $e));
$beta = asin($q/(1. - $aux));
$aux1 = ($equad/3. + 31.*$equad*$equad/180. + 517.*$equad*$equad*$equad/5040.)*sin(2.*$beta);
$aux2 = (23.*$equad*$equad/360. + 251.*$equad*$equad*$equad/3780.)*sin(4.*$beta);
$aux3 = (761.*$equad*$equad*$equad/45360.)*sin(6.*$beta);
$teta = abs(atan($x/($ro0 - $y)));
if ($sinal == 1)
{
if ($x < 0.)
$teta = -$teta;
}
if ($sinal == -1)
{
if ($x > 0.)
$teta *= $sinal;
}
$lat = $beta + $aux1 + $aux2 + $aux3;
$lon = $longO + ($teta/$n);
}
/********************************************************************
geo_to_alb
Autor : Julio Cesar Lima d'Alge mar-90
Conversao p/ Php : Luis Maurano 2008
Resumo : transforma coordenadas geodesicas em coordenadas Albers
Entradas :
lat coordenadas geodesicas (em radianos)
lon coordenadas geodesicas (em radianos)
lat0 paralelo origem (em radianos)
long0 meridiano origem (em radianos)
lat1 paralelo padrao 1 (em radianos)
lat2 paralelo padrao 2 (em radianos)
semi_eixo semi_eixo maior do elipsoide
achat achatamento do elipsoide
Saidas : x y coordenadas Albers (em metros)
********************************************************************/
function geo_to_alb($lat,$lon,$latO,$longO,$lat1,$lat2,$semi_eixo,$achat,&$x,&$y)
{
$equad = 2.*$achat - pow($achat,(double)2);
$e = sqrt($equad);
$m1 = cos($lat1)/sqrt((double)1-$equad*pow(sin($lat1),(double)2));
$m2 = cos($lat2)/sqrt((double)1-$equad*pow(sin($lat2),(double)2));
$aux1 = sin($lat)/((double)1-$equad*pow(sin($lat),(double)2));
$aux10 = sin($latO)/((double)1-$equad*pow(sin($latO),(double)2));
$aux11 = sin($lat1)/((double)1-$equad*pow(sin($lat1),(double)2));
$aux12 = sin($lat2)/((double)1-$equad*pow(sin($lat2),(double)2));
$aux2 = log((1. - $e*sin($lat))/(1. + $e*sin($lat)));
$aux20 = log((1. - $e*sin($latO))/(1. + $e*sin($latO)));
$aux21 = log((1. - $e*sin($lat1))/(1. + $e*sin($lat1)));
$aux22 = log((1. - $e*sin($lat2))/(1. + $e*sin($lat2)));
$q0 = (1. - $equad)*($aux10 - (1./(2.*$e))*$aux20);
$q1 = (1. - $equad)*($aux11 - (1./(2.*$e))*$aux21);
$q2 = (1. - $equad)*($aux12 - (1./(2.*$e))*$aux22);
$q = (1. - $equad)*($aux1 - (1./(2.*$e))*$aux2);
$n = ($m1*$m1 - $m2*$m2)/($q2 - $q1);
$c = $m1*$m1 + $n*$q1;
$ro0 = $semi_eixo*sqrt($c - $n*$q0)/$n;
$teta = $n*($lon - $longO);
$ro = $semi_eixo*sqrt($c - $n*$q)/$n;
$x = $ro*sin($teta);
$y = $ro0 - $ro*cos($teta);
}
/********************************************************************
geo_to_cilequi
Autor : Julio Cesar Lima d'Alge fev-90
Conversao p/ Php : Luis Maurano 2008
Resumo : transforma coordenadas geodesicas em coordenadas cilindricas
equidistantes
Entradas :
lat coordenadas geodesicas (em radianos)
lon coordenadas geodesicas (em radianos)
r raio da Terra (em metros)
lon_mc meridiano origem
lat_padrao paralelo padrao
Saidas : x coordenadas cilindricas equidistantes
y
********************************************************************/
function geo_to_cilequi($lat,$lon,$r,$lat1,$longO,&$x,&$y)
{
$x = $r*($lon - $longO)*cos($lat1);
$y = $r*$lat;
}
/********************************************************************
cilequi_to_geo
Autor : Julio Cesar Lima d'Alge fev-90
Conversao p/ Php : Luis Maurano 2008
Resumo : transforma coordenadas cilindricas equidistantes em
coordenadas geodesicas
Entradas :
x coordenadas cilindricas equidistantes (em metros)
y coordenadas cilindricas equidistantes (em metros)
r raio da Terra (6371000 metros)
long0 meridiano origem (em radianos)
lat1 paralelo padrao (em radianos)
Saidas : lat lon coordenadas geodesicas (em radianos)
********************************************************************/
function cilequi_to_geo($x,$y,$r,$lat1,$longO,&$lat,&$lon)
{
$lat = $y/$r;
$lon = $longO + $x/($r*cos($lat1));
}
/********************************************************************
miller_to_geo
Autor : Julio Cesar Lima d'Alge mar-90
Conversao p/ Php : Luis Maurano 2008
Resumo : transforma coordenadas Miller em coordenadas geodesicas
Entradas :
x coordenadas Miller (em metros)
y coordenadas Miller (em metros)
r raio da Terra (6371000 metros)
long0 meridiano origem (em radianos)
Saidas : lat lon coordenadas geodesicas (em radianos)
********************************************************************/
function miller_to_geo($x,$y,$r,$longO,&$lat,&$lon)
{
// define PI
$pi = pi();
$e = 2.718281828;
$sinal = (int)(abs($y)/$y);
$lon = $longO + ($x/$r);
if ($sinal == 1)
$lat = 2.5*atan(pow($e,(0.8*$y/$r))) - 5.*$pi/8.;
else
{
$y *= $sinal;
$lat = 2.5*atan(pow($e,(0.8*$y/$r))) - 5.*$pi/8.;
$lat *= $sinal;
}
}
/********************************************************************
geo_to_miller
Autor : Julio Cesar Lima d'Alge mar-90
Conversao p/ Php : Luis Maurano 2008
Resumo : transforma coordenadas geodesicas em coordenadas Miller
Entradas :
lat coordenadas geodesicas (em radianos)
lon coordenadas geodesicas (em radianos)
r raio da Terra (6371000 metros)
long0 meridiano origem (em radianos)
Saidas : x y coordenadas Miller (em metros)
********************************************************************/
function geo_to_miller($lat,$lon,$r,$longO,&$x,&$y)
{
// define PI
$pi = pi();
$sinal = (int)(abs($lat)/$lat);
$x = $r*($lon - $longO);
if ($sinal == 1)
$y = ($r/0.8)*(log(tan($pi/4. + 0.4*$lat)));
else
{
$lat *= $sinal;
$y = ($r/0.8)*(log(tan($pi/4. + 0.4*$lat)));
$y *= $sinal;
}
}
/********************************************************************
datum1_to_datum2
Autor : Julio Cesar Lima d'Alge out-88
Conversao p/ Php : Luis Maurano 2008
Resumo : transformacao de datum horizontal
Entradas :
semi_eixo semi_eixo maior do elipsoide (datum entrada)
achat achatamento do elipsoide (datum entrada)
lat1 coordenadas geodesicas (em radianos)
lon1 coordenadas geodesicas (em radianos)
h1 0
semi_eixo2 semi_eixo maior do elipsoide (datum saida)
achat2 achatamento do elipsoide (datum saida)
lat2 coordenadas geodesicas (em radianos)
lon2 coordenadas geodesicas (em radianos)
Saidas : lat2 lon2 (em radianos)
********************************************************************/
function datum1_to_datum2 ($semi_eixo,$achat,$deltax,$deltay,$deltaz,$lat1,$lon1,$h1,$semi_eixo2,$achat2,$deltax2,$deltay2,$deltaz2,&$lat2,&$lon2,&$h2)
{
/* calcula coordenadas geocentricas cartesianas no datum 1 */
$equad1 = 2.* $achat - pow((double)1 * $achat,(double)2);
$n1 = $semi_eixo/sqrt((double)1-$equad1*pow(sin($lat1),(double)2));
$x1 = ($n1+$h1)*cos($lat1)*cos($lon1);
$y1 = ($n1+$h1)*cos($lat1)*sin($lon1);
$z1 = ($n1*(1-$equad1)+$h1)*sin($lat1);
/* calcula coordenadas geocentricas cartesianas no datum 2 */
$x2 = $x1 + ($deltax - $deltax2);
$y2 = $y1 + ($deltay - $deltay2);
$z2 = $z1 + ($deltaz - $deltaz2);
//echo "x1 = $x1 y1 = $y1 z1 = $z1 <br>";
//echo "x2 = $x2 y2 = $y2 z2 = $z2 <br>";
/* calcula coordenadas geodesicas no datum 2 */
$equad2 = 2.* $achat2 - pow((double)1 * $achat2,(double)2);
$lat2 = $lat1;
do
{
$n2 = $semi_eixo2/sqrt((double)1-$equad2*pow(sin($lat2),(double)2));
$lat2 = atan(($z2 + $n2*$equad2*sin($lat2))/sqrt($x2*$x2 + $y2*$y2));
$d = $semi_eixo2/sqrt((double)1-$equad2*pow(sin($lat2),(double)2))- $n2;
}
while (abs($d) > 0.00000000001);
$lon2 = atan($y2/$x2);
$h2 = $h1;
//echo "($semi_eixo,$achat,$deltax,$deltay,$deltaz,$lat1,$lon1,$h1,$semi_eixo2,$achat2,$deltax2,$deltay2,$deltaz2,&$lat2,&$lon2,&$h2)";
}
/********************************************************************
define_mer_cent
Autor : Julio Cesar Lima d'Alge fev-90
Conversao p/ Php : Luis Maurano 2008
Resumo : calcula meridianos centrais dos fusos UTM e GAUSS
a partir da longitude geodesica
Entradas : lon (em radianos)
Saidas : mc1_utm,mc2_utm,mc1_gauss,mc2_gauss
Obs. : os dois valores para UTM e Gauss referem-se 'as
situacoes de bordas de fusos
********************************************************************/
function define_mer_cent($lon,&$mc1_utm,&$mc2_utm)
{
// calcula MC
if ($lon != 0.)
$sinal = (int)($lon/abs($lon));
else
{
$mc1_utm = 3.;
$mc2_utm = -3.;
}
for ($k = $ind1 = 0, $ind2 = 6; abs($lon) > (double)$ind2; $k++, $ind1 = 6*$k, $ind2 = $ind1 + 6);
$mc = $ind1 + 3.;
if ((abs($lon) < (double)$mc) && ($lon != 0.))
{
$mc1_utm = $mc2_utm = (double)$mc;
}
if (((abs($lon) > (double)$mc)) && (abs($lon) != ((double)$mc + 3.)))
{
$mc1_utm = $mc2_utm = (double)$mc;
}
if ((abs($lon) == ((double)$mc + 3.)) && (abs($lon) != 180.))
{
$mc1_utm = (double)$mc;
$mc2_utm = (double)$mc + 6.;
}
if ($lon == 180.)
{
$mc1_utm = 177.;
$mc2_utm = -177.;
}
if ($lon == -180.)
{
$mc1_utm = -177.;
$mc2_utm = 177.;
}
if (abs($lon) == (double)$mc)
{
$mc1_utm = $mc2_utm = (double)$mc;
}
if (($lon != 0.) && (abs($lon) != 180))
{
$mc1_utm *= $sinal;
$mc2_utm *= $sinal;
}
}
?>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment