Completed
Branch 3.0-dev (e499db)
by Doug
02:12
created

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

This check looks for @param annotations where the type inferred by our type inference engine differs from the declared type.

It makes a suggestion as to what type it considers more descriptive.

Most often this is a case of a parameter that can be null in addition to its declared types.

Loading history...
226
     * @return int
227
     */
228 1
    public function distance(self $to): int {
229
230 1
        if ($this->refEll != $to->refEll) {
231
            throw new \RuntimeException('Source and destination co-ordinates are not using the same ellipsoid');
232
        }
233
234
        //Because this is a 2D grid, we can use simple Pythagoras
235 1
        $distanceX = $to->getX()-$this->getX();
236 1
        $distanceY = $to->getY()-$this->getY();
237
238 1
        return (int) round(pow(pow($distanceX, 2) + pow($distanceY, 2), 0.5));
239
240
    }
241
}
242