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                    VAN DER GRINTEN
11
 
12
  PURPOSE:	Transforms input Easting and Northing to longitude and
13
  latitude for the Van der Grinten projection.  The
14
  Easting and Northing must be in meters.  The longitude
15
  and latitude values will be returned in radians.
16
 
17
  PROGRAMMER              DATE
18
  ----------              ----
19
  T. Mittan		March, 1993
20
 
21
  This function was adapted from the Van Der Grinten projection code
22
  (FORTRAN) in the General Cartographic Transformation Package software
23
  which is available from the U.S. Geological Survey National Mapping Division.
24
 
25
  ALGORITHM REFERENCES
26
 
27
  1.  "New Equal-Area Map Projections for Noncircular Regions", John P. Snyder,
28
  The American Cartographer, Vol 15, No. 4, October 1988, pp. 341-355.
29
 
30
  2.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological
31
  Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United
32
  State Government Printing Office, Washington D.C., 1987.
33
 
34
  3.  "Software Documentation for GCTP General Cartographic Transformation
35
  Package", U.S. Geological Survey National Mapping Division, May 1982.
36
 * ***************************************************************************** */
37
 
38
class Proj4phpProjVandg {
39
 
40
    /* Initialize the Van Der Grinten projection
41
      ---------------------------------------- */
42
    public function init() {
43
        $this->R = 6370997.0; //Radius of earth
44
    }
45
 
46
    /**
47
     *
48
     * @param type $p
49
     * @return type
50
     */
51
    public function forward( $p ) {
52
 
53
        $lon = $p->x;
54
        $lat = $p->y;
55
 
56
        /* Forward equations
57
          ----------------- */
58
        $dlon = Proj4php::$common->adjust_lon( $lon - $this->long0 );
59
        $x;
60
        $y;
61
 
62
        if( abs( $lat ) <= Proj4php::$common->EPSLN ) {
63
            $x = $this->x0 + $this->R * $dlon;
64
            $y = $this->y0;
65
        }
66
        $theta = Proj4php::$common . asinz( 2.0 * abs( $lat / Proj4php::$common->PI ) );
67
        if( (abs( $dlon ) <= Proj4php::$common->EPSLN) || (abs( abs( $lat ) - Proj4php::$common->HALF_PI ) <= Proj4php::$common->EPSLN) ) {
68
            $x = $this->x0;
69
            if( $lat >= 0 ) {
70
                $y = $this->y0 + Proj4php::$common->PI * $this->R * tan( .5 * $theta );
71
            } else {
72
                $y = $this->y0 + Proj4php::$common->PI * $this->R * - tan( .5 * $theta );
73
            }
74
            //  return(OK);
75
        }
76
        $al = .5 * abs( (Proj4php::$common->PI / $dlon) - ($dlon / Proj4php::$common->PI) );
77
        $asq = $al * $al;
78
        $sinth = sin( $theta );
79
        $costh = cos( $theta );
80
 
81
        $g = $costh / ($sinth + $costh - 1.0);
82
        $gsq = $g * $g;
83
        $m = $g * (2.0 / $sinth - 1.0);
84
        $msq = $m * $m;
85
        $con = Proj4php::$common->PI * $this->R * ($al * ($g - $msq) + sqrt( $asq * ($g - $sq) * ($g - $msq) - ($msq + $asq) * ($gsq - $msq) )) / ($msq + $asq);
86
        if( $dlon < 0 ) {
87
            $con = -$con;
88
        }
89
        $x = $this->x0 + $con;
90
        $con = abs( $con / (Proj4php::$common->PI * $this->R) );
91
        if( $lat >= 0 ) {
92
            $y = $this->y0 + Proj4php::$common->PI * $this->R * sqrt( 1.0 - $con * $con - 2.0 * $al * $con );
93
        } else {
94
            $y = $this->y0 - Proj4php::$common->PI * $this->R * sqrt( 1.0 - $con * $con - 2.0 * $al * $con );
95
        }
96
 
97
        $p->x = $x;
98
        $p->y = $y;
99
 
100
        return $p;
101
    }
102
 
103
    /* Van Der Grinten inverse equations--mapping x,y to lat/long
104
      --------------------------------------------------------- */
105
 
106
    public function inverse( $p ) {
107
 
108
        /*
109
        $dlon;
110
        $xx;
111
        $yy;
112
        $xys;
113
        $c1;
114
        $c2;
115
        $c3;
116
        $al;
117
        $asq;
118
        $a1;
119
        $m1;
120
        $con;
121
        $th1;
122
        $d;
123
        */
124
 
125
        /* inverse equations
126
          ----------------- */
127
        $p->x -= $this->x0;
128
        $p->y -= $this->y0;
129
        $con = Proj4php::$common->PI * $this->R;
130
        $xx = $p->x / $con;
131
        $yy = $p->y / $con;
132
        $xys = $xx * $xx + $yy * $yy;
133
        $c1 = -abs( $yy ) * (1.0 + $xys);
134
        $c2 = $c1 - 2.0 * $yy * $yy + $xx * $xx;
135
        $c3 = -2.0 * $c1 + 1.0 + 2.0 * $yy * $yy + $xys * $xys;
136
        $d = $yy * $yy / $c3 + (2.0 * $c2 * $c2 * $c2 / $c3 / $c3 / $c3 - 9.0 * $c1 * $c2 / $c3 / $c3) / 27.0;
137
        $a1 = ($c1 - $c2 * $c2 / 3.0 / $c3) / $c3;
138
        $m1 = 2.0 * sqrt( -$a1 / 3.0 );
139
        $con = ((3.0 * $d) / $a1) / $m1;
140
        if( abs( $con ) > 1.0 ) {
141
            if( $con >= 0.0 ) {
142
                $con = 1.0;
143
            } else {
144
                $con = -1.0;
145
            }
146
        }
147
        $th1 = acos( $con ) / 3.0;
148
        if( $p->$y >= 0 ) {
149
            $lat = (-$m1 * cos( $th1 + Proj4php::$common->PI / 3.0 ) - $c2 / 3.0 / $c3) * Proj4php::$common->PI;
150
        } else {
151
            $lat = -(-m1 * cos( $th1 + Proj4php::$common->PI / 3.0 ) - $c2 / 3.0 / $c3) * Proj4php::$common->PI;
152
        }
153
 
154
        if( abs( $xx ) < Proj4php::$common->EPSLN ) {
155
            $lon = $this->$long0;
156
        }
157
        $lon = Proj4php::$common->adjust_lon( $this->long0 + Proj4php::$common->PI * ($xys - 1.0 + sqrt( 1.0 + 2.0 * ($xx * $xx - $yy * $yy) + $xys * $xys )) / 2.0 / $xx );
158
 
159
        $p->x = $lon;
160
        $p->y = $lat;
161
 
162
        return $p;
163
    }
164
 
165
}
166
 
167
Proj4php::$proj['vandg'] = new Proj4phpProjVandg();