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                            TRANSVERSE MERCATOR
11
 
12
  PURPOSE:	Transforms input longitude and latitude to Easting and
13
  Northing for the Transverse Mercator projection.  The
14
  longitude and latitude must be in radians.  The Easting
15
  and Northing values will be returned in meters.
16
 
17
  ALGORITHM REFERENCES
18
 
19
  1.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
20
  Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
21
  State Government Printing Office, Washington D.C., 1987.
22
 
23
  2.  Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
24
  U.S. Geological Survey Professional Paper 1453 , United State Government
25
  Printing Office, Washington D.C., 1989.
26
*******************************************************************************/
27
 
28
/**
29
  Initialize Transverse Mercator projection
30
 */
31
class Proj4phpProjTmerc {
32
 
33
    private $e0, $e1, $e2, $e3, $ml0;
34
 
35
    /**
36
     *
37
     */
38
    public function init() {
39
 
40
        $this->e0 = Proj4php::$common->e0fn( $this->es );
41
        $this->e1 = Proj4php::$common->e1fn( $this->es );
42
        $this->e2 = Proj4php::$common->e2fn( $this->es );
43
        $this->e3 = Proj4php::$common->e3fn( $this->es );
44
        $this->ml0 = $this->a * Proj4php::$common->mlfn( $this->e0, $this->e1, $this->e2, $this->e3, $this->lat0 );
45
    }
46
 
47
    /**
48
      Transverse Mercator Forward  - long/lat to x/y
49
      long/lat in radians
50
     */
51
    public function forward( $p ) {
52
 
53
        $lon = $p->x;
54
        $lat = $p->y;
55
 
56
        $delta_lon = Proj4php::$common->adjust_lon( $lon - $this->long0 ); // Delta longitude
57
        #$con = 0;    // cone constant
58
        #$x = 0;
59
        #$y = 0;
60
        $sin_phi = sin( $lat );
61
        $cos_phi = cos( $lat );
62
 
63
        if( isset($this->sphere) && $this->sphere === true ) { /* spherical form */
64
            $b = $cos_phi * sin( $delta_lon );
65
            if( (abs( abs( $b ) - 1.0 )) < .0000000001 ) {
66
                Proj4php::reportError( "tmerc:forward: Point projects into infinity" );
67
                return(93);
68
            } else {
69
                $x = .5 * $this->a * $this->k0 * log( (1.0 + $b) / (1.0 - $b) );
70
                $con = acos( $cos_phi * cos( $delta_lon ) / sqrt( 1.0 - $b * $b ) );
71
                if( $lat < 0 )
72
                    $con = - $con;
73
                $y = $this->a * $this->k0 * ($con - $this->lat0);
74
            }
75
        } else {
76
            $al = $cos_phi * $delta_lon;
77
            $als = pow( $al, 2 );
78
            $c = $this->ep2 * pow( $cos_phi, 2 );
79
            $tq = tan( $lat );
80
            $t = pow( $tq, 2 );
81
            $con = 1.0 - $this->es * pow( $sin_phi, 2 );
82
            $n = $this->a / sqrt( $con );
83
 
84
            $ml = $this->a * Proj4php::$common->mlfn( $this->e0, $this->e1, $this->e2, $this->e3, $lat );
85
 
86
            $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;
87
            $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;
88
        }
89
 
90
        $p->x = $x;
91
        $p->y = $y;
92
 
93
        return $p;
94
    }
95
 
96
    /**
97
      Transverse Mercator Inverse  -  x/y to long/lat
98
     */
99
    public function inverse( $p ) {
100
 
101
        #$phi;  /* temporary angles       */
102
        #$delta_phi; /* difference between longitudes    */
103
        $max_iter = 6;      /* maximun number of iterations */
104
 
105
        if( isset($this->sphere) && $this->sphere === true ) { /* spherical form */
106
            $f = exp( $p->x / ($this->a * $this->k0) );
107
            $g = .5 * ($f - 1 / $f);
108
            $temp = $this->lat0 + $p->y / ($this->a * $this->k0);
109
            $h = cos( $temp );
110
            $con = sqrt( (1.0 - $h * $h) / (1.0 + $g * $g) );
111
            $lat = Proj4php::$common->asinz( $con );
112
            if( $temp < 0 )
113
                $lat = -$lat;
114
            if( ($g == 0) && ($h == 0) ) {
115
                $lon = $this->long0;
116
            } else {
117
                $lon = Proj4php::$common->adjust_lon( atan2( $g, $h ) + $this->long0 );
118
            }
119
        } else {    // ellipsoidal form
120
            $x = $p->x - $this->x0;
121
            $y = $p->y - $this->y0;
122
 
123
            $con = ($this->ml0 + $y / $this->k0) / $this->a;
124
            $phi = $con;
125
 
126
            for( $i = 0; true; $i++ ) {
127
                $delta_phi = (($con + $this->e1 * sin( 2.0 * $phi ) - $this->e2 * sin( 4.0 * $phi ) + $this->e3 * sin( 6.0 * $phi )) / $this->e0) - $phi;
128
                $phi += $delta_phi;
129
                if( abs( $delta_phi ) <= Proj4php::$common->EPSLN )
130
                    break;
131
                if( $i >= $max_iter ) {
132
                    Proj4php::reportError( "tmerc:inverse: Latitude failed to converge" );
133
                    return(95);
134
                }
135
            } // for()
136
            if( abs( $phi ) < Proj4php::$common->HALF_PI ) {
137
                // sincos(phi, &sin_phi, &cos_phi);
138
                $sin_phi = sin( $phi );
139
                $cos_phi = cos( $phi );
140
                $tan_phi = tan( $phi );
141
                $c = $this->ep2 * pow( $cos_phi, 2 );
142
                $cs = pow( $c, 2 );
143
                $t = pow( $tan_phi, 2 );
144
                $ts = pow( $t, 2 );
145
                $con = 1.0 - $this->es * pow( $sin_phi, 2 );
146
                $n = $this->a / sqrt( $con );
147
                $r = $n * (1.0 - $this->es) / $con;
148
                $d = $x / ($n * $this->k0);
149
                $ds = pow( $d, 2 );
150
                $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)));
151
                $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) );
152
            } else {
153
                $lat = Proj4php::$common->HALF_PI * Proj4php::$common->sign( $y );
154
                $lon = $this->long0;
155
            }
156
        }
157
 
158
        $p->x = $lon;
159
        $p->y = $lat;
160
 
161
        return $p;
162
    }
163
 
164
}
165
 
166
Proj4php::$proj['tmerc'] = new Proj4phpProjTmerc();