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