Completed
Push — master ( 1a8ca7...b2e639 )
by Doug
01:57
created

TransverseMercator::getRefEll()   A

Complexity

Conditions 1
Paths 1

Size

Total Lines 4
Code Lines 2

Duplication

Lines 0
Ratio 0 %

Code Coverage

Tests 2
CRAP Score 1

Importance

Changes 1
Bugs 0 Features 1
Metric Value
c 1
b 0
f 1
dl 0
loc 4
ccs 2
cts 2
cp 1
rs 10
cc 1
eloc 2
nc 1
nop 0
crap 1
1
<?php
2
/**
3
 * PHPCoord
4
 * @package PHPCoord
5
 * @author Jonathan Stott
6
 * @author Doug Wright
7
 */
8
namespace PHPCoord;
9
10
/**
11
 * Abstract class representing a Tranverse Mercator Projection
12
 * @author Doug Wright
13
 * @package PHPCoord
14
 */
15
abstract class TransverseMercator
16
{
17
18
    /**
19
     * X
20
     * @var float
21
     */
22
    protected $x;
23
24
    /**
25
     * Y
26
     * @var float
27
     */
28
    protected $y;
29
30
    /**
31
     * h
32
     * @var float
33
     */
34
    protected $h;
35
36
    /**
37
     * Reference ellipsoid used in this datum
38
     * @var RefEll
39
     */
40
    protected $refEll;
41
42
    /**
43
     * Cartesian constructor.
44
     * @param float $x
45
     * @param float $y
46
     * @param float $h
47
     * @param RefEll $refEll
48
     */
49 29
    public function __construct($x, $y, $h, RefEll $refEll)
50
    {
51 29
        $this->setX($x);
52 29
        $this->setY($y);
53 29
        $this->setH($h);
54 29
        $this->setRefEll($refEll);
55 29
    }
56
57
    /**
58
     * String version of coordinate.
59
     * @return string
60
     */
61
    public function __toString()
62
    {
63
        return "({$this->x}, {$this->y}, {$this->h})";
64
    }
65
66
    /**
67
     * @return float
68
     */
69 1
    public function getX()
70
    {
71 1
        return $this->x;
72
    }
73
74
    /**
75
     * @param float $x
76
     */
77 29
    public function setX($x)
78
    {
79 29
        $this->x = $x;
80 29
    }
81
82
    /**
83
     * @return float
84
     */
85 1
    public function getY()
86
    {
87 1
        return $this->y;
88
    }
89
90
    /**
91
     * @param float $y
92
     */
93 29
    public function setY($y)
94
    {
95 29
        $this->y = $y;
96 29
    }
97
98
    /**
99
     * @return float
100
     */
101 1
    public function getH()
102
    {
103 1
        return $this->h;
104
    }
105
106
    /**
107
     * @param float $h
108
     */
109 29
    public function setH($h)
110
    {
111 29
        $this->h = $h;
112 29
    }
113
114
    /**
115
     * @return RefEll
116
     */
117 1
    public function getRefEll()
118
    {
119 1
        return $this->refEll;
120
    }
121
122
    /**
123
     * @param RefEll $refEll
124
     */
125 29
    public function setRefEll(RefEll $refEll)
126
    {
127 29
        $this->refEll = $refEll;
128 29
    }
129
130
131
    /**
132
     * Reference ellipsoid used by this projection
133
     * @return RefEll
134
     */
135
    abstract public function getReferenceEllipsoid();
136
137
    /**
138
     * Scale factor at central meridian
139
     * @return float
140
     */
141
    abstract public function getScaleFactor();
142
143
    /**
144
     * Northing of true origin
145
     * @return float
146
     */
147
    abstract public function getOriginNorthing();
148
149
    /**
150
     * Easting of true origin
151
     * @return float
152
     */
153
    abstract public function getOriginEasting();
154
155
    /**
156
     * Latitude of true origin
157
     * @return float
158
     */
159
    abstract public function getOriginLatitude();
160
161
    /**
162
     * Longitude of true origin
163
     * @return float
164
     */
165
    abstract public function getOriginLongitude();
166
167
    /**
168
     * Convert this grid reference into a latitude and longitude
169
     * Formula for transformation is taken from OS document
170
     * "A Guide to Coordinate Systems in Great Britain"
171
     *
172
     * @param float $N map coordinate (northing) of point to convert
173
     * @param float $E map coordinate (easting) of point to convert
174
     * @param float $N0 map coordinate (northing) of true origin
175
     * @param float $E0 map coordinate (easting) of true origin
176
     * @param float $phi0 map coordinate (latitude) of true origin
177
     * @param float $lambda0 map coordinate (longitude) of true origin and central meridian
178
     * @return LatLng
179
     */
180 7
    public function convertToLatitudeLongitude($N, $E, $N0, $E0, $phi0, $lambda0)
181
    {
182
183 7
        $phi0 = deg2rad($phi0);
184 7
        $lambda0 = deg2rad($lambda0);
185
186 7
        $refEll = $this->getReferenceEllipsoid();
187 7
        $F0 = $this->getScaleFactor();
188
189 7
        $a = $refEll->getMaj();
190 7
        $b = $refEll->getMin();
191 7
        $eSquared = $refEll->getEcc();
192 7
        $n = ($a - $b) / ($a + $b);
193 7
        $phiPrime = (($N - $N0) / ($a * $F0)) + $phi0;
194
195
        do {
196
            $M =
197 7
                ($b * $F0)
198 7
                * (((1 + $n + ((5 / 4) * $n * $n) + ((5 / 4) * $n * $n * $n))
199 7
                        * ($phiPrime - $phi0))
200 7
                    - (((3 * $n) + (3 * $n * $n) + ((21 / 8) * $n * $n * $n))
201 7
                        * sin($phiPrime - $phi0)
202 7
                        * cos($phiPrime + $phi0))
203 7
                    + ((((15 / 8) * $n * $n) + ((15 / 8) * $n * $n * $n))
204 7
                        * sin(2 * ($phiPrime - $phi0))
205 7
                        * cos(2 * ($phiPrime + $phi0)))
206 7
                    - (((35 / 24) * $n * $n * $n)
207 7
                        * sin(3 * ($phiPrime - $phi0))
208 7
                        * cos(3 * ($phiPrime + $phi0))));
209 7
            $phiPrime += ($N - $N0 - $M) / ($a * $F0);
210 7
        } while (($N - $N0 - $M) >= 0.00001);
211 7
        $v = $a * $F0 * pow(1 - $eSquared * pow(sin($phiPrime), 2), -0.5);
212
        $rho =
213
            $a
214
            * $F0
215 7
            * (1 - $eSquared)
216 7
            * pow(1 - $eSquared * pow(sin($phiPrime), 2), -1.5);
217 7
        $etaSquared = ($v / $rho) - 1;
218 7
        $VII = tan($phiPrime) / (2 * $rho * $v);
219
        $VIII =
220 7
            (tan($phiPrime) / (24 * $rho * pow($v, 3)))
221
            * (5
222 7
                + (3 * pow(tan($phiPrime), 2))
223 7
                + $etaSquared
224 7
                - (9 * pow(tan($phiPrime), 2) * $etaSquared));
225
        $IX =
226 7
            (tan($phiPrime) / (720 * $rho * pow($v, 5)))
227
            * (61
228 7
                + (90 * pow(tan($phiPrime), 2))
229 7
                + (45 * pow(tan($phiPrime), 2) * pow(tan($phiPrime), 2)));
230 7
        $X = (1 / cos($phiPrime)) / $v;
231
        $XI =
232 7
            ((1 / cos($phiPrime)) / (6 * $v * $v * $v))
233 7
            * (($v / $rho) + (2 * pow(tan($phiPrime), 2)));
234
        $XII =
235 7
            ((1 / cos($phiPrime)) / (120 * pow($v, 5)))
236
            * (5
237 7
                + (28 * pow(tan($phiPrime), 2))
238 7
                + (24 * pow(tan($phiPrime), 4)));
239
        $XIIA =
240 7
            ((1 / cos($phiPrime)) / (5040 * pow($v, 7)))
241
            * (61
242 7
                + (662 * pow(tan($phiPrime), 2))
243 7
                + (1320 * pow(tan($phiPrime), 4))
244 7
                + (720
245 7
                    * pow(tan($phiPrime), 6)));
246
        $phi =
247
            $phiPrime
248 7
            - ($VII * pow($E - $E0, 2))
249 7
            + ($VIII * pow($E - $E0, 4))
250 7
            - ($IX * pow($E - $E0, 6));
251
        $lambda =
252
            $lambda0
253 7
            + ($X * ($E - $E0))
254 7
            - ($XI * pow($E - $E0, 3))
255 7
            + ($XII * pow($E - $E0, 5))
256 7
            - ($XIIA * pow($E - $E0, 7));
257
258 7
        return new LatLng(rad2deg($phi), rad2deg($lambda), 0, $refEll);
259
    }
260
}
261