Subversion Repositories eFlore/Projets.eflore-projets

Rev

Rev 1001 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
1001 delphine 1
<?php
2
/**
3
 * Author : Julien Moquet
4
 *
5
 * Inspired by Proj4php from Mike Adair madairATdmsolutions.ca
6
 *                      and Richard Greenwood rich@greenwoodma$p->com
7
 * License: LGPL as per: http://www.gnu.org/copyleft/lesser.html
8
 */
9
/*******************************************************************************
10
  NAME                       SWISS OBLIQUE MERCATOR
11
 
12
  PURPOSE:	Swiss projection.
13
  WARNING:  X and Y are inverted (weird) in the swiss coordinate system. Not
14
  here, since we want X to be horizontal and Y vertical.
15
 
16
  ALGORITHM REFERENCES
17
  1. "Formules et constantes pour le Calcul pour la
18
  projection cylindrique conforme à axe oblique et pour la transformation entre
19
  des systèmes de référence".
20
  http://www.swisstopo.admin.ch/internet/swisstopo/fr/home/topics/survey/sys/refsys/switzerland.parsysrelated1.31216.downloadList.77004.DownloadFile.tmp/swissprojectionfr.pdf
21
 
22
*******************************************************************************/
23
 
24
class Proj4phpProjSomerc {
25
 
26
    /**
27
     *
28
     */
29
    public function init() {
30
        $phy0 = $this->lat0;
31
        $this->lambda0 = $this->long0;
32
        $sinPhy0 = sin( $phy0 );
33
        $semiMajorAxis = $this->a;
34
        $invF = $this->rf;
35
        $flattening = 1 / $invF;
36
        $e2 = 2 * $flattening - pow( $flattening, 2 );
37
        $e = $this->e = sqrt( $e2 );
38
        $this->R = $this->k0 * $semiMajorAxis * sqrt( 1 - $e2 ) / (1 - $e2 * pow( $sinPhy0, 2.0 ));
39
        $this->alpha = sqrt( 1 + $e2 / (1 - $e2) * pow( cos( $phy0 ), 4.0 ) );
40
        $this->b0 = asin( $sinPhy0 / $this->alpha );
41
        $this->K = log( tan( $PI / 4.0 + $this->b0 / 2.0 ) )
42
                  - $this->alpha
43
                  * log( tan( $PI / 4.0 + $phy0 / 2.0 ) )
44
                  + $this->alpha
45
                  * $e / 2
46
                  * log( (1 + $e * $sinPhy0)
47
                            / (1 - $e * $sinPhy0) );
48
    }
49
 
50
    /**
51
     *
52
     * @param type $p
53
     * @return type
54
     */
55
    public function forward( $p ) {
56
        $Sa1 = log( tan( $PI / 4.0 - $p->y / 2.0 ) );
57
        $Sa2 = $this->e / 2.0
58
                  * log( (1 + $this->e * sin( $p->y ))
59
                            / (1 - $this->e * sin( $p->y )) );
60
        $S = -$this->alpha * ($Sa1 + $Sa2) + $this->K;
61
 
62
        // spheric latitude
63
        $b = 2.0 * (atan( exp( $S ) ) - proj4phpCommon::PI / 4.0);
64
 
65
        // spheric longitude
66
        $I = $this->alpha * ($p->x - $this->lambda0);
67
 
68
        // psoeudo equatorial rotation
69
        $rotI = atan( sin( $I )
70
                  / (sin( $this->b0 ) * tan( $b ) +
71
                  cos( $this->b0 ) * cos( $I )) );
72
 
73
        $rotB = asin( cos( $this->b0 ) * sin( $b ) -
74
                  sin( $this->b0 ) * cos( $b ) * cos( $I ) );
75
 
76
        $p->y = $this->R / 2.0
77
                  * log( (1 + sin( $rotB )) / (1 - sin( $rotB )) )
78
                  + $this->y0;
79
 
80
        $p->x = $this->R * $rotI + $this->x0;
81
 
82
        return $p;
83
    }
84
 
85
    /**
86
     *
87
     * @param type $p
88
     * @return type
89
     */
90
    public function inverse( $p ) {
91
 
92
        $Y = $p->x - $this->x0;
93
        $X = $p->y - $this->y0;
94
 
95
        $rotI = $Y / $this->R;
96
        $rotB = 2 * (atan( exp( $X / $this->R ) ) - $PI / 4.0);
97
 
98
        $b = asin( cos( $this->b0 ) * sin( $rotB )
99
                  + sin( $this->b0 ) * cos( $rotB ) * cos( $rotI ) );
100
        $I = atan( sin( $rotI )
101
                  / (cos( $this->b0 ) * cos( $rotI ) - sin( $this->b0 )
102
                  * tan( $rotB )) );
103
 
104
        $lambda = $this->lambda0 + $I / $this->alpha;
105
 
106
        $S = 0.0;
107
        $phy = $b;
108
        $prevPhy = -1000.0;
109
        $iteration = 0;
110
        while( abs( $phy - $prevPhy ) > 0.0000001 ) {
111
            if( ++$iteration > 20 ) {
112
                Proj4php::reportError( "omercFwdInfinity" );
113
                return;
114
            }
115
            //S = log(tan(PI / 4.0 + phy / 2.0));
116
            $S = 1.0
117
                      / $this->alpha
118
                      * (log( tan( $PI / 4.0 + $b / 2.0 ) ) - $this->K)
119
                      + $this->e
120
                      * log( tan( $PI / 4.0
121
                                          + asin( $this->e * sin( $phy ) )
122
                                          / 2.0 ) );
123
            $prevPhy = $phy;
124
            $phy = 2.0 * atan( exp( $S ) ) - $PI / 2.0;
125
        }
126
 
127
        $p->x = $lambda;
128
        $p->y = $phy;
129
 
130
        return $p;
131
    }
132
 
133
}
134
 
135
Proj4php::$proj['somerc'] = new Proj4phpProjSomerc();