Subversion Repositories eFlore/Projets.eflore-projets

Rev

Details | 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                       OBLIQUE MERCATOR (HOTINE)
11
 
12
  PURPOSE:	Transforms input longitude and latitude to Easting and
13
  Northing for the Oblique Mercator projection.  The
14
  longitude and latitude must be in radians.  The Easting
15
  and Northing values will be returned in meters.
16
 
17
  PROGRAMMER              DATE
18
  ----------              ----
19
  T. Mittan		Mar, 1993
20
 
21
  ALGORITHM REFERENCES
22
 
23
  1.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
24
  Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
25
  State Government Printing Office, Washington D.C., 1987.
26
 
27
  2.  Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",
28
  U.S. Geological Survey Professional Paper 1453 , United State Government
29
  Printing Office, Washington D.C., 1989.
30
 * ***************************************************************************** */
31
 
32
class Proj4phpProjOmerc {
33
    /* Initialize the Oblique Mercator  projection
34
      ------------------------------------------ */
35
 
36
    public function init() {
37
        if( !$this->mode )
38
            $this->mode = 0;
39
        if( !$this->lon1 ) {
40
            $this->lon1 = 0;
41
            $this->mode = 1;
42
        }
43
        if( !$this->lon2 )
44
            $this->lon2 = 0;
45
        if( !$this->lat2 )
46
            $this->lat2 = 0;
47
 
48
        /* Place parameters in static storage for common use
49
          ------------------------------------------------- */
50
        $temp = $this->b / $this->a;
51
        $es = 1.0 - pow( $temp, 2 );
52
        $e = sqrt( $es );
53
 
54
        $this->sin_p20 = sin( $this->lat0 );
55
        $this->cos_p20 = cos( $this->lat0 );
56
 
57
        $this->con = 1.0 - $this->es * $this->sin_p20 * $this->sin_p20;
58
        $this->com = sqrt( 1.0 - $es );
59
        $this->bl = sqrt( 1.0 + $this->es * pow( $this->cos_p20, 4.0 ) / (1.0 - $es) );
60
        $this->al = $this->a * $this->bl * $this->k0 * $this->com / $this->con;
61
        if( abs( $this->lat0 ) < Proj4php::$common->EPSLN ) {
62
            $this->ts = 1.0;
63
            $this->d = 1.0;
64
            $this->el = 1.0;
65
        } else {
66
            $this->ts = Proj4php::$common->tsfnz( $this->e, $this->lat0, $this->sin_p20 );
67
            $this->con = sqrt( $this->con );
68
            $this->d = $this->bl * $this->com / ($this->cos_p20 * $this->con);
69
            if( ($this->d * $this->d - 1.0) > 0.0 ) {
70
                if( $this->lat0 >= 0.0 ) {
71
                    $this->f = $this->d + sqrt( $this->d * $this->d - 1.0 );
72
                } else {
73
                    $this->f = $this->d - sqrt( $this->d * $this->d - 1.0 );
74
                }
75
            } else {
76
                $this->f = $this->d;
77
            }
78
            $this->el = $this->f * pow( $this->ts, $this->bl );
79
        }
80
 
81
        //$this->longc=52.60353916666667;
82
 
83
        if( $this->mode != 0 ) {
84
            $this->g = .5 * ($this->f - 1.0 / $this->f);
85
            $this->gama = Proj4php::$common->asinz( sin( $this->alpha ) / $this->d );
86
            $this->longc = $this->longc - Proj4php::$common->asinz( $this->g * tan( $this->gama ) ) / $this->bl;
87
 
88
            /* Report parameters common to format B
89
              ------------------------------------- */
90
            //genrpt(azimuth * R2D,"Azimuth of Central Line:    ");
91
            //cenlon(lon_origin);
92
            // cenlat(lat_origin);
93
 
94
            $this->con = abs( $this->lat0 );
95
            if( ($this->con > Proj4php::$common->EPSLN) && (abs( $this->con - Proj4php::$common->HALF_PI ) > Proj4php::$common->EPSLN) ) {
96
                $this->singam = sin( $this->gama );
97
                $this->cosgam = cos( $this->gama );
98
 
99
                $this->sinaz = sin( $this->alpha );
100
                $this->cosaz = cos( $this->alpha );
101
 
102
                if( $this->lat0 >= 0 ) {
103
                    $this->u = ($this->al / $this->bl) * atan( sqrt( $this->d * $this->d - 1.0 ) / $this->cosaz );
104
                } else {
105
                    $this->u = -($this->al / $this->bl) * atan( sqrt( $this->d * $this->d - 1.0 ) / $this->cosaz );
106
                }
107
            } else {
108
                Proj4php::reportError( "omerc:Init:DataError" );
109
            }
110
        } else {
111
            $this->sinphi = sin( $this->at1 );
112
            $this->ts1 = Proj4php::$common->tsfnz( $this->e, $this->lat1, $this->sinphi );
113
            $this->sinphi = sin( $this->lat2 );
114
            $this->ts2 = Proj4php::$common->tsfnz( $this->e, $this->lat2, $this->sinphi );
115
            $this->h = pow( $this->ts1, $this->bl );
116
            $this->l = pow( $this->ts2, $this->bl );
117
            $this->f = $this->el / $this->h;
118
            $this->g = .5 * ($this->f - 1.0 / $this->f);
119
            $this->j = ($this->el * $this->el - $this->l * $this->h) / ($this->el * $this->el + $this->l * $this->h);
120
            $this->p = ($this->l - $this->h) / ($this->l + $this->h);
121
            $this->dlon = $this->lon1 - $this->lon2;
122
            if( $this->dlon < -Proj4php::$common->PI )
123
                $this->lon2 = $this->lon2 - 2.0 * Proj4php::$common->PI;
124
            if( $this->dlon > Proj4php::$common->PI )
125
                $this->lon2 = $this->lon2 + 2.0 * Proj4php::$common->PI;
126
            $this->dlon = $this->lon1 - $this->lon2;
127
            $this->longc = .5 * ($this->lon1 + $this->lon2) - atan( $this->j * tan( .5 * $this->bl * $this->dlon ) / $this->p ) / $this->bl;
128
            $this->dlon = Proj4php::$common->adjust_lon( $this->lon1 - $this->longc );
129
            $this->gama = atan( sin( $this->bl * $this->dlon ) / $this->g );
130
            $this->alpha = Proj4php::$common->asinz( $this->d * sin( $this->gama ) );
131
 
132
            /* Report parameters common to format A
133
              ------------------------------------- */
134
            if( abs( $this->lat1 - $this->lat2 ) <= Proj4php::$common->EPSLN ) {
135
                Proj4php::reportError( "omercInitDataError" );
136
                //return(202);
137
            } else {
138
                $this->con = abs( $this->lat1 );
139
            }
140
            if( ($this->con <= Proj4php::$common->EPSLN) || (abs( $this->con - Proj4php::$common->HALF_PI ) <= Proj4php::$common->EPSLN) ) {
141
                Proj4php::reportError( "omercInitDataError" );
142
                //return(202);
143
            } else {
144
                if( abs( abs( $this->lat0 ) - Proj4php::$common->HALF_PI ) <= Proj4php::$common->EPSLN ) {
145
                    Proj4php::reportError( "omercInitDataError" );
146
                    //return(202);
147
                }
148
            }
149
 
150
            $this->singam = sin( $this->gam );
151
            $this->cosgam = cos( $this->gam );
152
 
153
            $this->sinaz = sin( $this->alpha );
154
            $this->cosaz = cos( $this->alpha );
155
 
156
 
157
            if( $this->lat0 >= 0 ) {
158
                $this->u = ($this->al / $this->bl) * atan( sqrt( $this->d * $this->d - 1.0 ) / $this->cosaz );
159
            } else {
160
                $this->u = -($this->al / $this->bl) * atan( sqrt( $this->d * $this->d - 1.0 ) / $this->cosaz );
161
            }
162
        }
163
    }
164
 
165
    /* Oblique Mercator forward equations--mapping lat,long to x,y
166
      ---------------------------------------------------------- */
167
    public function forward( $p ) {
168
 
169
        /*
170
        $theta;   // angle
171
        $sin_phi;
172
        $cos_phi;  // sin and cos value
173
        $b;  // temporary values
174
        $c;
175
        $t;
176
        $tq; // temporary values
177
        $con;
178
        $n;
179
        $ml; // cone constant, small m
180
        $q;
181
        $us;
182
        $vl;
183
        $ul;
184
        $vs;
185
        $s;
186
        $dlon;
187
        $ts1;
188
        */
189
 
190
        $lon = $p->x;
191
        $lat = $p->y;
192
 
193
        /* Forward equations
194
          ----------------- */
195
        $sin_phi = sin( $lat );
196
        $dlon = Proj4php::$common->adjust_lon( $lon - $this->longc );
197
        $vl = sin( $this->bl * $dlon );
198
        if( abs( abs( $lat ) - Proj4php::$common->HALF_PI ) > Proj4php::$common->EPSLN ) {
199
            $ts1 = Proj4php::$common->tsfnz( $this->e, $lat, $sin_phi );
200
            $q = $this->el / (pow( $ts1, $this->bl ));
201
            $s = .5 * ($q - 1.0 / $q);
202
            $t = .5 * ($q + 1.0 / $q);
203
            $ul = ($s * $this->singam - $vl * $this->cosgam) / $t;
204
            $con = cos( $this->bl * $dlon );
205
            if( abs( con ) < .0000001 ) {
206
                $us = $this->al * $this->bl * $dlon;
207
            } else {
208
                $us = $this->al * atan( ($s * $this->cosgam + $vl * $this->singam) / $con ) / $this->bl;
209
                if( $con < 0 )
210
                    $us = $us + Proj4php::$common->PI * $this->al / $this->bl;
211
            }
212
        } else {
213
            if( $lat >= 0 ) {
214
                $ul = $this->singam;
215
            } else {
216
                $ul = -$this->singam;
217
            }
218
            $us = $this->al * $lat / $this->bl;
219
        }
220
        if( abs( abs( $ul ) - 1.0 ) <= Proj4php::$common->EPSLN ) {
221
            //alert("Point projects into infinity","omer-for");
222
            Proj4php::reportError( "omercFwdInfinity" );
223
            //return(205);
224
        }
225
        $vs = .5 * $this->al * log( (1.0 - $ul) / (1.0 + $ul) ) / $this->bl;
226
        $us = $us - $this->u;
227
        $p->x = $this->x0 + $vs * $this->cosaz + $us * $this->sinaz;
228
        $p->y = $this->y0 + $us * $this->cosaz - $vs * $this->sinaz;
229
 
230
        return $p;
231
    }
232
 
233
    /**
234
     *
235
     * @param type $p
236
     * @return type
237
     */
238
    public function inverse( $p ) {
239
        /*
240
        $delta_lon; /* Delta longitude (Given longitude - center
241
        $theta;  /* angle
242
        $delta_theta; /* adjusted longitude
243
        $sin_phi;
244
        $cos_phi; /* sin and cos value
245
        $b;  /* temporary values
246
        $c;
247
        $t;
248
        $tq; /* temporary values
249
        $con;
250
        $n;
251
        $ml; /* cone constant, small m
252
        $vs;
253
        $us;
254
        $q;
255
        $s;
256
        $ts1;
257
        $vl;
258
        $ul;
259
        $bs;
260
        $dlon;
261
        $flag;
262
        */
263
 
264
        /* Inverse equations
265
          ----------------- */
266
        $p->x -= $this->x0;
267
        $p->y -= $this->y0;
268
        #$flag = 0;
269
        $vs = $p->x * $this->cosaz - $p->y * $this->sinaz;
270
        $us = $p->y * $this->cosaz + $p->x * $this->sinaz;
271
        $us = $us + $this->u;
272
        $q = exp( -$this->bl * $vs / $this->al );
273
        $s = .5 * ($q - 1.0 / $q);
274
        $t = .5 * ($q + 1.0 / $q);
275
        $vl = sin( $this->bl * $us / $this->al );
276
        $ul = ($vl * $this->cosgam + $s * $this->singam) / $t;
277
        if( abs( abs( $ul ) - 1.0 ) <= Proj4php::$common->EPSLN ) {
278
            $lon = $this->longc;
279
            if( ul >= 0.0 ) {
280
                $lat = Proj4php::$common->HALF_PI;
281
            } else {
282
                $lat = -Proj4php::$common->HALF_PI;
283
            }
284
        } else {
285
            $con = 1.0 / $this->bl;
286
            $ts1 = pow( ($this->el / sqrt( (1.0 + $ul) / (1.0 - $ul) ) ), $con );
287
            $lat = Proj4php::$common->phi2z( $this->e, $ts1 );
288
            //if ($flag != 0)
289
            //return($flag);
290
            //~ con = cos($this->bl * us /al);
291
            $theta = $this->longc - atan2( ($s * $this->cosgam - $vl * $this->singam ), $con ) / $this->bl;
292
            $lon = Proj4php::$common->adjust_lon( $theta );
293
        }
294
        $p->x = $lon;
295
        $p->y = $lat;
296
        return $p;
297
    }
298
 
299
}
300
 
301
Proj4php::$proj['omerc'] = new Proj4phpProjOmerc();