Completed
Push — master ( 4c6b8b...d48a1a )
by Doug
03:34
created

LatLng::transformDatum()   A

Complexity

Conditions 2
Paths 2

Size

Total Lines 10
Code Lines 5

Duplication

Lines 0
Ratio 0 %

Code Coverage

Tests 6
CRAP Score 2

Importance

Changes 0
Metric Value
cc 2
eloc 5
nc 2
nop 8
dl 0
loc 10
ccs 6
cts 6
cp 1
crap 2
rs 10
c 0
b 0
f 0

How to fix   Many Parameters   

Many Parameters

Methods with many parameters are not only hard to understand, but their parameters also often become inconsistent when you need more, or different data.

There are several approaches to avoid long parameter lists:

1
<?php
2
3
declare(strict_types=1);
4
5
namespace PHPCoord;
6
7
/**
8
 * Latitude/Longitude reference.
9
 *
10
 * @author Jonathan Stott
11
 * @author Doug Wright
12
 */
13
class LatLng
14
{
15
    /**
16
     * Latitude.
17
     *
18
     * @var float
19
     */
20
    protected $lat;
21
22
    /**
23
     * Longitude.
24
     *
25
     * @var float
26
     */
27
    protected $lng;
28
29
    /**
30
     * Height.
31
     *
32
     * @var float
33
     */
34
    protected $h;
35
36
    /**
37
     * Reference ellipsoid the co-ordinates are from.
38
     *
39
     * @var RefEll
40
     */
41
    protected $refEll;
42
43
    /**
44
     * Create a new LatLng object from the given latitude and longitude.
45
     *
46
     * @param float  $lat
47
     * @param float  $lng
48
     * @param int    $height
49
     * @param RefEll $refEll
50
     */
51 31
    public function __construct(float $lat, float $lng, int $height, RefEll $refEll)
52
    {
53 31
        $this->lat = $lat;
54 31
        $this->lng = $lng;
55 31
        $this->h = $height;
56 31
        $this->refEll = $refEll;
57 31
    }
58
59
    /**
60
     * Return a string representation of this LatLng object.
61
     *
62
     * @return string
63
     */
64 13
    public function __toString(): string
65
    {
66 13
        return "({$this->getLat()}, {$this->getLng()})";
67
    }
68
69
    /**
70
     * @return float
71
     */
72 24
    public function getLat(): float
73
    {
74 24
        return round($this->lat, 5);
75
    }
76
77
    /**
78
     * @return float
79
     */
80 24
    public function getLng(): float
81
    {
82 24
        return round($this->lng, 5);
83
    }
84
85
    /**
86
     * @return int
87
     */
88 10
    public function getH(): int
89
    {
90 10
        return $this->h;
91
    }
92
93
    /**
94
     * @return RefEll
95
     */
96 9
    public function getRefEll(): RefEll
97
    {
98 9
        return $this->refEll;
99
    }
100
101
    /**
102
     * Calculate the surface distance between this LatLng object and the one
103
     * passed in as a parameter.
104
     *
105
     * @param LatLng $to a LatLng object to measure the surface distance to
106
     *
107
     * @return int distance in metres
108
     */
109 2
    public function distance(self $to): int
110
    {
111 2
        if ($this->refEll != $to->refEll) {
112 1
            throw new \RuntimeException('Source and destination co-ordinates are not using the same ellipsoid');
113
        }
114
115
        //Mean radius definition taken from Wikipedia
116 1
        $er = ((2 * $this->refEll->getMaj()) + $this->refEll->getMin()) / 3;
117
118 1
        $latFrom = deg2rad($this->lat);
119 1
        $latTo = deg2rad($to->lat);
120 1
        $lngFrom = deg2rad($this->lng);
121 1
        $lngTo = deg2rad($to->lng);
122
123 1
        $d = acos(sin($latFrom) * sin($latTo) + cos($latFrom) * cos($latTo) * cos($lngTo - $lngFrom)) * $er;
124
125 1
        return (int) round($d);
126
    }
127
128
    /**
129
     * Convert this LatLng object to OSGB36 datum.
130
     * Reference values for transformation are taken from OS document
131
     * "A Guide to Coordinate Systems in Great Britain".
132
     *
133
     * @return self
134
     */
135 4
    public function toOSGB36(): self
136
    {
137 4
        if ($this->refEll == RefEll::airy1830()) {
138 2
            return $this;
139
        }
140
141 2
        $asWGS84 = $this->toWGS84();
142
143 1
        $tx = -446.448;
144 1
        $ty = 125.157;
145 1
        $tz = -542.060;
146 1
        $s = 0.0000204894;
147 1
        $rx = deg2rad(-0.1502 / 3600);
148 1
        $ry = deg2rad(-0.2470 / 3600);
149 1
        $rz = deg2rad(-0.8421 / 3600);
150
151 1
        return $asWGS84->transformDatum(RefEll::airy1830(), $tx, $ty, $tz, $s, $rx, $ry, $rz);
152
    }
153
154
    /**
155
     * Convert this LatLng object to ED50 datum.
156
     * Reference values for transformation are taken from http://www.globalmapper.com/helpv9/datum_list.htm .
157
     *
158
     * @return self
159
     */
160 1
    public function toED50(): self
161
    {
162 1
        if ($this->refEll == RefEll::international1924()) {
163
            return $this;
164
        }
165
166 1
        $asWGS84 = $this->toWGS84();
167
168 1
        $tx = 87;
169 1
        $ty = 98;
170 1
        $tz = 121;
171 1
        $s = 0;
172 1
        $rx = deg2rad(0);
173 1
        $ry = deg2rad(0);
174 1
        $rz = deg2rad(0);
175
176 1
        return $asWGS84->transformDatum(RefEll::international1924(), $tx, $ty, $tz, $s, $rx, $ry, $rz);
177
    }
178
179
    /**
180
     * Convert this LatLng object to NAD27 datum.
181
     * Reference values for transformation are taken from Wikipedia.
182
     *
183
     * @return self
184
     */
185 1
    public function toNAD27(): self
186
    {
187 1
        if ($this->refEll == RefEll::clarke1866()) {
188
            return $this;
189
        }
190
191 1
        $asWGS84 = $this->toWGS84();
192
193 1
        $tx = 8;
194 1
        $ty = -160;
195 1
        $tz = -176;
196 1
        $s = 0;
197 1
        $rx = deg2rad(0);
198 1
        $ry = deg2rad(0);
199 1
        $rz = deg2rad(0);
200
201 1
        return $asWGS84->transformDatum(RefEll::clarke1866(), $tx, $ty, $tz, $s, $rx, $ry, $rz);
202
    }
203
204
    /**
205
     * Convert this LatLng object to Ireland 1975 datum.
206
     * Reference values for transformation are taken from OSI document
207
     * "Making maps compatible with GPS".
208
     *
209
     * @return self
210
     */
211 1
    public function toIreland1975(): self
212
    {
213 1
        if ($this->refEll == RefEll::airyModified()) {
214
            return $this;
215
        }
216
217 1
        $asWGS84 = $this->toWGS84();
218
219 1
        $tx = -482.530;
220 1
        $ty = 130.596;
221 1
        $tz = -564.557;
222 1
        $s = -0.00000815;
223 1
        $rx = deg2rad(-1.042 / 3600);
224 1
        $ry = deg2rad(-0.214 / 3600);
225 1
        $rz = deg2rad(-0.631 / 3600);
226
227 1
        return $asWGS84->transformDatum(RefEll::airyModified(), $tx, $ty, $tz, $s, $rx, $ry, $rz);
228
    }
229
230
    /**
231
     * Convert this LatLng object to WGS84 datum.
232
     *
233
     * @return self
234
     */
235 15
    public function toWGS84(): self
236
    {
237 15
        switch ($this->refEll) {
238 15
            case RefEll::wgs84():
239 8
                return $this; //do nothing
240
241 7
            case RefEll::airy1830(): // values from OSGB document "A Guide to Coordinate Systems in Great Britain"
242 2
                $tx = 446.448;
243 2
                $ty = -125.157;
244 2
                $tz = 542.060;
245 2
                $s = -0.0000204894;
246 2
                $rx = deg2rad(0.1502 / 3600);
247 2
                $ry = deg2rad(0.2470 / 3600);
248 2
                $rz = deg2rad(0.8421 / 3600);
249 2
                break;
250
251 5
            case RefEll::airyModified(): // values from OSI document "Making maps compatible with GPS"
252 1
                $tx = 482.530;
253 1
                $ty = -130.596;
254 1
                $tz = 564.557;
255 1
                $s = 0.00000815;
256 1
                $rx = deg2rad(1.042 / 3600);
257 1
                $ry = deg2rad(0.214 / 3600);
258 1
                $rz = deg2rad(0.631 / 3600);
259 1
                break;
260
261 4
            case RefEll::clarke1866(): // assumes NAD27, values from Wikipedia
262 1
                $tx = -8;
263 1
                $ty = 160;
264 1
                $tz = 176;
265 1
                $s = 0;
266 1
                $rx = deg2rad(0);
267 1
                $ry = deg2rad(0);
268 1
                $rz = deg2rad(0);
269 1
                break;
270
271 3
            case RefEll::international1924(): // assumes ED50, values from http://www.globalmapper.com/helpv9/datum_list.htm
272 1
                $tx = -87;
273 1
                $ty = -98;
274 1
                $tz = -121;
275 1
                $s = 0;
276 1
                $rx = deg2rad(0);
277 1
                $ry = deg2rad(0);
278 1
                $rz = deg2rad(0);
279 1
                break;
280
281 2
            case RefEll::bessel1841(): // assumes Germany, values from Wikipedia
282
                $tx = 582;
283
                $ty = -105;
284
                $tz = -414;
285
                $s = 0.0000083;
286
                $rx = deg2rad(1.04 / 3600);
287
                $ry = deg2rad(0.35 / 3600);
288
                $rz = deg2rad(-3.08 / 3600);
289
                break;
290
291
            default:
292 2
                throw new \RuntimeException('Transform parameters not known for this ellipsoid');
293
        }
294
295 5
        return $this->transformDatum(RefEll::wgs84(), $tx, $ty, $tz, $s, $rx, $ry, $rz);
296
    }
297
298
    /**
299
     * Transform co-ordinates from one datum to another using a Helmert transformation.
300
     *
301
     * @param RefEll $toRefEll
302
     * @param float  $tranX
303
     * @param float  $tranY
304
     * @param float  $tranZ
305
     * @param float  $scale
306
     * @param float  $rotX     rotation about x-axis in seconds
307
     * @param float  $rotY     rotation about y-axis in seconds
308
     * @param float  $rotZ     rotation about z-axis in seconds
309
     *
310
     * @return self
311
     */
312 15
    public function transformDatum(RefEll $toRefEll, float $tranX, float $tranY, float $tranZ, float $scale, float $rotX, float $rotY, float $rotZ): self
313
    {
314 15
        if ($this->refEll == $toRefEll) {
315 6
            return $this;
316
        }
317
318 9
        $cartesian = Cartesian::fromLatLong($this);
319 9
        $newCartesian = $cartesian->transformDatum($toRefEll, $tranX, $tranY, $tranZ, $scale, $rotX, $rotY, $rotZ);
320
321 9
        return $newCartesian->toLatitudeLongitude();
322
    }
323
324
    /**
325
     * Convert this LatLng object into an OSGB grid reference. Note that this
326
     * function does not take into account the bounds of the OSGB grid -
327
     * beyond the bounds of the OSGB grid, the resulting OSRef object has no
328
     * meaning.
329
     *
330
     * Reference values for transformation are taken from OS document
331
     * "A Guide to Coordinate Systems in Great Britain".
332
     *
333
     * @return OSRef
334
     */
335 3
    public function toOSRef(): OSRef
336
    {
337 3
        $asOSGB36 = $this->toOSGB36();
338
339 2
        $OSGB = new OSRef(0, 0); //dummy to get reference data
340 2
        $scale = $OSGB->getScaleFactor();
341 2
        $N0 = $OSGB->getOriginNorthing();
342 2
        $E0 = $OSGB->getOriginEasting();
343 2
        $phi0 = $OSGB->getOriginLatitude();
344 2
        $lambda0 = $OSGB->getOriginLongitude();
345
346 2
        $coords = $asOSGB36->toTransverseMercatorEastingNorthing($scale, $E0, $N0, $phi0, $lambda0);
347
348 2
        return new OSRef($coords['E'], $coords['N'], $this->h);
0 ignored issues
show
Bug introduced by
$this->h of type double is incompatible with the type integer expected by parameter $z of PHPCoord\OSRef::__construct(). ( Ignorable by Annotation )

If this is a false-positive, you can also ignore this issue in your code via the ignore-type  annotation

348
        return new OSRef($coords['E'], $coords['N'], /** @scrutinizer ignore-type */ $this->h);
Loading history...
349
    }
350
351
    /**
352
     * Convert this LatLng object into an ITM grid reference.
353
     *
354
     * @return ITMRef
355
     */
356 1
    public function toITMRef()
357
    {
358 1
        $asWGS84 = $this->toWGS84();
359
360 1
        $ITM = new ITMRef(0, 0); //dummy to get reference data
361 1
        $scale = $ITM->getScaleFactor();
362 1
        $N0 = $ITM->getOriginNorthing();
363 1
        $E0 = $ITM->getOriginEasting();
364 1
        $phi0 = $ITM->getOriginLatitude();
365 1
        $lambda0 = $ITM->getOriginLongitude();
366
367 1
        $coords = $asWGS84->toTransverseMercatorEastingNorthing($scale, $E0, $N0, $phi0, $lambda0);
368
369 1
        return new ITMRef($coords['E'], $coords['N'], $this->h);
0 ignored issues
show
Bug introduced by
$this->h of type double is incompatible with the type integer expected by parameter $z of PHPCoord\ITMRef::__construct(). ( Ignorable by Annotation )

If this is a false-positive, you can also ignore this issue in your code via the ignore-type  annotation

369
        return new ITMRef($coords['E'], $coords['N'], /** @scrutinizer ignore-type */ $this->h);
Loading history...
370
    }
371
372
    /**
373
     * Convert a WGS84 latitude and longitude to an UTM reference.
374
     *
375
     * Reference values for transformation are taken from OS document
376
     * "A Guide to Coordinate Systems in Great Britain".
377
     *
378
     * @return UTMRef
379
     */
380 4
    public function toUTMRef(): UTMRef
381
    {
382 4
        $asWGS84 = $this->toWGS84();
383
384 4
        $longitudeZone = (int) (($asWGS84->getLng() + 180) / 6) + 1;
385
386
        // Special zone for Norway
387 4
        if ($asWGS84->getLat() >= 56 && $asWGS84->getLat() < 64 && $asWGS84->getLng() >= 3 && $asWGS84->getLng() < 12) {
388
            $longitudeZone = 32;
389 4
        } elseif ($asWGS84->getLat() >= 72 && $asWGS84->getLat() < 84) { // Special zones for Svalbard
390
            if ($asWGS84->getLng() >= 0 && $asWGS84->getLng() < 9) {
391
                $longitudeZone = 31;
392
            } elseif ($asWGS84->getLng() >= 9 && $asWGS84->getLng() < 21) {
393
                $longitudeZone = 33;
394
            } elseif ($asWGS84->getLng() >= 21 && $asWGS84->getLng() < 33) {
395
                $longitudeZone = 35;
396
            } elseif ($asWGS84->getLng() >= 33 && $asWGS84->getLng() < 42) {
397
                $longitudeZone = 37;
398
            }
399
        }
400
401 4
        $UTMZone = $this->getUTMLatitudeZoneLetter($asWGS84->getLat());
402
403 2
        $UTM = new UTMRef(0, 0, 0, $UTMZone, $longitudeZone); //dummy to get reference data
404 2
        $scale = $UTM->getScaleFactor();
405 2
        $N0 = $UTM->getOriginNorthing();
406 2
        $E0 = $UTM->getOriginEasting();
407 2
        $phi0 = $UTM->getOriginLatitude();
408 2
        $lambda0 = $UTM->getOriginLongitude();
409
410 2
        $coords = $asWGS84->toTransverseMercatorEastingNorthing($scale, $E0, $N0, $phi0, $lambda0);
411
412 2
        if ($asWGS84->getLat() < 0) {
413 1
            $coords['N'] += 10000000;
414
        }
415
416 2
        return new UTMRef($coords['E'], $coords['N'], $asWGS84->getH(), $UTMZone, $longitudeZone);
417
    }
418
419
    /**
420
     * Work out the UTM latitude zone from the latitude.
421
     *
422
     * @param float $latitude
423
     *
424
     * @return string
425
     */
426 4
    private function getUTMLatitudeZoneLetter($latitude): string
427
    {
428 4
        if ($latitude < -80 || $latitude > 84) {
429 2
            throw new \OutOfRangeException('UTM zones do not apply in polar regions');
430
        }
431
432 2
        $zones = 'CDEFGHJKLMNPQRSTUVWXX';
433 2
        $zoneIndex = (int) (($latitude + 80) / 8);
434
435 2
        return $zones[$zoneIndex];
436
    }
437
438
    /**
439
     * Convert a latitude and longitude to easting and northing using a Transverse Mercator projection
440
     * Formula for transformation is taken from OS document
441
     * "A Guide to Coordinate Systems in Great Britain".
442
     *
443
     * @param float $scale          scale factor on central meridian
444
     * @param float $originEasting  easting of true origin
445
     * @param float $originNorthing northing of true origin
446
     * @param float $originLat      latitude of true origin
447
     * @param float $originLong     longitude of true origin
448
     *
449
     * @return array
450
     */
451 5
    public function toTransverseMercatorEastingNorthing(
452
        $scale,
453
        $originEasting,
454
        $originNorthing,
455
        $originLat,
456
        $originLong
457
    ): array {
458 5
        $originLat = deg2rad($originLat);
459 5
        $originLong = deg2rad($originLong);
460
461 5
        $lat = deg2rad($this->lat);
462 5
        $sinLat = sin($lat);
463 5
        $cosLat = cos($lat);
464 5
        $tanLat = tan($lat);
465 5
        $tanLatSq = $tanLat ** 2;
466 5
        $long = deg2rad($this->lng);
467
468 5
        $n = ($this->refEll->getMaj() - $this->refEll->getMin()) / ($this->refEll->getMaj() + $this->refEll->getMin());
469 5
        $nSq = $n ** 2;
470 5
        $nCu = $n ** 3;
471
472 5
        $v = $this->refEll->getMaj() * $scale * ((1 - $this->refEll->getEcc() * ($sinLat ** 2)) ** -0.5);
473 5
        $p = $this->refEll->getMaj() * $scale * (1 - $this->refEll->getEcc()) * ((1 - $this->refEll->getEcc() * ($sinLat ** 2)) ** -1.5);
474 5
        $hSq = (($v / $p) - 1);
475
476 5
        $latPlusOrigin = $lat + $originLat;
477 5
        $latMinusOrigin = $lat - $originLat;
478
479 5
        $longMinusOrigin = $long - $originLong;
480
481 5
        $M = $this->refEll->getMin() * $scale
482 5
            * ((1 + $n + 1.25 * ($nSq + $nCu)) * $latMinusOrigin
483 5
                - (3 * ($n + $nSq) + 2.625 * $nCu) * sin($latMinusOrigin) * cos($latPlusOrigin)
484 5
                + 1.875 * ($nSq + $nCu) * sin(2 * $latMinusOrigin) * cos(2 * $latPlusOrigin)
485 5
                - (35 / 24 * $nCu * sin(3 * $latMinusOrigin) * cos(3 * $latPlusOrigin)));
486
487 5
        $I = $M + $originNorthing;
488 5
        $II = $v / 2 * $sinLat * $cosLat;
489 5
        $III = $v / 24 * $sinLat * ($cosLat ** 3) * (5 - $tanLatSq + 9 * $hSq);
490 5
        $IIIA = $v / 720 * $sinLat * ($cosLat ** 5) * (61 - 58 * $tanLatSq + ($tanLatSq ** 2));
491 5
        $IV = $v * $cosLat;
492 5
        $V = $v / 6 * ($cosLat ** 3) * ($v / $p - $tanLatSq);
493 5
        $VI = $v / 120 * ($cosLat ** 5) * (5 - 18 * $tanLatSq + ($tanLatSq ** 2) + 14 * $hSq - 58 * $tanLatSq * $hSq);
494
495 5
        $E = (int) round($originEasting + $IV * $longMinusOrigin + $V * ($longMinusOrigin ** 3) + $VI * ($longMinusOrigin ** 5));
496 5
        $N = (int) round($I + $II * ($longMinusOrigin ** 2) + $III * ($longMinusOrigin ** 4) + $IIIA * ($longMinusOrigin ** 6));
497
498 5
        return array('E' => $E, 'N' => $N);
499
    }
500
}
501