Created
January 16, 2020 19:27
-
-
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.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| <?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