| 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(); |