Passed
Push — master ( fd93b5...48e916 )
by Doug
40:26 queued 29:39
created

GeographicPoint::coordinateFrameRotation()   A

Complexity

Conditions 1
Paths 1

Size

Total Lines 22
Code Lines 12

Duplication

Lines 0
Ratio 0 %

Code Coverage

Tests 6
CRAP Score 1

Importance

Changes 0
Metric Value
eloc 12
c 0
b 0
f 0
dl 0
loc 22
ccs 6
cts 6
cp 1
rs 9.8666
cc 1
nc 1
nop 8
crap 1

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
 * PHPCoord.
4
 *
5
 * @author Doug Wright
6
 */
7
declare(strict_types=1);
8
9
namespace PHPCoord;
10
11
use function abs;
12
use function asinh;
13
use function atan;
14
use function atan2;
15
use function atanh;
16
use function cos;
17
use function cosh;
18
use DateTime;
19
use DateTimeImmutable;
20
use DateTimeInterface;
21
use function get_class;
22
use function hypot;
23
use function implode;
24
use function is_nan;
25
use function log;
26
use const M_E;
27
use const M_PI;
28
use function max;
29
use PHPCoord\CoordinateOperation\AutoConversion;
30
use PHPCoord\CoordinateOperation\ComplexNumber;
31
use PHPCoord\CoordinateOperation\ConvertiblePoint;
32
use PHPCoord\CoordinateOperation\GeocentricValue;
33
use PHPCoord\CoordinateOperation\GeographicGeoidHeightGrid;
34
use PHPCoord\CoordinateOperation\GeographicGrid;
35
use PHPCoord\CoordinateOperation\GeographicValue;
36
use PHPCoord\CoordinateOperation\NADCON5Grid;
37
use PHPCoord\CoordinateOperation\NADCON5Grids;
38
use PHPCoord\CoordinateOperation\OSTNOSGM15Grid;
39
use PHPCoord\CoordinateReferenceSystem\Compound;
40
use PHPCoord\CoordinateReferenceSystem\Geocentric;
41
use PHPCoord\CoordinateReferenceSystem\Geographic;
42
use PHPCoord\CoordinateReferenceSystem\Geographic2D;
43
use PHPCoord\CoordinateReferenceSystem\Geographic3D;
44
use PHPCoord\CoordinateReferenceSystem\Projected;
0 ignored issues
show
Bug introduced by
The type PHPCoord\CoordinateReferenceSystem\Projected was not found. Maybe you did not declare it correctly or list all dependencies?

The issue could also be caused by a filter entry in the build configuration. If the path has been excluded in your configuration, e.g. excluded_paths: ["lib/*"], you can move it to the dependency path list as follows:

filter:
    dependency_paths: ["lib/*"]

For further information see https://scrutinizer-ci.com/docs/tools/php/php-scrutinizer/#list-dependency-paths

Loading history...
45
use PHPCoord\CoordinateReferenceSystem\Vertical;
46
use PHPCoord\CoordinateSystem\Axis;
47
use PHPCoord\CoordinateSystem\Cartesian;
48
use PHPCoord\Datum\Datum;
49
use PHPCoord\Datum\Ellipsoid;
50
use PHPCoord\Exception\InvalidCoordinateReferenceSystemException;
51
use PHPCoord\Exception\UnknownAxisException;
52
use PHPCoord\Geometry\BoundingArea;
53
use PHPCoord\UnitOfMeasure\Angle\Angle;
54
use PHPCoord\UnitOfMeasure\Angle\ArcSecond;
55
use PHPCoord\UnitOfMeasure\Angle\Degree;
56
use PHPCoord\UnitOfMeasure\Angle\Radian;
57
use PHPCoord\UnitOfMeasure\Length\Length;
58
use PHPCoord\UnitOfMeasure\Length\Metre;
59
use PHPCoord\UnitOfMeasure\Scale\Coefficient;
60
use PHPCoord\UnitOfMeasure\Scale\Scale;
61
use PHPCoord\UnitOfMeasure\Scale\Unity;
62
use function sin;
63
use function sinh;
64
use function sprintf;
65
use function sqrt;
66
use function tan;
67
use TypeError;
68
69
/**
70
 * Coordinate representing a point on an ellipsoid.
71
 */
72
class GeographicPoint extends Point implements ConvertiblePoint
73
{
74
    use AutoConversion;
75
76
    /**
77
     * Latitude.
78
     */
79
    protected Angle $latitude;
80
81
    /**
82
     * Longitude.
83
     */
84
    protected Angle $longitude;
85
86
    /**
87
     * Height above ellipsoid (N.B. *not* height above ground, sea-level or anything else tangible).
88
     */
89
    protected ?Length $height;
90
91
    /**
92
     * Coordinate reference system.
93
     */
94
    protected Geographic $crs;
95
96
    /**
97
     * Coordinate epoch (date for which the specified coordinates represented this point).
98
     */
99
    protected ?DateTimeImmutable $epoch;
100
101 1499
    protected function __construct(Angle $latitude, Angle $longitude, ?Length $height, Geographic $crs, ?DateTimeInterface $epoch = null)
102
    {
103 1499
        if (!$crs instanceof Geographic2D && !$crs instanceof Geographic3D) {
104
            throw new TypeError(sprintf("A geographic point must be associated with a geographic CRS, but a '%s' CRS was given", get_class($crs)));
105
        }
106
107 1499
        if ($crs instanceof Geographic2D && $height !== null) {
108 9
            throw new InvalidCoordinateReferenceSystemException('A 2D geographic point must not include a height');
109
        }
110
111 1490
        if ($crs instanceof Geographic3D && $height === null) {
112 9
            throw new InvalidCoordinateReferenceSystemException('A 3D geographic point must include a height, none given');
113
        }
114
115 1481
        $this->crs = $crs;
116
117 1481
        $latitude = $this->normaliseLatitude($latitude);
118 1481
        $longitude = $this->normaliseLongitude($longitude);
119
120 1481
        $this->latitude = $latitude::convert($latitude, $this->crs->getCoordinateSystem()->getAxisByName(Axis::GEODETIC_LATITUDE)->getUnitOfMeasureId());
121 1481
        $this->longitude = $longitude::convert($longitude, $this->crs->getCoordinateSystem()->getAxisByName(Axis::GEODETIC_LONGITUDE)->getUnitOfMeasureId());
122
123 1481
        if ($height) {
124 191
            $this->height = $height::convert($height, $this->crs->getCoordinateSystem()->getAxisByName(Axis::ELLIPSOIDAL_HEIGHT)->getUnitOfMeasureId());
125
        } else {
126 1369
            $this->height = null;
127
        }
128
129 1481
        if ($epoch instanceof DateTime) {
130 9
            $epoch = DateTimeImmutable::createFromMutable($epoch);
131
        }
132 1481
        $this->epoch = $epoch;
133 1481
    }
134
135
    /**
136
     * @param Angle   $latitude  refer to CRS for preferred unit of measure, but any angle unit accepted
137
     * @param Angle   $longitude refer to CRS for preferred unit of measure, but any angle unit accepted
138
     * @param ?Length $height    refer to CRS for preferred unit of measure, but any length unit accepted
139
     */
140 1499
    public static function create(Angle $latitude, Angle $longitude, ?Length $height = null, Geographic $crs, ?DateTimeInterface $epoch = null): self
141
    {
142 1499
        return new static($latitude, $longitude, $height, $crs, $epoch);
143
    }
144
145 804
    public function getLatitude(): Angle
146
    {
147 804
        return $this->latitude;
148
    }
149
150 768
    public function getLongitude(): Angle
151
    {
152 768
        return $this->longitude;
153
    }
154
155 348
    public function getHeight(): ?Length
156
    {
157 348
        return $this->height;
158
    }
159
160 398
    public function getCRS(): Geographic
161
    {
162 398
        return $this->crs;
163
    }
164
165 45
    public function getCoordinateEpoch(): ?DateTimeImmutable
166
    {
167 45
        return $this->epoch;
168
    }
169
170 1481
    protected function normaliseLatitude(Angle $latitude): Angle
171
    {
172 1481
        if ($latitude->asDegrees()->getValue() > 90) {
173
            return new Degree(90);
174
        }
175 1481
        if ($latitude->asDegrees()->getValue() < -90) {
176
            return new Degree(-90);
177
        }
178
179 1481
        return $latitude;
180
    }
181
182 1481
    protected function normaliseLongitude(Angle $longitude): Angle
183
    {
184 1481
        while ($longitude->asDegrees()->getValue() > 180) {
185 9
            $longitude = $longitude->subtract(new Degree(360));
186
        }
187 1481
        while ($longitude->asDegrees()->getValue() <= -180) {
188
            $longitude = $longitude->add(new Degree(360));
189
        }
190
191 1481
        return $longitude;
192
    }
193
194
    /**
195
     * Calculate surface distance between two points.
196
     */
197 153
    public function calculateDistance(Point $to): Length
198
    {
199
        try {
200 153
            if ($to instanceof ConvertiblePoint) {
201 144
                $to = $to->convert($this->crs);
202
            }
203
        } finally {
204 153
            if ($to->getCRS()->getSRID() !== $this->crs->getSRID()) {
205 9
                throw new InvalidCoordinateReferenceSystemException('Can only calculate distances between two points in the same CRS');
206
            }
207
208
            /* @var GeographicPoint $to */
209 144
            return static::vincenty($this->asGeographicValue(), $to->asGeographicValue(), $this->getCRS()->getDatum()->getEllipsoid());
210
        }
211
    }
212
213 36
    public function __toString(): string
214
    {
215 36
        $values = [];
216 36
        foreach ($this->getCRS()->getCoordinateSystem()->getAxes() as $axis) {
217 36
            if ($axis->getName() === Axis::GEODETIC_LATITUDE) {
218 36
                $values[] = $this->latitude;
219 36
            } elseif ($axis->getName() === Axis::GEODETIC_LONGITUDE) {
220 36
                $values[] = $this->longitude;
221 9
            } elseif ($axis->getName() === Axis::ELLIPSOIDAL_HEIGHT) {
222 9
                $values[] = $this->height;
223
            } else {
224
                throw new UnknownAxisException(); // @codeCoverageIgnore
225
            }
226
        }
227
228 36
        return '(' . implode(', ', $values) . ')';
229
    }
230
231
    /**
232
     * Geographic/geocentric conversions
233
     * In applications it is often concatenated with the 3- 7- or 10-parameter transformations 9603, 9606, 9607 or
234
     * 9636 to form a geographic to geographic transformation.
235
     */
236 9
    public function geographicGeocentric(
237
        Geocentric $to
238
    ): GeocentricPoint {
239 9
        $geographicValue = new GeographicValue($this->latitude, $this->longitude, $this->height, $this->crs->getDatum());
240 9
        $asGeocentric = $geographicValue->asGeocentricValue();
241
242 9
        return GeocentricPoint::create($asGeocentric->getX(), $asGeocentric->getY(), $asGeocentric->getZ(), $to, $this->epoch);
243
    }
244
245
    /**
246
     * Coordinate Frame rotation (geog2D/geog3D domain)
247
     * Note the analogy with the Position Vector tfm (codes 9606/1037) but beware of the differences!  The Position Vector
248
     * convention is used by IAG and recommended by ISO 19111. See methods 1032/1038/9607 for similar tfms operating
249
     * between other CRS types.
250
     */
251 72
    public function coordinateFrameRotation(
252
        Geographic $to,
253
        Length $xAxisTranslation,
254
        Length $yAxisTranslation,
255
        Length $zAxisTranslation,
256
        Angle $xAxisRotation,
257
        Angle $yAxisRotation,
258
        Angle $zAxisRotation,
259
        Scale $scaleDifference
260
    ): self {
261 72
        return $this->coordinateFrameMolodenskyBadekas(
262 72
            $to,
263
            $xAxisTranslation,
264
            $yAxisTranslation,
265
            $zAxisTranslation,
266
            $xAxisRotation,
267
            $yAxisRotation,
268
            $zAxisRotation,
269
            $scaleDifference,
270 72
            new Metre(0),
271 72
            new Metre(0),
272 72
            new Metre(0)
273
        );
274
    }
275
276
    /**
277
     * Molodensky-Badekas (CF geog2D/geog3D domain)
278
     * See method codes 1034 and 1039/9636 for this operation in other coordinate domains and method code 1062/1063 for the
279
     * opposite rotation convention in geographic 2D domain.
280
     */
281 90
    public function coordinateFrameMolodenskyBadekas(
282
        Geographic $to,
283
        Length $xAxisTranslation,
284
        Length $yAxisTranslation,
285
        Length $zAxisTranslation,
286
        Angle $xAxisRotation,
287
        Angle $yAxisRotation,
288
        Angle $zAxisRotation,
289
        Scale $scaleDifference,
290
        Length $ordinate1OfEvaluationPoint,
291
        Length $ordinate2OfEvaluationPoint,
292
        Length $ordinate3OfEvaluationPoint
293
    ): self {
294 90
        $geographicValue = new GeographicValue($this->latitude, $this->longitude, $this->height, $this->crs->getDatum());
295 90
        $asGeocentric = $geographicValue->asGeocentricValue();
296
297 90
        $xs = $asGeocentric->getX()->asMetres()->getValue();
298 90
        $ys = $asGeocentric->getY()->asMetres()->getValue();
299 90
        $zs = $asGeocentric->getZ()->asMetres()->getValue();
300 90
        $tx = $xAxisTranslation->asMetres()->getValue();
301 90
        $ty = $yAxisTranslation->asMetres()->getValue();
302 90
        $tz = $zAxisTranslation->asMetres()->getValue();
303 90
        $rx = $xAxisRotation->asRadians()->getValue();
304 90
        $ry = $yAxisRotation->asRadians()->getValue();
305 90
        $rz = $zAxisRotation->asRadians()->getValue();
306 90
        $M = 1 + $scaleDifference->asUnity()->getValue();
307 90
        $xp = $ordinate1OfEvaluationPoint->asMetres()->getValue();
308 90
        $yp = $ordinate2OfEvaluationPoint->asMetres()->getValue();
309 90
        $zp = $ordinate3OfEvaluationPoint->asMetres()->getValue();
310
311 90
        $xt = $M * ((($xs - $xp) * 1) + (($ys - $yp) * $rz) + (($zs - $zp) * -$ry)) + $tx + $xp;
312 90
        $yt = $M * ((($xs - $xp) * -$rz) + (($ys - $yp) * 1) + (($zs - $zp) * $rx)) + $ty + $yp;
313 90
        $zt = $M * ((($xs - $xp) * $ry) + (($ys - $yp) * -$rx) + (($zs - $zp) * 1)) + $tz + $zp;
314 90
        $newGeocentric = new GeocentricValue(new Metre($xt), new Metre($yt), new Metre($zt), $to->getDatum());
315 90
        $newGeographic = $newGeocentric->asGeographicValue();
316
317 90
        return static::create($newGeographic->getLatitude(), $newGeographic->getLongitude(), $to instanceof Geographic3D ? $newGeographic->getHeight() : null, $to, $this->epoch);
318
    }
319
320
    /**
321
     * Position Vector transformation (geog2D/geog3D domain)
322
     * Note the analogy with the Coordinate Frame rotation (code 9607/1038) but beware of the differences!  The Position
323
     * Vector convention is used by IAG and recommended by ISO 19111. See methods 1033/1037/9606 for similar tfms
324
     * operating between other CRS types.
325
     */
326 180
    public function positionVectorTransformation(
327
        Geographic $to,
328
        Length $xAxisTranslation,
329
        Length $yAxisTranslation,
330
        Length $zAxisTranslation,
331
        Angle $xAxisRotation,
332
        Angle $yAxisRotation,
333
        Angle $zAxisRotation,
334
        Scale $scaleDifference
335
    ): self {
336 180
        return $this->positionVectorMolodenskyBadekas(
337 180
            $to,
338
            $xAxisTranslation,
339
            $yAxisTranslation,
340
            $zAxisTranslation,
341
            $xAxisRotation,
342
            $yAxisRotation,
343
            $zAxisRotation,
344
            $scaleDifference,
345 180
            new Metre(0),
346 180
            new Metre(0),
347 180
            new Metre(0)
348
        );
349
    }
350
351
    /**
352
     * Molodensky-Badekas (PV geog2D/geog3D domain)
353
     * See method codes 1061 and 1062/1063 for this operation in other coordinate domains and method code 1039/9636 for opposite
354
     * rotation in geographic 2D/3D domain.
355
     */
356 198
    public function positionVectorMolodenskyBadekas(
357
        Geographic $to,
358
        Length $xAxisTranslation,
359
        Length $yAxisTranslation,
360
        Length $zAxisTranslation,
361
        Angle $xAxisRotation,
362
        Angle $yAxisRotation,
363
        Angle $zAxisRotation,
364
        Scale $scaleDifference,
365
        Length $ordinate1OfEvaluationPoint,
366
        Length $ordinate2OfEvaluationPoint,
367
        Length $ordinate3OfEvaluationPoint
368
    ): self {
369 198
        $geographicValue = new GeographicValue($this->latitude, $this->longitude, $this->height, $this->crs->getDatum());
370 198
        $asGeocentric = $geographicValue->asGeocentricValue();
371
372 198
        $xs = $asGeocentric->getX()->asMetres()->getValue();
373 198
        $ys = $asGeocentric->getY()->asMetres()->getValue();
374 198
        $zs = $asGeocentric->getZ()->asMetres()->getValue();
375 198
        $tx = $xAxisTranslation->asMetres()->getValue();
376 198
        $ty = $yAxisTranslation->asMetres()->getValue();
377 198
        $tz = $zAxisTranslation->asMetres()->getValue();
378 198
        $rx = $xAxisRotation->asRadians()->getValue();
379 198
        $ry = $yAxisRotation->asRadians()->getValue();
380 198
        $rz = $zAxisRotation->asRadians()->getValue();
381 198
        $M = 1 + $scaleDifference->asUnity()->getValue();
382 198
        $xp = $ordinate1OfEvaluationPoint->asMetres()->getValue();
383 198
        $yp = $ordinate2OfEvaluationPoint->asMetres()->getValue();
384 198
        $zp = $ordinate3OfEvaluationPoint->asMetres()->getValue();
385
386 198
        $xt = $M * ((($xs - $xp) * 1) + (($ys - $yp) * -$rz) + (($zs - $zp) * $ry)) + $tx + $xp;
387 198
        $yt = $M * ((($xs - $xp) * $rz) + (($ys - $yp) * 1) + (($zs - $zp) * -$rx)) + $ty + $yp;
388 198
        $zt = $M * ((($xs - $xp) * -$ry) + (($ys - $yp) * $rx) + (($zs - $zp) * 1)) + $tz + $zp;
389 198
        $newGeocentric = new GeocentricValue(new Metre($xt), new Metre($yt), new Metre($zt), $to->getDatum());
390 198
        $newGeographic = $newGeocentric->asGeographicValue();
391
392 198
        return static::create($newGeographic->getLatitude(), $newGeographic->getLongitude(), $to instanceof Geographic3D ? $newGeographic->getHeight() : null, $to, $this->epoch);
393
    }
394
395
    /**
396
     * Geocentric translations
397
     * This method allows calculation of geocentric coords in the target system by adding the parameter values to the
398
     * corresponding coordinates of the point in the source system. See methods 1031 and 1035 for similar tfms
399
     * operating between other CRSs types.
400
     */
401 91
    public function geocentricTranslation(
402
        Geographic $to,
403
        Length $xAxisTranslation,
404
        Length $yAxisTranslation,
405
        Length $zAxisTranslation
406
    ): self {
407 91
        return $this->positionVectorTransformation(
408 91
            $to,
409
            $xAxisTranslation,
410
            $yAxisTranslation,
411
            $zAxisTranslation,
412 91
            new Radian(0),
413 91
            new Radian(0),
414 91
            new Radian(0),
415 91
            new Unity(0)
416
        );
417
    }
418
419
    /**
420
     * Abridged Molodensky
421
     * This transformation is a truncated Taylor series expansion of a transformation between two geographic coordinate
422
     * systems, modelled as a set of geocentric translations.
423
     */
424 18
    public function abridgedMolodensky(
425
        Geographic $to,
426
        Length $xAxisTranslation,
427
        Length $yAxisTranslation,
428
        Length $zAxisTranslation,
429
        Length $differenceInSemiMajorAxis,
430
        Scale $differenceInFlattening
431
    ): self {
432 18
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
433 18
        $latitude = $this->latitude->asRadians()->getValue();
434 18
        $longitude = $this->longitude->asRadians()->getValue();
435 18
        $fromHeight = $this->height ? $this->height->asMetres()->getValue() : 0;
436 18
        $tx = $xAxisTranslation->asMetres()->getValue();
437 18
        $ty = $yAxisTranslation->asMetres()->getValue();
438 18
        $tz = $zAxisTranslation->asMetres()->getValue();
439 18
        $da = $differenceInSemiMajorAxis->asMetres()->getValue();
440 18
        $df = $differenceInFlattening->asUnity()->getValue();
441
442 18
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
443 18
        $e2 = $ellipsoid->getEccentricitySquared();
444
445 18
        $rho = $a * (1 - $e2) / (1 - $e2 * sin($latitude) ** 2) ** (3 / 2);
446 18
        $nu = $a / sqrt(1 - $e2 * (sin($latitude) ** 2));
447
448 18
        $f = $ellipsoid->getInverseFlattening();
449
450 18
        $dLatitude = ((-$tx * sin($latitude) * cos($longitude)) - ($ty * sin($latitude) * sin($longitude)) + ($tz * cos($latitude)) + ((($a * $df) + ($ellipsoid->getInverseFlattening() * $da)) * sin(2 * $latitude))) / ($rho * sin((new ArcSecond(1))->asRadians()->getValue()));
451 18
        $dLongitude = (-$tx * sin($longitude) + $ty * cos($longitude)) / (($nu * cos($latitude)) * sin((new ArcSecond(1))->asRadians()->getValue()));
452 18
        $dHeight = ($tx * cos($latitude) * cos($longitude)) + ($ty * cos($latitude) * sin($longitude)) + ($tz * sin($latitude)) + (($a * $df + $f * $da) * (sin($latitude) ** 2)) - $da;
453
454 18
        $toLatitude = $latitude + (new ArcSecond($dLatitude))->asRadians()->getValue();
455 18
        $toLongitude = $longitude + (new ArcSecond($dLongitude))->asRadians()->getValue();
456 18
        $toHeight = $fromHeight + $dHeight;
457
458 18
        return static::create(new Radian($toLatitude), new Radian($toLongitude), $to instanceof Geographic3D ? new Metre($toHeight) : null, $to, $this->epoch);
459
    }
460
461
    /**
462
     * Molodensky
463
     * See Abridged Molodensky.
464
     */
465 18
    public function molodensky(
466
        Geographic $to,
467
        Length $xAxisTranslation,
468
        Length $yAxisTranslation,
469
        Length $zAxisTranslation,
470
        Length $differenceInSemiMajorAxis,
471
        Scale $differenceInFlattening
472
    ): self {
473 18
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
474 18
        $latitude = $this->latitude->asRadians()->getValue();
475 18
        $longitude = $this->longitude->asRadians()->getValue();
476 18
        $fromHeight = $this->height ? $this->height->asMetres()->getValue() : 0;
477 18
        $tx = $xAxisTranslation->asMetres()->getValue();
478 18
        $ty = $yAxisTranslation->asMetres()->getValue();
479 18
        $tz = $zAxisTranslation->asMetres()->getValue();
480 18
        $da = $differenceInSemiMajorAxis->asMetres()->getValue();
481 18
        $df = $differenceInFlattening->asUnity()->getValue();
482
483 18
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
484 18
        $b = $ellipsoid->getSemiMinorAxis()->asMetres()->getValue();
485 18
        $e2 = $ellipsoid->getEccentricitySquared();
486
487 18
        $rho = $a * (1 - $e2) / (1 - $e2 * sin($latitude) ** 2) ** (3 / 2);
488 18
        $nu = $a / sqrt(1 - $e2 * (sin($latitude) ** 2));
489
490 18
        $f = $ellipsoid->getInverseFlattening();
0 ignored issues
show
Unused Code introduced by
The assignment to $f is dead and can be removed.
Loading history...
491
492 18
        $dLatitude = ((-$tx * sin($latitude) * cos($longitude)) - ($ty * sin($latitude) * sin($longitude)) + ($tz * cos($latitude)) + ($da * ($nu * $e2 * sin($latitude) * cos($latitude)) / $a + $df * ($rho * ($a / $b) + $nu * ($b / $a)) * sin($latitude) * cos($latitude))) / (($rho + $fromHeight) * sin((new ArcSecond(1))->asRadians()->getValue()));
493 18
        $dLongitude = (-$tx * sin($longitude) + $ty * cos($longitude)) / ((($nu + $fromHeight) * cos($latitude)) * sin((new ArcSecond(1))->asRadians()->getValue()));
494 18
        $dHeight = ($tx * cos($latitude) * cos($longitude)) + ($ty * cos($latitude) * sin($longitude)) + ($tz * sin($latitude)) - $da * $a / $nu + $df * $b / $a * $nu * sin($latitude) ** 2;
495
496 18
        $toLatitude = $latitude + (new ArcSecond($dLatitude))->asRadians()->getValue();
497 18
        $toLongitude = $longitude + (new ArcSecond($dLongitude))->asRadians()->getValue();
498 18
        $toHeight = $fromHeight + $dHeight;
499
500 18
        return static::create(new Radian($toLatitude), new Radian($toLongitude), $to instanceof Geographic3D ? new Metre($toHeight) : null, $to, $this->epoch);
501
    }
502
503
    /**
504
     * Albers Equal Area.
505
     */
506 18
    public function albersEqualArea(
507
        Projected $to,
508
        Angle $latitudeOfFalseOrigin,
509
        Angle $longitudeOfFalseOrigin,
510
        Angle $latitudeOf1stStandardParallel,
511
        Angle $latitudeOf2ndStandardParallel,
512
        Length $eastingAtFalseOrigin,
513
        Length $northingAtFalseOrigin
514
    ): ProjectedPoint {
515 18
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
516 18
        $latitude = $this->latitude->asRadians()->getValue();
517 18
        $longitude = $this->longitude->asRadians()->getValue();
0 ignored issues
show
Unused Code introduced by
The assignment to $longitude is dead and can be removed.
Loading history...
518 18
        $phiOrigin = $latitudeOfFalseOrigin->asRadians()->getValue();
519 18
        $phi1 = $latitudeOf1stStandardParallel->asRadians()->getValue();
520 18
        $phi2 = $latitudeOf2ndStandardParallel->asRadians()->getValue();
521 18
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
522 18
        $e = $ellipsoid->getEccentricity();
523 18
        $e2 = $ellipsoid->getEccentricitySquared();
524
525 18
        $centralMeridianFirstParallel = cos($phi1) / sqrt(1 - ($e2 * sin($phi1) ** 2));
526 18
        $centralMeridianSecondParallel = cos($phi2) / sqrt(1 - ($e2 * sin($phi2) ** 2));
527
528 18
        $alpha = (1 - $e2) * (sin($latitude) / (1 - $e2 * sin($latitude) ** 2) - (1 / 2 / $e) * log((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))));
529 18
        $alphaOrigin = (1 - $e2) * (sin($phiOrigin) / (1 - $e2 * sin($phiOrigin) ** 2) - (1 / 2 / $e) * log((1 - $e * sin($phiOrigin)) / (1 + $e * sin($phiOrigin))));
530 18
        $alphaFirstParallel = (1 - $e2) * (sin($phi1) / (1 - $e2 * sin($phi1) ** 2) - (1 / 2 / $e) * log((1 - $e * sin($phi1)) / (1 + $e * sin($phi1))));
531 18
        $alphaSecondParallel = (1 - $e2) * (sin($phi2) / (1 - $e2 * sin($phi2) ** 2) - (1 / 2 / $e) * log((1 - $e * sin($phi2)) / (1 + $e * sin($phi2))));
532
533 18
        $n = ($centralMeridianFirstParallel ** 2 - $centralMeridianSecondParallel ** 2) / ($alphaSecondParallel - $alphaFirstParallel);
534 18
        $C = $centralMeridianFirstParallel ** 2 + $n * $alphaFirstParallel;
535 18
        $theta = $n * $this->normaliseLongitude($this->longitude->subtract($longitudeOfFalseOrigin))->asRadians()->getValue();
536 18
        $rho = $a * sqrt($C - $n * $alpha) / $n;
537 18
        $rhoOrigin = ($a * sqrt($C - $n * $alphaOrigin)) / $n;
538
539 18
        $easting = $eastingAtFalseOrigin->asMetres()->getValue() + ($rho * sin($theta));
540 18
        $northing = $northingAtFalseOrigin->asMetres()->getValue() + $rhoOrigin - ($rho * cos($theta));
541
542 18
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
543
    }
544
545
    /**
546
     * American Polyconic.
547
     */
548 9
    public function americanPolyconic(
549
        Projected $to,
550
        Angle $latitudeOfNaturalOrigin,
551
        Angle $longitudeOfNaturalOrigin,
552
        Length $falseEasting,
553
        Length $falseNorthing
554
    ): ProjectedPoint {
555 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
556 9
        $latitude = $this->latitude->asRadians()->getValue();
557 9
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
558 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
559 9
        $e = $ellipsoid->getEccentricity();
560 9
        $e2 = $ellipsoid->getEccentricitySquared();
561 9
        $e4 = $e ** 4;
562 9
        $e6 = $e ** 6;
563
564 9
        $M = $a * ((1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256) * $latitude - (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024) * sin(2 * $latitude) + (15 * $e4 / 256 + 45 * $e6 / 1024) * sin(4 * $latitude) - (35 * $e6 / 3072) * sin(6 * $latitude));
565 9
        $MO = $a * ((1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256) * $latitudeOrigin - (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024) * sin(2 * $latitudeOrigin) + (15 * $e4 / 256 + 45 * $e6 / 1024) * sin(4 * $latitudeOrigin) - (35 * $e6 / 3072) * sin(6 * $latitudeOrigin));
566
567 9
        if ($latitude === 0.0) {
0 ignored issues
show
introduced by
The condition $latitude === 0.0 is always false.
Loading history...
568
            $easting = $falseEasting->asMetres()->getValue() + $a * $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue();
569
            $northing = $falseNorthing->asMetres()->getValue() - $MO;
570
        } else {
571 9
            $L = $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue() * sin($latitude);
572 9
            $nu = $a / sqrt(1 - $e2 * sin($latitude) ** 2);
573
574 9
            $easting = $falseEasting->asMetres()->getValue() + $nu * 1 / tan($latitude) * sin($L);
575 9
            $northing = $falseNorthing->asMetres()->getValue() + $M - $MO + $nu * 1 / tan($latitude) * (1 - cos($L));
576
        }
577
578 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
579
    }
580
581
    /**
582
     * Bonne.
583
     */
584 9
    public function bonne(
585
        Projected $to,
586
        Angle $latitudeOfNaturalOrigin,
587
        Angle $longitudeOfNaturalOrigin,
588
        Length $falseEasting,
589
        Length $falseNorthing
590
    ): ProjectedPoint {
591 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
592 9
        $latitude = $this->latitude->asRadians()->getValue();
593 9
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
594 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
595 9
        $e = $ellipsoid->getEccentricity();
596 9
        $e2 = $ellipsoid->getEccentricitySquared();
597 9
        $e4 = $e ** 4;
598 9
        $e6 = $e ** 6;
599
600 9
        $m = cos($latitude) / sqrt(1 - $e2 * sin($latitude) ** 2);
601 9
        $mO = cos($latitudeOrigin) / sqrt(1 - $e2 * sin($latitudeOrigin) ** 2);
602
603 9
        $M = $a * ((1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256) * $latitude - (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024) * sin(2 * $latitude) + (15 * $e4 / 256 + 45 * $e6 / 1024) * sin(4 * $latitude) - (35 * $e6 / 3072) * sin(6 * $latitude));
604 9
        $MO = $a * ((1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256) * $latitudeOrigin - (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024) * sin(2 * $latitudeOrigin) + (15 * $e4 / 256 + 45 * $e6 / 1024) * sin(4 * $latitudeOrigin) - (35 * $e6 / 3072) * sin(6 * $latitudeOrigin));
605
606 9
        $rho = $a * $mO / sin($latitudeOrigin) + $MO - $M;
607 9
        $tau = $a * $m * $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue() / $rho;
608
609 9
        $easting = $falseEasting->asMetres()->getValue() + ($rho * sin($tau));
610 9
        $northing = $falseNorthing->asMetres()->getValue() + (($a * $mO / sin($latitudeOrigin) - $rho * cos($tau)));
611
612 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
613
    }
614
615
    /**
616
     * Bonne South Orientated.
617
     */
618 9
    public function bonneSouthOrientated(
619
        Projected $to,
620
        Angle $latitudeOfNaturalOrigin,
621
        Angle $longitudeOfNaturalOrigin,
622
        Length $falseEasting,
623
        Length $falseNorthing
624
    ): ProjectedPoint {
625 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
626 9
        $latitude = $this->latitude->asRadians()->getValue();
627 9
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
628 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
629 9
        $e = $ellipsoid->getEccentricity();
630 9
        $e2 = $ellipsoid->getEccentricitySquared();
631 9
        $e4 = $e ** 4;
632 9
        $e6 = $e ** 6;
633
634 9
        $m = cos($latitude) / sqrt(1 - $e2 * sin($latitude) ** 2);
635 9
        $mO = cos($latitudeOrigin) / sqrt(1 - $e2 * sin($latitudeOrigin) ** 2);
636
637 9
        $M = $a * ((1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256) * $latitude - (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024) * sin(2 * $latitude) + (15 * $e4 / 256 + 45 * $e6 / 1024) * sin(4 * $latitude) - (35 * $e6 / 3072) * sin(6 * $latitude));
638 9
        $MO = $a * ((1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256) * $latitudeOrigin - (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024) * sin(2 * $latitudeOrigin) + (15 * $e4 / 256 + 45 * $e6 / 1024) * sin(4 * $latitudeOrigin) - (35 * $e6 / 3072) * sin(6 * $latitudeOrigin));
639
640 9
        $rho = $a * $mO / sin($latitudeOrigin) + $MO - $M;
641 9
        $tau = $a * $m * $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue() / $rho;
642
643 9
        $westing = $falseEasting->asMetres()->getValue() - ($rho * sin($tau));
644 9
        $southing = $falseNorthing->asMetres()->getValue() - (($a * $mO / sin($latitudeOrigin) - $rho * cos($tau)));
645
646 9
        return ProjectedPoint::create(new Metre(-$westing), new Metre(-$southing), new Metre($westing), new Metre($southing), $to, $this->epoch);
647
    }
648
649
    /**
650
     * Cassini-Soldner.
651
     */
652 9
    public function cassiniSoldner(
653
        Projected $to,
654
        Angle $latitudeOfNaturalOrigin,
655
        Angle $longitudeOfNaturalOrigin,
656
        Length $falseEasting,
657
        Length $falseNorthing
658
    ): ProjectedPoint {
659 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
660 9
        $latitude = $this->latitude->asRadians()->getValue();
661 9
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
662 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
663 9
        $e = $ellipsoid->getEccentricity();
664 9
        $e2 = $ellipsoid->getEccentricitySquared();
665 9
        $e4 = $e ** 4;
666 9
        $e6 = $e ** 6;
667
668 9
        $M = $a * ((1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256) * $latitude - (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024) * sin(2 * $latitude) + (15 * $e4 / 256 + 45 * $e6 / 1024) * sin(4 * $latitude) - (35 * $e6 / 3072) * sin(6 * $latitude));
669 9
        $MO = $a * ((1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256) * $latitudeOrigin - (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024) * sin(2 * $latitudeOrigin) + (15 * $e4 / 256 + 45 * $e6 / 1024) * sin(4 * $latitudeOrigin) - (35 * $e6 / 3072) * sin(6 * $latitudeOrigin));
670
671 9
        $A = $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue() * cos($latitude);
672 9
        $T = tan($latitude) ** 2;
673 9
        $C = $e2 * cos($latitude) ** 2 / (1 - $e2);
674 9
        $nu = $a / sqrt(1 - $e2 * (sin($latitude) ** 2));
675 9
        $X = $M - $MO + $nu * tan($latitude) * ($A ** 2 / 2 + (5 - $T + 6 * $C) * $A ** 4 / 24);
676
677 9
        $easting = $falseEasting->asMetres()->getValue() + $nu * ($A - $T * $A ** 3 / 6 - (8 - $T + 8 * $C) * $T * $A ** 5 / 120);
678 9
        $northing = $falseNorthing->asMetres()->getValue() + $X;
679
680 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
681
    }
682
683
    /**
684
     * Hyperbolic Cassini-Soldner.
685
     */
686 18
    public function hyperbolicCassiniSoldner(
687
        Projected $to,
688
        Angle $latitudeOfNaturalOrigin,
689
        Angle $longitudeOfNaturalOrigin,
690
        Length $falseEasting,
691
        Length $falseNorthing
692
    ): ProjectedPoint {
693 18
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
694 18
        $latitude = $this->latitude->asRadians()->getValue();
695 18
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
696 18
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
697 18
        $e = $ellipsoid->getEccentricity();
698 18
        $e2 = $ellipsoid->getEccentricitySquared();
699 18
        $e4 = $e ** 4;
700 18
        $e6 = $e ** 6;
701
702 18
        $M = $a * ((1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256) * $latitude - (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024) * sin(2 * $latitude) + (15 * $e4 / 256 + 45 * $e6 / 1024) * sin(4 * $latitude) - (35 * $e6 / 3072) * sin(6 * $latitude));
703 18
        $MO = $a * ((1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256) * $latitudeOrigin - (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024) * sin(2 * $latitudeOrigin) + (15 * $e4 / 256 + 45 * $e6 / 1024) * sin(4 * $latitudeOrigin) - (35 * $e6 / 3072) * sin(6 * $latitudeOrigin));
704
705 18
        $A = $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue() * cos($latitude);
706 18
        $T = tan($latitude) ** 2;
707 18
        $C = $e2 * cos($latitude) ** 2 / (1 - $e2);
708 18
        $nu = $a / sqrt(1 - $e2 * (sin($latitude) ** 2));
709 18
        $rho = $a * (1 - $e2) / (1 - $e2 * sin($latitude) ** 2) ** (3 / 2);
710 18
        $X = $M - $MO + $nu * tan($latitude) * ($A ** 2 / 2 + (5 - $T + 6 * $C) * $A ** 4 / 24);
711
712 18
        $easting = $falseEasting->asMetres()->getValue() + $nu * ($A - $T * $A ** 3 / 6 - (8 - $T + 8 * $C) * $T * $A ** 5 / 120);
713 18
        $northing = $falseNorthing->asMetres()->getValue() + $X - ($X ** 3 / (6 * $rho * $nu));
714
715 18
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
716
    }
717
718
    /**
719
     * Colombia Urban.
720
     */
721 9
    public function columbiaUrban(
722
        Projected $to,
723
        Angle $latitudeOfNaturalOrigin,
724
        Angle $longitudeOfNaturalOrigin,
725
        Length $falseEasting,
726
        Length $falseNorthing,
727
        Length $projectionPlaneOriginHeight
728
    ): ProjectedPoint {
729 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
730 9
        $latitude = $this->latitude->asRadians()->getValue();
731 9
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
732 9
        $heightOrigin = $projectionPlaneOriginHeight->asMetres()->getValue();
733 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
734 9
        $e2 = $ellipsoid->getEccentricitySquared();
735
736 9
        $rhoOrigin = $a * (1 - $e2) / (1 - $e2 * sin($latitudeOrigin) ** 2) ** (3 / 2);
737 9
        $rhoMid = $a * (1 - $e2) / (1 - $e2 * sin(($latitude + $latitudeOrigin) / 2) ** 2) ** (3 / 2);
738
739 9
        $nu = $a / sqrt(1 - $e2 * (sin($latitude) ** 2));
740 9
        $nuOrigin = $a / sqrt(1 - $e2 * (sin($latitudeOrigin) ** 2));
741
742 9
        $A = 1 + $heightOrigin / $nuOrigin;
743 9
        $B = tan($latitudeOrigin) / (2 * $rhoOrigin * $nuOrigin);
744 9
        $G = 1 + $heightOrigin / $rhoMid;
745
746 9
        $easting = $falseEasting->asMetres()->getValue() + $A * $nu * cos($latitude) * $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue();
747 9
        $northing = $falseNorthing->asMetres()->getValue() + $G * $rhoOrigin * (($latitude - $latitudeOrigin) + ($B * $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue() ** 2 * $nu ** 2 * cos($latitude) ** 2));
748
749 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
750
    }
751
752
    /**
753
     * Equal Earth.
754
     */
755 9
    public function equalEarth(
756
        Projected $to,
757
        Angle $longitudeOfNaturalOrigin,
758
        Length $falseEasting,
759
        Length $falseNorthing
760
    ): ProjectedPoint {
761 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
762 9
        $latitude = $this->latitude->asRadians()->getValue();
763 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
764 9
        $e = $ellipsoid->getEccentricity();
765 9
        $e2 = $ellipsoid->getEccentricitySquared();
766
767 9
        $q = (1 - $e2) * ((sin($latitude) / (1 - $e2 * sin($latitude) ** 2)) - (1 / (2 * $e) * log((1 - $e * sin($latitude)) / (1 + $e * sin($latitude)))));
768 9
        $qP = (1 - $e2) * ((1 / (1 - $e2)) - (1 / (2 * $e) * log((1 - $e) / (1 + $e))));
769 9
        $beta = self::asin($q / $qP);
770 9
        $theta = self::asin(sin($beta) * sqrt(3) / 2);
771 9
        $Rq = $a * sqrt($qP / 2);
772
773 9
        $easting = $falseEasting->asMetres()->getValue() + ($Rq * 2 * $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue() * cos($theta)) / (sqrt(3) * (1.340264 - 0.243318 * $theta ** 2 + $theta ** 6 * (0.006251 + 0.034164 * $theta ** 2)));
774 9
        $northing = $falseNorthing->asMetres()->getValue() + $Rq * $theta * (1.340264 - 0.081106 * $theta ** 2 + $theta ** 6 * (0.000893 + 0.003796 * $theta ** 2));
775
776 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
777
    }
778
779
    /**
780
     * Equidistant Cylindrical
781
     * See method code 1029 for spherical development. See also Pseudo Plate Carree, method code 9825.
782
     */
783 9
    public function equidistantCylindrical(
784
        Projected $to,
785
        Angle $latitudeOf1stStandardParallel,
786
        Angle $longitudeOfNaturalOrigin,
787
        Length $falseEasting,
788
        Length $falseNorthing
789
    ): ProjectedPoint {
790 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
791 9
        $latitude = $this->latitude->asRadians()->getValue();
792 9
        $latitudeFirstParallel = $latitudeOf1stStandardParallel->asRadians()->getValue();
793 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
794 9
        $e = $ellipsoid->getEccentricity();
795 9
        $e2 = $ellipsoid->getEccentricitySquared();
796 9
        $e4 = $e ** 4;
797 9
        $e6 = $e ** 6;
798 9
        $e8 = $e ** 8;
799 9
        $e10 = $e ** 10;
800 9
        $e12 = $e ** 12;
801 9
        $e14 = $e ** 14;
802
803 9
        $nu1 = $a / sqrt(1 - $e2 * sin($latitudeFirstParallel) ** 2);
804
805 9
        $M = $a * (
806 9
            (1 - 1 / 4 * $e2 - 3 / 64 * $e4 - 5 / 256 * $e6 - 175 / 16384 * $e8 - 441 / 65536 * $e10 - 4851 / 1048576 * $e12 - 14157 / 4194304 * $e14) * $latitude +
807 9
            (-3 / 8 * $e2 - 3 / 32 * $e4 - 45 / 1024 * $e6 - 105 / 4096 * $e8 - 2205 / 131072 * $e10 - 6237 / 524288 * $e12 - 297297 / 33554432 * $e14) * sin(2 * $latitude) +
808 9
            (15 / 256 * $e4 + 45 / 1024 * $e ** 6 + 525 / 16384 * $e ** 8 + 1575 / 65536 * $e10 + 155925 / 8388608 * $e12 + 495495 / 33554432 * $e14) * sin(4 * $latitude) +
809 9
            (-35 / 3072 * $e6 - 175 / 12288 * $e8 - 3675 / 262144 * $e10 - 13475 / 1048576 * $e12 - 385385 / 33554432 * $e14) * sin(6 * $latitude) +
810 9
            (315 / 131072 * $e8 + 2205 / 524288 * $e10 + 43659 / 8388608 * $e12 + 189189 / 33554432 * $e14) * sin(8 * $latitude) +
811 9
            (-693 / 1310720 * $e10 - 6537 / 5242880 * $e12 - 297297 / 167772160 * $e14) * sin(10 * $latitude) +
812 9
            (1001 / 8388608 * $e12 + 11011 / 33554432 * $e14) * sin(12 * $latitude) +
813 9
            (-6435 / 234881024 * $e ** 14) * sin(14 * $latitude)
814
        );
815
816 9
        $easting = $falseEasting->asMetres()->getValue() + $nu1 * cos($latitudeFirstParallel) * $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue();
817 9
        $northing = $falseNorthing->asMetres()->getValue() + $M;
818
819 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
820
    }
821
822
    /**
823
     * Guam Projection
824
     * Simplified form of Oblique Azimuthal Equidistant projection method.
825
     */
826 9
    public function guamProjection(
827
        Projected $to,
828
        Angle $latitudeOfNaturalOrigin,
829
        Angle $longitudeOfNaturalOrigin,
830
        Length $falseEasting,
831
        Length $falseNorthing
832
    ): ProjectedPoint {
833 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
834 9
        $latitude = $this->latitude->asRadians()->getValue();
835 9
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
836 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
837 9
        $e = $ellipsoid->getEccentricity();
838 9
        $e2 = $ellipsoid->getEccentricitySquared();
839 9
        $e4 = $e ** 4;
840 9
        $e6 = $e ** 6;
841
842 9
        $M = $a * ((1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256) * $latitude - (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024) * sin(2 * $latitude) + (15 * $e4 / 256 + 45 * $e6 / 1024) * sin(4 * $latitude) - (35 * $e6 / 3072) * sin(6 * $latitude));
843 9
        $MO = $a * ((1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256) * $latitudeOrigin - (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024) * sin(2 * $latitudeOrigin) + (15 * $e4 / 256 + 45 * $e6 / 1024) * sin(4 * $latitudeOrigin) - (35 * $e6 / 3072) * sin(6 * $latitudeOrigin));
844 9
        $x = ($a * $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue() * cos($latitude)) / sqrt(1 - $e2 * sin($latitude) ** 2);
845
846 9
        $easting = $falseEasting->asMetres()->getValue() + $x;
847 9
        $northing = $falseNorthing->asMetres()->getValue() + $M - $MO + ($x ** 2 * tan($latitude) * sqrt(1 - $e2 * sin($latitude) ** 2) / (2 * $a));
848
849 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
850
    }
851
852
    /**
853
     * Krovak.
854
     */
855 36
    public function krovak(
856
        Projected $to,
857
        Angle $latitudeOfProjectionCentre,
858
        Angle $longitudeOfOrigin,
859
        Angle $coLatitudeOfConeAxis,
860
        Angle $latitudeOfPseudoStandardParallel,
861
        Scale $scaleFactorOnPseudoStandardParallel,
862
        Length $falseEasting,
863
        Length $falseNorthing
864
    ): ProjectedPoint {
865 36
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
866 36
        $longitudeOffset = $to->getDatum()->getPrimeMeridian()->getGreenwichLongitude()->asRadians()->getValue() - $this->getCRS()->getDatum()->getPrimeMeridian()->getGreenwichLongitude()->asRadians()->getValue();
867 36
        $latitude = $this->latitude->asRadians()->getValue();
868 36
        $longitude = $this->longitude->asRadians()->getValue() - $longitudeOffset;
869 36
        $latitudeC = $latitudeOfProjectionCentre->asRadians()->getValue();
870 36
        $longitudeO = $longitudeOfOrigin->asRadians()->getValue();
871 36
        $alphaC = $coLatitudeOfConeAxis->asRadians()->getValue();
872 36
        $latitudeP = $latitudeOfPseudoStandardParallel->asRadians()->getValue();
873 36
        $kP = $scaleFactorOnPseudoStandardParallel->asUnity()->getValue();
874 36
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
875 36
        $e = $ellipsoid->getEccentricity();
876 36
        $e2 = $ellipsoid->getEccentricitySquared();
877
878 36
        $A = $a * sqrt(1 - $e2) / (1 - $e2 * sin($latitudeC) ** 2);
879 36
        $B = sqrt(1 + $e2 * cos($latitudeC) ** 4 / (1 - $e2));
880 36
        $upsilonO = self::asin(sin($latitudeC) / $B);
881 36
        $tO = tan(M_PI / 4 + $upsilonO / 2) * ((1 + $e * sin($latitudeC)) / (1 - $e * sin($latitudeC))) ** ($e * $B / 2) / (tan(M_PI / 4 + $latitudeC / 2) ** $B);
882 36
        $n = sin($latitudeP);
883 36
        $rO = $kP * $A / tan($latitudeP);
884
885 36
        $U = 2 * (atan($tO * tan($latitude / 2 + M_PI / 4) ** $B / ((1 + $e * sin($latitude)) / (1 - $e * sin($latitude))) ** ($e * $B / 2)) - M_PI / 4);
886 36
        $V = $B * ($longitudeO - $longitude);
887 36
        $T = self::asin(cos($alphaC) * sin($U) + sin($alphaC) * cos($U) * cos($V));
888 36
        $D = atan2(cos($U) * sin($V) / cos($T), ((cos($alphaC) * sin($T) - sin($U)) / (sin($alphaC) * cos($T))));
889 36
        $theta = $n * $D;
890 36
        $r = $rO * tan(M_PI / 4 + $latitudeP / 2) ** $n / tan($T / 2 + M_PI / 4) ** $n;
891 36
        $X = $r * cos($theta);
892 36
        $Y = $r * sin($theta);
893
894 36
        $westing = $Y + $falseEasting->asMetres()->getValue();
895 36
        $southing = $X + $falseNorthing->asMetres()->getValue();
896
897 36
        return ProjectedPoint::create(new Metre(-$westing), new Metre(-$southing), new Metre($westing), new Metre($southing), $to, $this->epoch);
898
    }
899
900
    /**
901
     * Krovak Modified
902
     * Incorporates a polynomial transformation which is defined to be exact and for practical purposes is considered
903
     * to be a map projection.
904
     */
905 18
    public function krovakModified(
906
        Projected $to,
907
        Angle $latitudeOfProjectionCentre,
908
        Angle $longitudeOfOrigin,
909
        Angle $coLatitudeOfConeAxis,
910
        Angle $latitudeOfPseudoStandardParallel,
911
        Scale $scaleFactorOnPseudoStandardParallel,
912
        Length $falseEasting,
913
        Length $falseNorthing,
914
        Length $ordinate1OfEvaluationPoint,
915
        Length $ordinate2OfEvaluationPoint,
916
        Coefficient $C1,
917
        Coefficient $C2,
918
        Coefficient $C3,
919
        Coefficient $C4,
920
        Coefficient $C5,
921
        Coefficient $C6,
922
        Coefficient $C7,
923
        Coefficient $C8,
924
        Coefficient $C9,
925
        Coefficient $C10
926
    ): ProjectedPoint {
927 18
        $asKrovak = $this->krovak($to, $latitudeOfProjectionCentre, $longitudeOfOrigin, $coLatitudeOfConeAxis, $latitudeOfPseudoStandardParallel, $scaleFactorOnPseudoStandardParallel, new Metre(0), new Metre(0));
928
929 18
        $westing = $asKrovak->getWesting()->asMetres()->getValue();
930 18
        $southing = $asKrovak->getSouthing()->asMetres()->getValue();
931 18
        $c1 = $C1->asUnity()->getValue();
932 18
        $c2 = $C2->asUnity()->getValue();
933 18
        $c3 = $C3->asUnity()->getValue();
934 18
        $c4 = $C4->asUnity()->getValue();
935 18
        $c5 = $C5->asUnity()->getValue();
936 18
        $c6 = $C6->asUnity()->getValue();
937 18
        $c7 = $C7->asUnity()->getValue();
938 18
        $c8 = $C8->asUnity()->getValue();
939 18
        $c9 = $C9->asUnity()->getValue();
940 18
        $c10 = $C10->asUnity()->getValue();
941
942 18
        $Xr = $southing - $ordinate1OfEvaluationPoint->asMetres()->getValue();
943 18
        $Yr = $westing - $ordinate2OfEvaluationPoint->asMetres()->getValue();
944
945 18
        $dX = $c1 + $c3 * $Xr - $c4 * $Yr - 2 * $c6 * $Xr * $Yr + $c5 * ($Xr ** 2 - $Yr ** 2) + $c7 * $Xr * ($Xr ** 2 - 3 * $Yr ** 2) - $c8 * $Yr * (3 * $Xr ** 2 - $Yr ** 2) + 4 * $c9 * $Xr * $Yr * ($Xr ** 2 - $Yr ** 2) + $c10 * ($Xr ** 4 + $Yr ** 4 - 6 * $Xr ** 2 * $Yr ** 2);
946 18
        $dY = $c2 + $c3 * $Yr + $c4 * $Xr + 2 * $c5 * $Xr * $Yr + $c6 * ($Xr ** 2 - $Yr ** 2) + $c8 * $Xr * ($Xr ** 2 - 3 * $Yr ** 2) + $c7 * $Yr * (3 * $Xr ** 2 - $Yr ** 2) - 4 * $c10 * $Xr * $Yr * ($Xr ** 2 - $Yr ** 2) + $c9 * ($Xr ** 4 + $Yr ** 4 - 6 * $Xr ** 2 * $Yr ** 2);
947
948 18
        $westing += $falseEasting->asMetres()->getValue() - $dY;
949 18
        $southing += $falseNorthing->asMetres()->getValue() - $dX;
950
951 18
        return ProjectedPoint::create(new Metre(-$westing), new Metre(-$southing), new Metre($westing), new Metre($southing), $to, $this->epoch);
952
    }
953
954
    /**
955
     * Lambert Azimuthal Equal Area
956
     * This is the ellipsoidal form of the projection.
957
     */
958 9
    public function lambertAzimuthalEqualArea(
959
        Projected $to,
960
        Angle $latitudeOfNaturalOrigin,
961
        Angle $longitudeOfNaturalOrigin,
962
        Length $falseEasting,
963
        Length $falseNorthing
964
    ): ProjectedPoint {
965 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
966 9
        $latitude = $this->latitude->asRadians()->getValue();
967 9
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
968 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
969 9
        $e = $ellipsoid->getEccentricity();
970 9
        $e2 = $ellipsoid->getEccentricitySquared();
971
972 9
        $q = (1 - $e2) * ((sin($latitude) / (1 - $e2 * sin($latitude) ** 2)) - ((1 / (2 * $e)) * log((1 - $e * sin($latitude)) / (1 + $e * sin($latitude)))));
973 9
        $qO = (1 - $e2) * ((sin($latitudeOrigin) / (1 - $e2 * sin($latitudeOrigin) ** 2)) - ((1 / (2 * $e)) * log((1 - $e * sin($latitudeOrigin)) / (1 + $e * sin($latitudeOrigin)))));
974 9
        $qP = (1 - $e2) * ((1 / (1 - $e2)) - ((1 / (2 * $e)) * log((1 - $e) / (1 + $e))));
975 9
        $beta = self::asin($q / $qP);
976 9
        $betaO = self::asin($qO / $qP);
977 9
        $Rq = $a * sqrt($qP / 2);
978 9
        $B = $Rq * sqrt(2 / (1 + sin($betaO) * sin($beta) + (cos($betaO) * cos($beta) * cos($this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue()))));
979 9
        $D = $a * (cos($latitudeOrigin) / sqrt(1 - $e2 * sin($latitudeOrigin) ** 2)) / ($Rq * cos($betaO));
980
981 9
        $easting = $falseEasting->asMetres()->getValue() + (($B * $D) * (cos($beta) * sin($this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue())));
982 9
        $northing = $falseNorthing->asMetres()->getValue() + ($B / $D) * ((cos($betaO) * sin($beta)) - (sin($betaO) * cos($beta) * cos($this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue())));
983
984 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
985
    }
986
987
    /**
988
     * Lambert Azimuthal Equal Area (Spherical)
989
     * This is the spherical form of the projection.  See coordinate operation method Lambert Azimuthal Equal Area
990
     * (code 9820) for ellipsoidal form.  Differences of several tens of metres result from comparison of the two
991
     * methods.
992
     */
993 9
    public function lambertAzimuthalEqualAreaSpherical(
994
        Projected $to,
995
        Angle $latitudeOfNaturalOrigin,
996
        Angle $longitudeOfNaturalOrigin,
997
        Length $falseEasting,
998
        Length $falseNorthing
999
    ): ProjectedPoint {
1000 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1001 9
        $latitude = $this->latitude->asRadians()->getValue();
1002 9
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
1003 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1004
1005 9
        $k = sqrt(2 / (1 + sin($latitudeOrigin) * sin($latitude) + cos($latitudeOrigin) * cos($latitude) * cos($this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue())));
1006
1007 9
        $easting = $falseEasting->asMetres()->getValue() + ($a * $k * cos($latitude) * sin($this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue()));
1008 9
        $northing = $falseNorthing->asMetres()->getValue() + ($a * $k * (cos($latitudeOrigin) * sin($latitude) - sin($latitudeOrigin) * cos($latitude) * cos($this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue())));
1009
1010 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1011
    }
1012
1013
    /**
1014
     * Lambert Conic Conformal (1SP).
1015
     */
1016 9
    public function lambertConicConformal1SP(
1017
        Projected $to,
1018
        Angle $latitudeOfNaturalOrigin,
1019
        Angle $longitudeOfNaturalOrigin,
1020
        Scale $scaleFactorAtNaturalOrigin,
1021
        Length $falseEasting,
1022
        Length $falseNorthing
1023
    ): ProjectedPoint {
1024 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1025 9
        $latitude = $this->latitude->asRadians()->getValue();
1026 9
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
1027 9
        $kO = $scaleFactorAtNaturalOrigin->asUnity()->getValue();
1028 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1029 9
        $e = $ellipsoid->getEccentricity();
1030 9
        $e2 = $ellipsoid->getEccentricitySquared();
1031
1032 9
        $mO = cos($latitudeOrigin) / sqrt(1 - $e2 * sin($latitudeOrigin) ** 2);
1033 9
        $tO = tan(M_PI / 4 - $latitudeOrigin / 2) / ((1 - $e * sin($latitudeOrigin)) / (1 + $e * sin($latitudeOrigin))) ** ($e / 2);
1034 9
        $t = tan(M_PI / 4 - $latitude / 2) / ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2);
1035 9
        $n = sin($latitudeOrigin);
1036 9
        $F = $mO / ($n * $tO ** $n);
1037 9
        $rO = $a * $F * $tO ** $n * $kO;
1038 9
        $r = $a * $F * $t ** $n * $kO;
1039 9
        $theta = $n * $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue();
1040
1041 9
        $easting = $falseEasting->asMetres()->getValue() + $r * sin($theta);
1042 9
        $northing = $falseNorthing->asMetres()->getValue() + $rO - $r * cos($theta);
1043
1044 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1045
    }
1046
1047
    /**
1048
     * Lambert Conic Conformal (1SP) Variant B.
1049
     */
1050
    public function lambertConicConformal1SPVariantB(
1051
        Projected $to,
1052
        Angle $latitudeOfNaturalOrigin,
1053
        Scale $scaleFactorAtNaturalOrigin,
1054
        Angle $latitudeOfFalseOrigin,
1055
        Angle $longitudeOfFalseOrigin,
1056
        Length $eastingAtFalseOrigin,
1057
        Length $northingAtFalseOrigin
1058
    ): ProjectedPoint {
1059
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1060
        $latitude = $this->latitude->asRadians()->getValue();
1061
        $latitudeNaturalOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
1062
        $latitudeFalseOrigin = $latitudeOfFalseOrigin->asRadians()->getValue();
1063
        $kO = $scaleFactorAtNaturalOrigin->asUnity()->getValue();
1064
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1065
        $e = $ellipsoid->getEccentricity();
1066
        $e2 = $ellipsoid->getEccentricitySquared();
1067
1068
        $mO = cos($latitudeNaturalOrigin) / sqrt(1 - $e2 * sin($latitudeNaturalOrigin) ** 2);
1069
        $tO = tan(M_PI / 4 - $latitudeNaturalOrigin / 2) / ((1 - $e * sin($latitudeNaturalOrigin)) / (1 + $e * sin($latitudeNaturalOrigin))) ** ($e / 2);
1070
        $tF = tan(M_PI / 4 - $latitudeFalseOrigin / 2) / ((1 - $e * sin($latitudeFalseOrigin)) / (1 + $e * sin($latitudeFalseOrigin))) ** ($e / 2);
1071
        $t = tan(M_PI / 4 - $latitude / 2) / ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2);
1072
        $n = sin($latitudeNaturalOrigin);
1073
        $F = $mO / ($n * $tO ** $n);
1074
        $rF = $a * $F * $tF ** $n * $kO;
1075
        $r = $a * $F * $t ** $n * $kO;
1076
        $theta = $n * $this->normaliseLongitude($this->longitude->subtract($longitudeOfFalseOrigin))->asRadians()->getValue();
1077
1078
        $easting = $eastingAtFalseOrigin->asMetres()->getValue() + $r * sin($theta);
1079
        $northing = $northingAtFalseOrigin->asMetres()->getValue() + $rF - $r * cos($theta);
1080
1081
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1082
    }
1083
1084
    /**
1085
     * Lambert Conic Conformal (2SP Belgium)
1086
     * In 2000 this modification was replaced through use of the regular Lambert Conic Conformal (2SP) method [9802]
1087
     * with appropriately modified parameter values.
1088
     */
1089 9
    public function lambertConicConformal2SPBelgium(
1090
        Projected $to,
1091
        Angle $latitudeOfFalseOrigin,
1092
        Angle $longitudeOfFalseOrigin,
1093
        Angle $latitudeOf1stStandardParallel,
1094
        Angle $latitudeOf2ndStandardParallel,
1095
        Length $eastingAtFalseOrigin,
1096
        Length $northingAtFalseOrigin
1097
    ): ProjectedPoint {
1098 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1099 9
        $latitude = $this->latitude->asRadians()->getValue();
1100 9
        $phiF = $latitudeOfFalseOrigin->asRadians()->getValue();
1101 9
        $phi1 = $latitudeOf1stStandardParallel->asRadians()->getValue();
1102 9
        $phi2 = $latitudeOf2ndStandardParallel->asRadians()->getValue();
1103 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1104 9
        $e = $ellipsoid->getEccentricity();
1105 9
        $e2 = $ellipsoid->getEccentricitySquared();
1106
1107 9
        $m1 = cos($phi1) / sqrt(1 - $e2 * sin($phi1) ** 2);
1108 9
        $m2 = cos($phi2) / sqrt(1 - $e2 * sin($phi2) ** 2);
1109 9
        $t = tan(M_PI / 4 - $latitude / 2) / ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2);
1110 9
        $t1 = tan(M_PI / 4 - $phi1 / 2) / ((1 - $e * sin($phi1)) / (1 + $e * sin($phi1))) ** ($e / 2);
1111 9
        $t2 = tan(M_PI / 4 - $phi2 / 2) / ((1 - $e * sin($phi2)) / (1 + $e * sin($phi2))) ** ($e / 2);
1112 9
        $tF = tan(M_PI / 4 - $phiF / 2) / ((1 - $e * sin($phiF)) / (1 + $e * sin($phiF))) ** ($e / 2);
1113 9
        $n = (log($m1) - log($m2)) / (log($t1) - log($t2));
1114 9
        $F = $m1 / ($n * $t1 ** $n);
1115 9
        $r = $a * $F * $t ** $n;
1116 9
        $rF = $a * $F * $tF ** $n;
1117 9
        if (is_nan($rF)) {
1118 9
            $rF = 0;
1119
        }
1120 9
        $theta = ($n * $this->normaliseLongitude($this->longitude->subtract($longitudeOfFalseOrigin))->asRadians()->getValue()) - (new ArcSecond(29.2985))->asRadians()->getValue();
1121
1122 9
        $easting = $eastingAtFalseOrigin->asMetres()->getValue() + $r * sin($theta);
1123 9
        $northing = $northingAtFalseOrigin->asMetres()->getValue() + $rF - $r * cos($theta);
1124
1125 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1126
    }
1127
1128
    /**
1129
     * Lambert Conic Conformal (2SP Michigan).
1130
     */
1131 9
    public function lambertConicConformal2SPMichigan(
1132
        Projected $to,
1133
        Angle $latitudeOfFalseOrigin,
1134
        Angle $longitudeOfFalseOrigin,
1135
        Angle $latitudeOf1stStandardParallel,
1136
        Angle $latitudeOf2ndStandardParallel,
1137
        Length $eastingAtFalseOrigin,
1138
        Length $northingAtFalseOrigin,
1139
        Scale $ellipsoidScalingFactor
1140
    ): ProjectedPoint {
1141 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1142 9
        $latitude = $this->latitude->asRadians()->getValue();
1143 9
        $phiF = $latitudeOfFalseOrigin->asRadians()->getValue();
1144 9
        $phi1 = $latitudeOf1stStandardParallel->asRadians()->getValue();
1145 9
        $phi2 = $latitudeOf2ndStandardParallel->asRadians()->getValue();
1146 9
        $K = $ellipsoidScalingFactor->asUnity()->getValue();
1147 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1148 9
        $e = $ellipsoid->getEccentricity();
1149 9
        $e2 = $ellipsoid->getEccentricitySquared();
1150
1151 9
        $m1 = cos($phi1) / sqrt(1 - $e2 * sin($phi1) ** 2);
1152 9
        $m2 = cos($phi2) / sqrt(1 - $e2 * sin($phi2) ** 2);
1153 9
        $t = tan(M_PI / 4 - $latitude / 2) / ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2);
1154 9
        $t1 = tan(M_PI / 4 - $phi1 / 2) / ((1 - $e * sin($phi1)) / (1 + $e * sin($phi1))) ** ($e / 2);
1155 9
        $t2 = tan(M_PI / 4 - $phi2 / 2) / ((1 - $e * sin($phi2)) / (1 + $e * sin($phi2))) ** ($e / 2);
1156 9
        $tF = tan(M_PI / 4 - $phiF / 2) / ((1 - $e * sin($phiF)) / (1 + $e * sin($phiF))) ** ($e / 2);
1157 9
        $n = (log($m1) - log($m2)) / (log($t1) - log($t2));
1158 9
        $F = $m1 / ($n * $t1 ** $n);
1159 9
        $r = $a * $K * $F * $t ** $n;
1160 9
        $rF = $a * $K * $F * $tF ** $n;
1161 9
        $theta = $n * $this->normaliseLongitude($this->longitude->subtract($longitudeOfFalseOrigin))->asRadians()->getValue();
1162
1163 9
        $easting = $eastingAtFalseOrigin->asMetres()->getValue() + $r * sin($theta);
1164 9
        $northing = $northingAtFalseOrigin->asMetres()->getValue() + $rF - $r * cos($theta);
1165
1166 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1167
    }
1168
1169
    /**
1170
     * Lambert Conic Conformal (2SP).
1171
     */
1172 9
    public function lambertConicConformal2SP(
1173
        Projected $to,
1174
        Angle $latitudeOfFalseOrigin,
1175
        Angle $longitudeOfFalseOrigin,
1176
        Angle $latitudeOf1stStandardParallel,
1177
        Angle $latitudeOf2ndStandardParallel,
1178
        Length $eastingAtFalseOrigin,
1179
        Length $northingAtFalseOrigin
1180
    ): ProjectedPoint {
1181 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1182 9
        $latitude = $this->latitude->asRadians()->getValue();
1183 9
        $phiF = $latitudeOfFalseOrigin->asRadians()->getValue();
1184 9
        $phi1 = $latitudeOf1stStandardParallel->asRadians()->getValue();
1185 9
        $phi2 = $latitudeOf2ndStandardParallel->asRadians()->getValue();
1186 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1187 9
        $e = $ellipsoid->getEccentricity();
1188 9
        $e2 = $ellipsoid->getEccentricitySquared();
1189
1190 9
        $m1 = cos($phi1) / sqrt(1 - $e2 * sin($phi1) ** 2);
1191 9
        $m2 = cos($phi2) / sqrt(1 - $e2 * sin($phi2) ** 2);
1192 9
        $t = tan(M_PI / 4 - $latitude / 2) / ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2);
1193 9
        $t1 = tan(M_PI / 4 - $phi1 / 2) / ((1 - $e * sin($phi1)) / (1 + $e * sin($phi1))) ** ($e / 2);
1194 9
        $t2 = tan(M_PI / 4 - $phi2 / 2) / ((1 - $e * sin($phi2)) / (1 + $e * sin($phi2))) ** ($e / 2);
1195 9
        $tF = tan(M_PI / 4 - $phiF / 2) / ((1 - $e * sin($phiF)) / (1 + $e * sin($phiF))) ** ($e / 2);
1196 9
        $n = (log($m1) - log($m2)) / (log($t1) - log($t2));
1197 9
        $F = $m1 / ($n * $t1 ** $n);
1198 9
        $r = $a * $F * $t ** $n;
1199 9
        $rF = $a * $F * $tF ** $n;
1200 9
        $theta = $n * $this->normaliseLongitude($this->longitude->subtract($longitudeOfFalseOrigin))->asRadians()->getValue();
1201
1202 9
        $easting = $eastingAtFalseOrigin->asMetres()->getValue() + $r * sin($theta);
1203 9
        $northing = $northingAtFalseOrigin->asMetres()->getValue() + $rF - $r * cos($theta);
1204
1205 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1206
    }
1207
1208
    /**
1209
     * Lambert Conic Conformal (West Orientated).
1210
     */
1211
    public function lambertConicConformalWestOrientated(
1212
        Projected $to,
1213
        Angle $latitudeOfNaturalOrigin,
1214
        Angle $longitudeOfNaturalOrigin,
1215
        Scale $scaleFactorAtNaturalOrigin,
1216
        Length $falseEasting,
1217
        Length $falseNorthing
1218
    ): ProjectedPoint {
1219
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1220
        $latitude = $this->latitude->asRadians()->getValue();
1221
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
1222
        $kO = $scaleFactorAtNaturalOrigin->asUnity()->getValue();
1223
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1224
        $e = $ellipsoid->getEccentricity();
1225
        $e2 = $ellipsoid->getEccentricitySquared();
1226
1227
        $mO = cos($latitudeOrigin) / sqrt(1 - $e2 * sin($latitudeOrigin) ** 2);
1228
        $tO = tan(M_PI / 4 - $latitudeOrigin / 2) / ((1 - $e * sin($latitudeOrigin)) / (1 + $e * sin($latitudeOrigin))) ** ($e / 2);
1229
        $t = tan(M_PI / 4 - $latitude / 2) / ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2);
1230
        $n = sin($latitudeOrigin);
1231
        $F = $mO / ($n * $tO ** $n);
1232
        $rO = $a * $F * $tO ** $n ** $kO;
1233
        $r = $a * $F * $t ** $n ** $kO;
1234
        $theta = $n * $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue();
1235
1236
        $westing = $falseEasting->asMetres()->getValue() - $r * sin($theta);
1237
        $northing = $falseNorthing->asMetres()->getValue() + $rO - $r * cos($theta);
1238
1239
        return ProjectedPoint::create(new Metre(-$westing), new Metre($northing), new Metre($westing), new Metre(-$northing), $to, $this->epoch);
1240
    }
1241
1242
    /**
1243
     * Lambert Conic Near-Conformal
1244
     * The Lambert Near-Conformal projection is derived from the Lambert Conformal Conic projection by truncating the
1245
     * series expansion of the projection formulae.
1246
     */
1247 9
    public function lambertConicNearConformal(
1248
        Projected $to,
1249
        Angle $latitudeOfNaturalOrigin,
1250
        Angle $longitudeOfNaturalOrigin,
1251
        Scale $scaleFactorAtNaturalOrigin,
1252
        Length $falseEasting,
1253
        Length $falseNorthing
1254
    ): ProjectedPoint {
1255 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1256 9
        $latitude = $this->latitude->asRadians()->getValue();
1257 9
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
1258 9
        $kO = $scaleFactorAtNaturalOrigin->asUnity()->getValue();
1259 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1260 9
        $e2 = $ellipsoid->getEccentricitySquared();
1261 9
        $f = $ellipsoid->getInverseFlattening();
1262
1263 9
        $n = $f / (2 - $f);
1264 9
        $rhoO = $a * (1 - $e2) / (1 - $e2 * sin($latitudeOrigin) ** 2) ** (3 / 2);
1265 9
        $nuO = $a / sqrt(1 - $e2 * (sin($latitudeOrigin) ** 2));
1266 9
        $A = 1 / (6 * $rhoO * $nuO);
1267 9
        $APrime = $a * (1 - $n + 5 * ($n ** 2 - $n ** 3) / 4 + 81 * ($n ** 4 - $n ** 5) / 64);
1268 9
        $BPrime = 3 * $a * ($n - $n ** 2 + 7 * ($n ** 3 - $n ** 4) / 8 + 55 * $n ** 5 / 64) / 2;
1269 9
        $CPrime = 15 * $a * ($n ** 2 - $n ** 3 + 3 * ($n ** 4 - $n ** 5) / 4) / 16;
1270 9
        $DPrime = 35 * $a * ($n ** 3 - $n ** 4 + 11 * $n ** 5 / 16) / 48;
1271 9
        $EPrime = 315 * $a * ($n ** 4 - $n ** 5) / 512;
1272 9
        $rO = $kO * $nuO / tan($latitudeOrigin);
1273 9
        $sO = $APrime * $latitudeOrigin - $BPrime * sin(2 * $latitudeOrigin) + $CPrime * sin(4 * $latitudeOrigin) - $DPrime * sin(6 * $latitudeOrigin) + $EPrime * sin(8 * $latitudeOrigin);
1274 9
        $s = $APrime * $latitude - $BPrime * sin(2 * $latitude) + $CPrime * sin(4 * $latitude) - $DPrime * sin(6 * $latitude) + $EPrime * sin(8 * $latitude);
1275 9
        $m = $s - $sO;
1276 9
        $M = $kO * ($m + $A * $m ** 3);
1277 9
        $r = $rO - $M;
1278 9
        $theta = $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue() * sin($latitudeOrigin);
1279
1280 9
        $easting = $falseEasting->asMetres()->getValue() + $r * sin($theta);
1281 9
        $northing = $falseNorthing->asMetres()->getValue() + $M + $r * sin($theta) * tan($theta / 2);
1282
1283 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1284
    }
1285
1286
    /**
1287
     * Lambert Cylindrical Equal Area
1288
     * This is the ellipsoidal form of the projection.
1289
     */
1290 9
    public function lambertCylindricalEqualArea(
1291
        Projected $to,
1292
        Angle $latitudeOf1stStandardParallel,
1293
        Angle $longitudeOfNaturalOrigin,
1294
        Length $falseEasting,
1295
        Length $falseNorthing
1296
    ): ProjectedPoint {
1297 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1298 9
        $latitude = $this->latitude->asRadians()->getValue();
1299 9
        $latitudeFirstParallel = $latitudeOf1stStandardParallel->asRadians()->getValue();
1300 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1301 9
        $e = $ellipsoid->getEccentricity();
1302 9
        $e2 = $ellipsoid->getEccentricitySquared();
1303
1304 9
        $k = cos($latitudeFirstParallel) / sqrt(1 - $e2 * sin($latitudeFirstParallel) ** 2);
1305 9
        $q = (1 - $e2) * ((sin($latitude) / (1 - $e2 * sin($latitude) ** 2)) - (1 / (2 * $e)) * log((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))));
1306
1307 9
        $x = $a * $k * $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue();
1308 9
        $y = $a * $q / (2 * $k);
1309
1310 9
        $easting = $falseEasting->asMetres()->getValue() + $x;
1311 9
        $northing = $falseNorthing->asMetres()->getValue() + $y;
1312
1313 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1314
    }
1315
1316
    /**
1317
     * Modified Azimuthal Equidistant
1318
     * Modified form of Oblique Azimuthal Equidistant projection method developed for Polynesian islands. For the
1319
     * distances over which these projections are used (under 800km) this modification introduces no significant error.
1320
     */
1321 9
    public function modifiedAzimuthalEquidistant(
1322
        Projected $to,
1323
        Angle $latitudeOfNaturalOrigin,
1324
        Angle $longitudeOfNaturalOrigin,
1325
        Length $falseEasting,
1326
        Length $falseNorthing
1327
    ): ProjectedPoint {
1328 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1329 9
        $latitude = $this->latitude->asRadians()->getValue();
1330 9
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
1331 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1332 9
        $e = $ellipsoid->getEccentricity();
1333 9
        $e2 = $ellipsoid->getEccentricitySquared();
1334
1335 9
        $nuO = $a / sqrt(1 - $e2 * sin($latitudeOrigin) ** 2);
1336 9
        $nu = $a / sqrt(1 - $e2 * sin($latitude) ** 2);
1337 9
        $psi = atan((1 - $e2) * tan($latitude) + ($e2 * $nuO * sin($latitudeOrigin)) / ($nu * cos($latitude)));
1338 9
        $alpha = atan2(sin($this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue()), (cos($latitudeOrigin) * tan($psi) - sin($latitudeOrigin) * cos($this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue())));
1339 9
        $G = $e * sin($latitudeOrigin) / sqrt(1 - $e2);
1340 9
        $H = $e * cos($latitudeOrigin) * cos($alpha) / sqrt(1 - $e2);
1341
1342 9
        if (sin($alpha) === 0.0) {
1343
            $s = self::asin(cos($latitudeOrigin) * sin($psi) - sin($latitudeOrigin) * cos($alpha)) * cos($alpha) / abs(cos($alpha));
1344
        } else {
1345 9
            $s = self::asin(sin($this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue()) * cos($psi) / sin($alpha));
1346
        }
1347
1348 9
        $c = $nuO * $s * ((1 - $s ** 2 * $H ** 2 * (1 - $H ** 2) / 6) + (($s ** 3 / 8) * $G * $H * (1 - 2 * $H ** 2)) + ($s ** 4 / 120) * ($H ** 2 * (4 - 7 * $H ** 2) - 3 * $G ** 2 * (1 - 7 * $H ** 2)) - (($s ** 5 / 48) * $G * $H));
1349
1350 9
        $easting = $falseEasting->asMetres()->getValue() + $c * sin($alpha);
1351 9
        $northing = $falseNorthing->asMetres()->getValue() + $c * cos($alpha);
1352
1353 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1354
    }
1355
1356
    /**
1357
     * Oblique Stereographic
1358
     * This is not the same as the projection method of the same name in USGS Professional Paper no. 1395, "Map
1359
     * Projections - A Working Manual" by John P. Snyder.
1360
     */
1361 9
    public function obliqueStereographic(
1362
        Projected $to,
1363
        Angle $latitudeOfNaturalOrigin,
1364
        Angle $longitudeOfNaturalOrigin,
1365
        Scale $scaleFactorAtNaturalOrigin,
1366
        Length $falseEasting,
1367
        Length $falseNorthing
1368
    ): ProjectedPoint {
1369 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1370 9
        $latitude = $this->latitude->asRadians()->getValue();
1371 9
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
1372 9
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
1373 9
        $kO = $scaleFactorAtNaturalOrigin->asUnity()->getValue();
1374 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1375 9
        $e = $ellipsoid->getEccentricity();
1376 9
        $e2 = $ellipsoid->getEccentricitySquared();
1377
1378 9
        $rhoOrigin = $a * (1 - $e2) / (1 - $e2 * sin($latitudeOrigin) ** 2) ** (3 / 2);
1379 9
        $nuOrigin = $a / sqrt(1 - $e2 * (sin($latitudeOrigin) ** 2));
1380 9
        $R = sqrt($rhoOrigin * $nuOrigin);
1381
1382 9
        $n = sqrt(1 + ($e2 * cos($latitudeOrigin) ** 4 / (1 - $e2)));
1383 9
        $S1 = (1 + sin($latitudeOrigin)) / (1 - sin($latitudeOrigin));
1384 9
        $S2 = (1 - $e * sin($latitudeOrigin)) / (1 + $e * sin($latitudeOrigin));
1385 9
        $w1 = ($S1 * ($S2 ** $e)) ** $n;
1386 9
        $c = (($n + sin($latitudeOrigin)) * (1 - ($w1 - 1) / ($w1 + 1))) / (($n - sin($latitudeOrigin)) * (1 + ($w1 - 1) / ($w1 + 1)));
1387 9
        $w2 = $c * $w1;
1388 9
        $chiOrigin = self::asin(($w2 - 1) / ($w2 + 1));
1389
1390 9
        $lambda = $n * $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue() + $longitudeOrigin;
1391
1392 9
        $Sa = (1 + sin($latitude)) / (1 - sin($latitude));
1393 9
        $Sb = (1 - $e * sin($latitude)) / (1 + $e * sin($latitude));
1394 9
        $w = $c * ($Sa * ($Sb ** $e)) ** $n;
1395 9
        $chi = self::asin(($w - 1) / ($w + 1));
1396
1397 9
        $B = (1 + sin($chi) * sin($chiOrigin) + cos($chi) * cos($chiOrigin) * cos($lambda - $longitudeOrigin));
1398
1399 9
        $easting = $falseEasting->asMetres()->getValue() + 2 * $R * $kO * cos($chi) * sin($lambda - $longitudeOrigin) / $B;
1400 9
        $northing = $falseNorthing->asMetres()->getValue() + 2 * $R * $kO * (sin($chi) * cos($chiOrigin) - cos($chi) * sin($chiOrigin) * cos($lambda - $longitudeOrigin)) / $B;
1401
1402 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1403
    }
1404
1405
    /**
1406
     * Polar Stereographic (variant A)
1407
     * Latitude of natural origin must be either 90 degrees or -90 degrees (or equivalent in alternative angle unit).
1408
     */
1409 9
    public function polarStereographicVariantA(
1410
        Projected $to,
1411
        Angle $latitudeOfNaturalOrigin,
1412
        Angle $longitudeOfNaturalOrigin,
1413
        Scale $scaleFactorAtNaturalOrigin,
1414
        Length $falseEasting,
1415
        Length $falseNorthing
1416
    ): ProjectedPoint {
1417 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1418 9
        $latitude = $this->latitude->asRadians()->getValue();
1419 9
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
1420 9
        $kO = $scaleFactorAtNaturalOrigin->asUnity()->getValue();
1421 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1422 9
        $e = $ellipsoid->getEccentricity();
1423
1424 9
        if ($latitudeOrigin < 0) {
1425
            $t = tan(M_PI / 4 + $latitude / 2) / (((1 + $e * sin($latitude)) / (1 - $e * sin($latitude))) ** ($e / 2));
1426
        } else {
1427 9
            $t = tan(M_PI / 4 - $latitude / 2) * (((1 + $e * sin($latitude)) / (1 - $e * sin($latitude))) ** ($e / 2));
1428
        }
1429 9
        $rho = 2 * $a * $kO * $t / sqrt((1 + $e) ** (1 + $e) * (1 - $e) ** (1 - $e));
1430
1431 9
        $theta = $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue();
1432 9
        $dE = $rho * sin($theta);
1433 9
        $dN = $rho * cos($theta);
1434
1435 9
        $easting = $falseEasting->asMetres()->getValue() + $dE;
1436 9
        if ($latitudeOrigin < 0) {
1437
            $northing = $falseNorthing->asMetres()->getValue() + $dN;
1438
        } else {
1439 9
            $northing = $falseNorthing->asMetres()->getValue() - $dN;
1440
        }
1441
1442 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1443
    }
1444
1445
    /**
1446
     * Polar Stereographic (variant B).
1447
     */
1448 9
    public function polarStereographicVariantB(
1449
        Projected $to,
1450
        Angle $latitudeOfStandardParallel,
1451
        Angle $longitudeOfOrigin,
1452
        Length $falseEasting,
1453
        Length $falseNorthing
1454
    ): ProjectedPoint {
1455 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1456 9
        $latitude = $this->latitude->asRadians()->getValue();
1457 9
        $firstStandardParallel = $latitudeOfStandardParallel->asRadians()->getValue();
1458 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1459 9
        $e = $ellipsoid->getEccentricity();
1460 9
        $e2 = $ellipsoid->getEccentricitySquared();
1461
1462 9
        if ($firstStandardParallel < 0) {
1463 9
            $tF = tan(M_PI / 4 + $firstStandardParallel / 2) / (((1 + $e * sin($firstStandardParallel)) / (1 - $e * sin($firstStandardParallel))) ** ($e / 2));
1464 9
            $t = tan(M_PI / 4 + $latitude / 2) / (((1 + $e * sin($latitude)) / (1 - $e * sin($latitude))) ** ($e / 2));
1465
        } else {
1466
            $tF = tan(M_PI / 4 - $firstStandardParallel / 2) * (((1 + $e * sin($firstStandardParallel)) / (1 - $e * sin($firstStandardParallel))) ** ($e / 2));
1467
            $t = tan(M_PI / 4 - $latitude / 2) * (((1 + $e * sin($latitude)) / (1 - $e * sin($latitude))) ** ($e / 2));
1468
        }
1469 9
        $mF = cos($firstStandardParallel) / sqrt(1 - $e2 * sin($firstStandardParallel) ** 2);
1470 9
        $kO = $mF * sqrt((1 + $e) ** (1 + $e) * (1 - $e) ** (1 - $e)) / (2 * $tF);
1471
1472 9
        $rho = 2 * $a * $kO * $t / sqrt((1 + $e) ** (1 + $e) * (1 - $e) ** (1 - $e));
1473
1474 9
        $theta = $this->normaliseLongitude($this->longitude->subtract($longitudeOfOrigin))->asRadians()->getValue();
1475 9
        $dE = $rho * sin($theta);
1476 9
        $dN = $rho * cos($theta);
1477
1478 9
        $easting = $falseEasting->asMetres()->getValue() + $dE;
1479 9
        if ($firstStandardParallel < 0) {
1480 9
            $northing = $falseNorthing->asMetres()->getValue() + $dN;
1481
        } else {
1482
            $northing = $falseNorthing->asMetres()->getValue() - $dN;
1483
        }
1484
1485 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1486
    }
1487
1488
    /**
1489
     * Polar Stereographic (variant C).
1490
     */
1491 9
    public function polarStereographicVariantC(
1492
        Projected $to,
1493
        Angle $latitudeOfStandardParallel,
1494
        Angle $longitudeOfOrigin,
1495
        Length $eastingAtFalseOrigin,
1496
        Length $northingAtFalseOrigin
1497
    ): ProjectedPoint {
1498 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1499 9
        $latitude = $this->latitude->asRadians()->getValue();
1500 9
        $firstStandardParallel = $latitudeOfStandardParallel->asRadians()->getValue();
1501 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1502 9
        $e = $ellipsoid->getEccentricity();
1503 9
        $e2 = $ellipsoid->getEccentricitySquared();
1504
1505 9
        if ($firstStandardParallel < 0) {
1506 9
            $tF = tan(M_PI / 4 + $firstStandardParallel / 2) / (((1 + $e * sin($firstStandardParallel)) / (1 - $e * sin($firstStandardParallel))) ** ($e / 2));
1507 9
            $t = tan(M_PI / 4 + $latitude / 2) / (((1 + $e * sin($latitude)) / (1 - $e * sin($latitude))) ** ($e / 2));
1508
        } else {
1509
            $tF = tan(M_PI / 4 - $firstStandardParallel / 2) * (((1 + $e * sin($firstStandardParallel)) / (1 - $e * sin($firstStandardParallel))) ** ($e / 2));
1510
            $t = tan(M_PI / 4 - $latitude / 2) * (((1 + $e * sin($latitude)) / (1 - $e * sin($latitude))) ** ($e / 2));
1511
        }
1512 9
        $mF = cos($firstStandardParallel) / sqrt(1 - $e2 * sin($firstStandardParallel) ** 2);
1513
1514 9
        $rhoF = $a * $mF;
1515 9
        $rho = $rhoF * $t / $tF;
1516
1517 9
        $theta = $this->normaliseLongitude($this->longitude->subtract($longitudeOfOrigin))->asRadians()->getValue();
1518 9
        $dE = $rho * sin($theta);
1519 9
        $dN = $rho * cos($theta);
1520
1521 9
        $easting = $eastingAtFalseOrigin->asMetres()->getValue() + $dE;
1522 9
        if ($firstStandardParallel < 0) {
1523 9
            $northing = $northingAtFalseOrigin->asMetres()->getValue() - $rhoF + $dN;
1524
        } else {
1525
            $northing = $northingAtFalseOrigin->asMetres()->getValue() + $rhoF - $dN;
1526
        }
1527
1528 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1529
    }
1530
1531
    /**
1532
     * Popular Visualisation Pseudo Mercator
1533
     * Applies spherical formulas to the ellipsoid. As such does not have the properties of a true Mercator projection.
1534
     */
1535 9
    public function popularVisualisationPseudoMercator(
1536
        Projected $to,
1537
        Angle $latitudeOfNaturalOrigin,
0 ignored issues
show
Unused Code introduced by
The parameter $latitudeOfNaturalOrigin is not used and could be removed. ( Ignorable by Annotation )

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

1537
        /** @scrutinizer ignore-unused */ Angle $latitudeOfNaturalOrigin,

This check looks for parameters that have been defined for a function or method, but which are not used in the method body.

Loading history...
1538
        Angle $longitudeOfNaturalOrigin,
1539
        Length $falseEasting,
1540
        Length $falseNorthing
1541
    ): ProjectedPoint {
1542 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1543 9
        $latitude = $this->latitude->asRadians()->getValue();
1544 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1545
1546 9
        $easting = $falseEasting->asMetres()->getValue() + $a * ($this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue());
1547 9
        $northing = $falseNorthing->asMetres()->getValue() + $a * log(tan(M_PI / 4 + $latitude / 2));
1548
1549 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1550
    }
1551
1552
    /**
1553
     * Mercator (variant A)
1554
     * Note that in these formulas the parameter latitude of natural origin (latO) is not used. However for this
1555
     * Mercator (variant A) method the EPSG dataset includes this parameter, which must have a value of zero, for
1556
     * completeness in CRS labelling.
1557
     */
1558 18
    public function mercatorVariantA(
1559
        Projected $to,
1560
        Angle $latitudeOfNaturalOrigin,
0 ignored issues
show
Unused Code introduced by
The parameter $latitudeOfNaturalOrigin is not used and could be removed. ( Ignorable by Annotation )

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

1560
        /** @scrutinizer ignore-unused */ Angle $latitudeOfNaturalOrigin,

This check looks for parameters that have been defined for a function or method, but which are not used in the method body.

Loading history...
1561
        Angle $longitudeOfNaturalOrigin,
1562
        Scale $scaleFactorAtNaturalOrigin,
1563
        Length $falseEasting,
1564
        Length $falseNorthing
1565
    ): ProjectedPoint {
1566 18
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1567 18
        $latitude = $this->latitude->asRadians()->getValue();
1568 18
        $kO = $scaleFactorAtNaturalOrigin->asUnity()->getValue();
1569
1570 18
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1571 18
        $e = $ellipsoid->getEccentricity();
1572
1573 18
        $easting = $falseEasting->asMetres()->getValue() + $a * $kO * $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue();
1574 18
        $northing = $falseNorthing->asMetres()->getValue() + $a * $kO * log(tan(M_PI / 4 + $latitude / 2) * ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2));
1575
1576 18
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1577
    }
1578
1579
    /**
1580
     * Mercator (variant B)
1581
     * Used for most nautical charts.
1582
     */
1583 9
    public function mercatorVariantB(
1584
        Projected $to,
1585
        Angle $latitudeOf1stStandardParallel,
1586
        Angle $longitudeOfNaturalOrigin,
1587
        Length $falseEasting,
1588
        Length $falseNorthing
1589
    ): ProjectedPoint {
1590 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1591 9
        $latitude = $this->latitude->asRadians()->getValue();
1592 9
        $firstStandardParallel = $latitudeOf1stStandardParallel->asRadians()->getValue();
1593 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1594 9
        $e = $ellipsoid->getEccentricity();
1595 9
        $e2 = $ellipsoid->getEccentricitySquared();
1596
1597 9
        $kO = cos($firstStandardParallel) / sqrt(1 - $e2 * sin($firstStandardParallel) ** 2);
1598
1599 9
        $easting = $falseEasting->asMetres()->getValue() + $a * $kO * $this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue();
1600 9
        $northing = $falseNorthing->asMetres()->getValue() + $a * $kO * log(tan(M_PI / 4 + $latitude / 2) * ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2));
1601
1602 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1603
    }
1604
1605
    /**
1606
     * Longitude rotation
1607
     * This transformation allows calculation of the longitude of a point in the target system by adding the parameter
1608
     * value to the longitude value of the point in the source system.
1609
     */
1610 9
    public function longitudeRotation(
1611
        Geographic $to,
1612
        Angle $longitudeOffset
1613
    ): self {
1614 9
        $newLongitude = $this->longitude->add($longitudeOffset);
1615
1616 9
        return static::create($this->latitude, $newLongitude, $this->height, $to, $this->epoch);
1617
    }
1618
1619
    /**
1620
     * Hotine Oblique Mercator (variant A).
1621
     */
1622 9
    public function obliqueMercatorHotineVariantA(
1623
        Projected $to,
1624
        Angle $latitudeOfProjectionCentre,
1625
        Angle $longitudeOfProjectionCentre,
1626
        Angle $azimuthOfInitialLine,
1627
        Angle $angleFromRectifiedToSkewGrid,
1628
        Scale $scaleFactorOnInitialLine,
1629
        Length $falseEasting,
1630
        Length $falseNorthing
1631
    ): ProjectedPoint {
1632 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1633 9
        $latitude = $this->latitude->asRadians()->getValue();
1634 9
        $longitude = $this->longitude->asRadians()->getValue();
1635 9
        $latC = $latitudeOfProjectionCentre->asRadians()->getValue();
1636 9
        $lonC = $longitudeOfProjectionCentre->asRadians()->getValue();
1637 9
        $alphaC = $azimuthOfInitialLine->asRadians()->getValue();
1638 9
        $kC = $scaleFactorOnInitialLine->asUnity()->getValue();
1639 9
        $gammaC = $angleFromRectifiedToSkewGrid->asRadians()->getValue();
1640 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1641 9
        $e = $ellipsoid->getEccentricity();
1642 9
        $e2 = $ellipsoid->getEccentricitySquared();
1643
1644 9
        $B = sqrt(1 + ($e2 * cos($latC) ** 4 / (1 - $e2)));
1645 9
        $A = $a * $B * $kC * sqrt(1 - $e2) / (1 - $e2 * sin($latC) ** 2);
1646 9
        $tO = tan(M_PI / 4 - $latC / 2) / ((1 - $e * sin($latC)) / (1 + $e * sin($latC))) ** ($e / 2);
1647 9
        $D = $B * sqrt((1 - $e2)) / (cos($latC) * sqrt(1 - $e2 * sin($latC) ** 2));
1648 9
        $DD = max(1, $D ** 2);
1649 9
        $F = $D + sqrt($DD - 1) * static::sign($latC);
1650 9
        $H = $F * ($tO) ** $B;
1651 9
        $G = ($F - 1 / $F) / 2;
1652 9
        $gammaO = self::asin(sin($alphaC) / $D);
1653 9
        $lonO = $lonC - (self::asin($G * tan($gammaO))) / $B;
1654
1655 9
        $t = tan(M_PI / 4 - $latitude / 2) / ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2);
1656 9
        $Q = $H / $t ** $B;
1657 9
        $S = ($Q - 1 / $Q) / 2;
1658 9
        $T = ($Q + 1 / $Q) / 2;
1659 9
        $V = sin($B * ($longitude - $lonO));
1660 9
        $U = (-$V * cos($gammaO) + $S * sin($gammaO)) / $T;
1661 9
        $v = $A * log((1 - $U) / (1 + $U)) / (2 * $B);
1662 9
        $u = $A * atan2(($S * cos($gammaO) + $V * sin($gammaO)), cos($B * ($longitude - $lonO))) / $B;
1663
1664 9
        $easting = $v * cos($gammaC) + $u * sin($gammaC) + $falseEasting->asMetres()->getValue();
1665 9
        $northing = $u * cos($gammaC) - $v * sin($gammaC) + $falseNorthing->asMetres()->getValue();
1666
1667 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1668
    }
1669
1670
    /**
1671
     * Hotine Oblique Mercator (variant B).
1672
     */
1673 9
    public function obliqueMercatorHotineVariantB(
1674
        Projected $to,
1675
        Angle $latitudeOfProjectionCentre,
1676
        Angle $longitudeOfProjectionCentre,
1677
        Angle $azimuthOfInitialLine,
1678
        Angle $angleFromRectifiedToSkewGrid,
1679
        Scale $scaleFactorOnInitialLine,
1680
        Length $eastingAtProjectionCentre,
1681
        Length $northingAtProjectionCentre
1682
    ): ProjectedPoint {
1683 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1684 9
        $latitude = $this->latitude->asRadians()->getValue();
1685 9
        $longitude = $this->longitude->asRadians()->getValue();
1686 9
        $latC = $latitudeOfProjectionCentre->asRadians()->getValue();
1687 9
        $lonC = $longitudeOfProjectionCentre->asRadians()->getValue();
1688 9
        $alphaC = $azimuthOfInitialLine->asRadians()->getValue();
1689 9
        $kC = $scaleFactorOnInitialLine->asUnity()->getValue();
1690 9
        $gammaC = $angleFromRectifiedToSkewGrid->asRadians()->getValue();
1691 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1692 9
        $e = $ellipsoid->getEccentricity();
1693 9
        $e2 = $ellipsoid->getEccentricitySquared();
1694
1695 9
        $B = sqrt(1 + ($e2 * cos($latC) ** 4 / (1 - $e2)));
1696 9
        $A = $a * $B * $kC * sqrt(1 - $e2) / (1 - $e2 * sin($latC) ** 2);
1697 9
        $tO = tan(M_PI / 4 - $latC / 2) / ((1 - $e * sin($latC)) / (1 + $e * sin($latC))) ** ($e / 2);
1698 9
        $D = $B * sqrt((1 - $e2)) / (cos($latC) * sqrt(1 - $e2 * sin($latC) ** 2));
1699 9
        $DD = max(1, $D ** 2);
1700 9
        $F = $D + sqrt($DD - 1) * static::sign($latC);
1701 9
        $H = $F * ($tO) ** $B;
1702 9
        $G = ($F - 1 / $F) / 2;
1703 9
        $gammaO = self::asin(sin($alphaC) / $D);
1704 9
        $lonO = $lonC - (self::asin($G * tan($gammaO))) / $B;
1705 9
        $vC = 0;
0 ignored issues
show
Unused Code introduced by
The assignment to $vC is dead and can be removed.
Loading history...
1706 9
        if ($alphaC === M_PI / 2) {
1707
            $uC = $A * ($lonC - $lonO);
1708
        } else {
1709 9
            $uC = ($A / $B) * atan2(sqrt($DD - 1), cos($alphaC)) * static::sign($latC);
1710
        }
1711
1712 9
        $t = tan(M_PI / 4 - $latitude / 2) / ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2);
1713 9
        $Q = $H / $t ** $B;
1714 9
        $S = ($Q - 1 / $Q) / 2;
1715 9
        $T = ($Q + 1 / $Q) / 2;
1716 9
        $V = sin($B * ($longitude - $lonO));
1717 9
        $U = (-$V * cos($gammaO) + $S * sin($gammaO)) / $T;
1718 9
        $v = $A * log((1 - $U) / (1 + $U)) / (2 * $B);
1719
1720 9
        if ($alphaC === M_PI / 2) {
1721
            if ($longitude === $lonC) {
1722
                $u = 0;
1723
            } else {
1724
                $u = ($A * atan(($S * cos($gammaO) + $V * sin($gammaO)) / cos($B * ($longitude - $lonO))) / $B) - (abs($uC) * static::sign($latC) * static::sign($lonC - $longitude));
1725
            }
1726
        } else {
1727 9
            $u = ($A * atan2(($S * cos($gammaO) + $V * sin($gammaO)), cos($B * ($longitude - $lonO))) / $B) - (abs($uC) * static::sign($latC));
1728
        }
1729
1730 9
        $easting = $v * cos($gammaC) + $u * sin($gammaC) + $eastingAtProjectionCentre->asMetres()->getValue();
1731 9
        $northing = $u * cos($gammaC) - $v * sin($gammaC) + $northingAtProjectionCentre->asMetres()->getValue();
1732
1733 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1734
    }
1735
1736
    /**
1737
     * Laborde Oblique Mercator.
1738
     */
1739 9
    public function obliqueMercatorLaborde(
1740
        Projected $to,
1741
        Angle $latitudeOfProjectionCentre,
1742
        Angle $longitudeOfProjectionCentre,
1743
        Angle $azimuthOfInitialLine,
1744
        Scale $scaleFactorOnInitialLine,
1745
        Length $falseEasting,
1746
        Length $falseNorthing
1747
    ): ProjectedPoint {
1748 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1749 9
        $latitude = $this->latitude->asRadians()->getValue();
1750 9
        $latC = $latitudeOfProjectionCentre->asRadians()->getValue();
1751 9
        $alphaC = $azimuthOfInitialLine->asRadians()->getValue();
1752 9
        $kC = $scaleFactorOnInitialLine->asUnity()->getValue();
1753 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1754 9
        $e = $ellipsoid->getEccentricity();
1755 9
        $e2 = $ellipsoid->getEccentricitySquared();
1756
1757 9
        $B = sqrt(1 + ($e2 * cos($latC) ** 4 / (1 - $e2)));
1758 9
        $latS = self::asin(sin($latC) / $B);
1759 9
        $R = $a * $kC * (sqrt(1 - $e2) / (1 - $e2 * sin($latC) ** 2));
1760 9
        $C = log(tan(M_PI / 4 + $latS / 2)) - $B * log(tan(M_PI / 4 + $latC / 2) * ((1 - $e * sin($latC)) / (1 + $e * sin($latC))) ** ($e / 2));
1761
1762 9
        $L = $B * $this->normaliseLongitude($this->longitude->subtract($longitudeOfProjectionCentre))->asRadians()->getValue();
1763 9
        $q = $C + $B * log(tan(M_PI / 4 + $latitude / 2) * ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2));
1764 9
        $P = 2 * atan(M_E ** $q) - M_PI / 2;
1765 9
        $U = cos($P) * cos($L) * cos($latS) + sin($P) * sin($latS);
1766 9
        $V = cos($P) * cos($L) * sin($latS) - sin($P) * cos($latS);
1767 9
        $W = cos($P) * sin($L);
1768 9
        $d = hypot($U, $V);
1769 9
        if ($d === 0.0) {
1770
            $LPrime = 0;
1771
            $PPrime = static::sign($W) * M_PI / 2;
1772
        } else {
1773 9
            $LPrime = 2 * atan($V / ($U + $d));
1774 9
            $PPrime = atan($W / $d);
1775
        }
1776 9
        $H = new ComplexNumber(-$LPrime, log(tan(M_PI / 4 + $PPrime / 2)));
1777 9
        $G = (new ComplexNumber(1 - cos(2 * $alphaC), sin(2 * $alphaC)))->divide(new ComplexNumber(12, 0));
1778
1779 9
        $easting = $falseEasting->asMetres()->getValue() + $R * $H->pow(3)->multiply($G)->add($H)->getImaginary();
1780 9
        $northing = $falseNorthing->asMetres()->getValue() + $R * $H->pow(3)->multiply($G)->add($H)->getReal();
1781
1782 9
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1783
    }
1784
1785
    /**
1786
     * Transverse Mercator.
1787
     */
1788 112
    public function transverseMercator(
1789
        Projected $to,
1790
        Angle $latitudeOfNaturalOrigin,
1791
        Angle $longitudeOfNaturalOrigin,
1792
        Scale $scaleFactorAtNaturalOrigin,
1793
        Length $falseEasting,
1794
        Length $falseNorthing
1795
    ): ProjectedPoint {
1796 112
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1797 112
        $latitude = $this->latitude->asRadians()->getValue();
1798 112
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
1799 112
        $kO = $scaleFactorAtNaturalOrigin->asUnity()->getValue();
1800 112
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1801 112
        $e = $ellipsoid->getEccentricity();
1802 112
        $f = $ellipsoid->getInverseFlattening();
1803
1804 112
        $n = $f / (2 - $f);
1805 112
        $B = ($a / (1 + $n)) * (1 + $n ** 2 / 4 + $n ** 4 / 64);
1806
1807 112
        $h1 = $n / 2 - (2 / 3) * $n ** 2 + (5 / 16) * $n ** 3 + (41 / 180) * $n ** 4;
1808 112
        $h2 = (13 / 48) * $n ** 2 - (3 / 5) * $n ** 3 + (557 / 1440) * $n ** 4;
1809 112
        $h3 = (61 / 240) * $n ** 3 - (103 / 140) * $n ** 4;
1810 112
        $h4 = (49561 / 161280) * $n ** 4;
1811
1812 112
        if ($latitudeOrigin === 0.0) {
0 ignored issues
show
introduced by
The condition $latitudeOrigin === 0.0 is always false.
Loading history...
1813 81
            $mO = 0;
1814 31
        } elseif ($latitudeOrigin === M_PI / 2) {
1815
            $mO = $B * M_PI / 2;
1816 31
        } elseif ($latitudeOrigin === -M_PI / 2) {
1817
            $mO = $B * -M_PI / 2;
1818
        } else {
1819 31
            $qO = asinh(tan($latitudeOrigin)) - ($e * atanh($e * sin($latitudeOrigin)));
1820 31
            $betaO = atan(sinh($qO));
1821 31
            $xiO0 = self::asin(sin($betaO));
1822 31
            $xiO1 = $h1 * sin(2 * $xiO0);
1823 31
            $xiO2 = $h2 * sin(4 * $xiO0);
1824 31
            $xiO3 = $h3 * sin(6 * $xiO0);
1825 31
            $xiO4 = $h4 * sin(8 * $xiO0);
1826 31
            $xiO = $xiO0 + $xiO1 + $xiO2 + $xiO3 + $xiO4;
1827 31
            $mO = $B * $xiO;
1828
        }
1829
1830 112
        $Q = asinh(tan($latitude)) - ($e * atanh($e * sin($latitude)));
1831 112
        $beta = atan(sinh($Q));
1832 112
        $eta0 = atanh(cos($beta) * sin($this->normaliseLongitude($this->longitude->subtract($longitudeOfNaturalOrigin))->asRadians()->getValue()));
1833 112
        $xi0 = self::asin(sin($beta) * cosh($eta0));
1834 112
        $xi1 = $h1 * sin(2 * $xi0) * cosh(2 * $eta0);
1835 112
        $eta1 = $h1 * cos(2 * $xi0) * sinh(2 * $eta0);
1836 112
        $xi2 = $h2 * sin(4 * $xi0) * cosh(4 * $eta0);
1837 112
        $eta2 = $h2 * cos(4 * $xi0) * sinh(4 * $eta0);
1838 112
        $xi3 = $h3 * sin(6 * $xi0) * cosh(6 * $eta0);
1839 112
        $eta3 = $h3 * cos(6 * $xi0) * sinh(6 * $eta0);
1840 112
        $xi4 = $h4 * sin(8 * $xi0) * cosh(8 * $eta0);
1841 112
        $eta4 = $h4 * cos(8 * $xi0) * sinh(8 * $eta0);
1842 112
        $xi = $xi0 + $xi1 + $xi2 + $xi3 + $xi4;
1843 112
        $eta = $eta0 + $eta1 + $eta2 + $eta3 + $eta4;
1844
1845 112
        $easting = $falseEasting->asMetres()->getValue() + $kO * $B * $eta;
1846 112
        $northing = $falseNorthing->asMetres()->getValue() + $kO * ($B * $xi - $mO);
1847
1848 112
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1849
    }
1850
1851
    /**
1852
     * Transverse Mercator Zoned Grid System
1853
     * If locations fall outwith the fixed zones the general Transverse Mercator method (code 9807) must be used for
1854
     * each zone.
1855
     */
1856 36
    public function transverseMercatorZonedGrid(
1857
        Projected $to,
1858
        Angle $latitudeOfNaturalOrigin,
1859
        Angle $initialLongitude,
1860
        Angle $zoneWidth,
1861
        Scale $scaleFactorAtNaturalOrigin,
1862
        Length $falseEasting,
1863
        Length $falseNorthing
1864
    ): ProjectedPoint {
1865 36
        $W = $zoneWidth->asDegrees()->getValue();
1866 36
        $Z = (int) ($this->longitude->subtract($initialLongitude)->asDegrees()->getValue() / $W) % (int) (360 / $W) + 1;
1867
1868 36
        $longitudeOrigin = $initialLongitude->add(new Degree($Z * $W - $W / 2));
1869 36
        $falseEasting = $falseEasting->add(new Metre($Z * 1000000));
1870
1871 36
        return $this->transverseMercator($to, $latitudeOfNaturalOrigin, $longitudeOrigin, $scaleFactorAtNaturalOrigin, $falseEasting, $falseNorthing);
1872
    }
1873
1874
    /**
1875
     * New Zealand Map Grid.
1876
     */
1877 27
    public function newZealandMapGrid(
1878
        Projected $to,
1879
        Angle $latitudeOfNaturalOrigin,
1880
        Angle $longitudeOfNaturalOrigin,
1881
        Length $falseEasting,
1882
        Length $falseNorthing
1883
    ): ProjectedPoint {
1884 27
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1885 27
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1886
1887 27
        $deltaLatitudeToOrigin = Angle::convert($this->latitude->subtract($latitudeOfNaturalOrigin), Angle::EPSG_ARC_SECOND)->getValue();
1888 27
        $deltaLongitudeToOrigin = $this->longitude->subtract($longitudeOfNaturalOrigin)->asRadians();
1889
1890 27
        $deltaPsi = 0;
1891 27
        $deltaPsi += 0.6399175073 * ($deltaLatitudeToOrigin * 0.00001) ** 1;
1892 27
        $deltaPsi += -0.1358797613 * ($deltaLatitudeToOrigin * 0.00001) ** 2;
1893 27
        $deltaPsi += 0.063294409 * ($deltaLatitudeToOrigin * 0.00001) ** 3;
1894 27
        $deltaPsi += -0.02526853 * ($deltaLatitudeToOrigin * 0.00001) ** 4;
1895 27
        $deltaPsi += 0.0117879 * ($deltaLatitudeToOrigin * 0.00001) ** 5;
1896 27
        $deltaPsi += -0.0055161 * ($deltaLatitudeToOrigin * 0.00001) ** 6;
1897 27
        $deltaPsi += 0.0026906 * ($deltaLatitudeToOrigin * 0.00001) ** 7;
1898 27
        $deltaPsi += -0.001333 * ($deltaLatitudeToOrigin * 0.00001) ** 8;
1899 27
        $deltaPsi += 0.00067 * ($deltaLatitudeToOrigin * 0.00001) ** 9;
1900 27
        $deltaPsi += -0.00034 * ($deltaLatitudeToOrigin * 0.00001) ** 10;
1901
1902 27
        $zeta = new ComplexNumber($deltaPsi, $deltaLongitudeToOrigin->getValue());
1903
1904 27
        $B1 = new ComplexNumber(0.7557853228, 0.0);
1905 27
        $B2 = new ComplexNumber(0.249204646, 0.003371507);
1906 27
        $B3 = new ComplexNumber(-0.001541739, 0.041058560);
1907 27
        $B4 = new ComplexNumber(-0.10162907, 0.01727609);
1908 27
        $B5 = new ComplexNumber(-0.26623489, -0.36249218);
1909 27
        $B6 = new ComplexNumber(-0.6870983, -1.1651967);
1910 27
        $z = new ComplexNumber(0, 0);
1911 27
        $z = $z->add($B1->multiply($zeta->pow(1)));
1912 27
        $z = $z->add($B2->multiply($zeta->pow(2)));
1913 27
        $z = $z->add($B3->multiply($zeta->pow(3)));
1914 27
        $z = $z->add($B4->multiply($zeta->pow(4)));
1915 27
        $z = $z->add($B5->multiply($zeta->pow(5)));
1916 27
        $z = $z->add($B6->multiply($zeta->pow(6)));
1917
1918 27
        $easting = $falseEasting->asMetres()->getValue() + $z->getImaginary() * $a;
1919 27
        $northing = $falseNorthing->asMetres()->getValue() + $z->getReal() * $a;
1920
1921 27
        return ProjectedPoint::create(new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $to, $this->epoch);
1922
    }
1923
1924
    /**
1925
     * Madrid to ED50 polynomial.
1926
     */
1927 9
    public function madridToED50Polynomial(
1928
        Geographic2D $to,
1929
        Scale $A0,
1930
        Scale $A1,
1931
        Scale $A2,
1932
        Scale $A3,
1933
        Angle $B00,
1934
        Scale $B0,
1935
        Scale $B1,
1936
        Scale $B2,
1937
        Scale $B3
1938
    ): self {
1939 9
        $dLatitude = new ArcSecond($A0->add($A1->multiply($this->latitude->getValue()))->add($A2->multiply($this->longitude->getValue()))->add($A3->multiply($this->height ? $this->height->getValue() : 0))->getValue());
1940 9
        $dLongitude = $B00->add(new ArcSecond($B0->add($B1->multiply($this->latitude->getValue()))->add($B2->multiply($this->longitude->getValue()))->add($B3->multiply($this->height ? $this->height->getValue() : 0))->getValue()));
1941
1942 9
        return self::create($this->latitude->add($dLatitude), $this->longitude->add($dLongitude), null, $to, $this->epoch);
1943
    }
1944
1945
    /**
1946
     * Geographic3D to 2D conversion.
1947
     */
1948 46
    public function threeDToTwoD(
1949
        Geographic $to
1950
    ): self {
1951 46
        if ($to instanceof Geographic2D) {
1952 46
            return static::create($this->latitude, $this->longitude, null, $to, $this->epoch);
1953
        }
1954
1955
        return static::create($this->latitude, $this->longitude, new Metre(0), $to, $this->epoch);
1956
    }
1957
1958
    /**
1959
     * Geographic2D offsets.
1960
     * This transformation allows calculation of coordinates in the target system by adding the parameter value to the
1961
     * coordinate values of the point in the source system.
1962
     */
1963 9
    public function geographic2DOffsets(
1964
        Geographic $to,
1965
        Angle $latitudeOffset,
1966
        Angle $longitudeOffset
1967
    ): self {
1968 9
        $toLatitude = $this->latitude->add($latitudeOffset);
1969 9
        $toLongitude = $this->longitude->add($longitudeOffset);
1970
1971 9
        return static::create($toLatitude, $toLongitude, null, $to, $this->epoch);
1972
    }
1973
1974
    /*
1975
     * Geographic2D with Height Offsets.
1976
     * This transformation allows calculation of coordinates in the target system by adding the parameter value to the
1977
     * coordinate values of the point in the source system.
1978
     */
1979
    public function geographic2DWithHeightOffsets(
1980
        Compound $to,
1981
        Angle $latitudeOffset,
1982
        Angle $longitudeOffset,
1983
        Length $geoidUndulation
1984
    ): CompoundPoint {
1985
        $toLatitude = $this->latitude->add($latitudeOffset);
1986
        $toLongitude = $this->longitude->add($longitudeOffset);
1987
        $toHeight = $this->height->add($geoidUndulation);
0 ignored issues
show
Bug introduced by
The method add() does not exist on null. ( Ignorable by Annotation )

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

1987
        /** @scrutinizer ignore-call */ 
1988
        $toHeight = $this->height->add($geoidUndulation);

This check looks for calls to methods that do not seem to exist on a given type. It looks for the method on the type itself as well as in inherited classes or implemented interfaces.

This is most likely a typographical error or the method has been renamed.

Loading history...
1988
1989
        $horizontal = static::create($toLatitude, $toLongitude, null, $to->getHorizontal(), $this->epoch);
1990
        $vertical = VerticalPoint::create($toHeight, $to->getVertical(), $this->epoch);
0 ignored issues
show
Bug introduced by
It seems like $toHeight can also be of type null; however, parameter $height of PHPCoord\VerticalPoint::create() does only seem to accept PHPCoord\UnitOfMeasure\Length\Length, maybe add an additional type check? ( Ignorable by Annotation )

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

1990
        $vertical = VerticalPoint::create(/** @scrutinizer ignore-type */ $toHeight, $to->getVertical(), $this->epoch);
Loading history...
1991
1992
        return CompoundPoint::create($horizontal, $vertical, $to, $this->epoch);
1993
    }
1994
1995
    /**
1996
     * General polynomial.
1997
     * @param Coefficient[] $powerCoefficients
1998
     */
1999 18
    public function generalPolynomial(
2000
        Geographic $to,
2001
        Angle $ordinate1OfEvaluationPointInSourceCRS,
2002
        Angle $ordinate2OfEvaluationPointInSourceCRS,
2003
        Angle $ordinate1OfEvaluationPointInTargetCRS,
2004
        Angle $ordinate2OfEvaluationPointInTargetCRS,
2005
        Scale $scalingFactorForSourceCRSCoordDifferences,
2006
        Scale $scalingFactorForTargetCRSCoordDifferences,
2007
        Scale $A0,
2008
        Scale $B0,
2009
        array $powerCoefficients
2010
    ): self {
2011 18
        $xs = $this->latitude->getValue();
2012 18
        $ys = $this->longitude->getValue();
2013
2014 18
        $t = $this->generalPolynomialUnitless(
2015 18
            $xs,
2016
            $ys,
2017
            $ordinate1OfEvaluationPointInSourceCRS,
2018
            $ordinate2OfEvaluationPointInSourceCRS,
2019
            $ordinate1OfEvaluationPointInTargetCRS,
2020
            $ordinate2OfEvaluationPointInTargetCRS,
2021
            $scalingFactorForSourceCRSCoordDifferences,
2022
            $scalingFactorForTargetCRSCoordDifferences,
2023
            $A0,
2024
            $B0,
2025
            $powerCoefficients
2026
        );
2027
2028 18
        $xtUnit = $to->getCoordinateSystem()->getAxes()[0]->getUnitOfMeasureId();
2029 18
        $ytUnit = $to->getCoordinateSystem()->getAxes()[1]->getUnitOfMeasureId();
2030
2031 18
        return static::create(
2032 18
            Angle::makeUnit($t['xt'], $xtUnit),
2033 18
            Angle::makeUnit($t['yt'], $ytUnit),
2034 18
            $this->height,
2035
            $to,
2036 18
            $this->epoch
2037
        );
2038
    }
2039
2040
    /**
2041
     * Reversible polynomial.
2042
     * @param Coefficient[] $powerCoefficients
2043
     */
2044 36
    public function reversiblePolynomial(
2045
        Geographic $to,
2046
        Angle $ordinate1OfEvaluationPoint,
2047
        Angle $ordinate2OfEvaluationPoint,
2048
        Scale $scalingFactorForCoordDifferences,
2049
        Scale $A0,
2050
        Scale $B0,
2051
        $powerCoefficients
2052
    ): self {
2053 36
        $xs = $this->latitude->getValue();
2054 36
        $ys = $this->longitude->getValue();
2055
2056 36
        $t = $this->reversiblePolynomialUnitless(
2057 36
            $xs,
2058
            $ys,
2059
            $ordinate1OfEvaluationPoint,
2060
            $ordinate2OfEvaluationPoint,
2061
            $scalingFactorForCoordDifferences,
2062
            $A0,
2063
            $B0,
2064
            $powerCoefficients
2065
        );
2066
2067 36
        $xtUnit = $to->getCoordinateSystem()->getAxes()[0]->getUnitOfMeasureId();
2068 36
        $ytUnit = $to->getCoordinateSystem()->getAxes()[1]->getUnitOfMeasureId();
2069
2070 36
        return static::create(
2071 36
            Angle::makeUnit($t['xt'], $xtUnit),
2072 36
            Angle::makeUnit($t['yt'], $ytUnit),
2073 36
            $this->height,
2074
            $to,
2075 36
            $this->epoch
2076
        );
2077
    }
2078
2079
    /**
2080
     * Axis Order Reversal.
2081
     */
2082
    public function axisReversal(
2083
        Geographic $to
2084
    ): self {
2085
        // axes are read in from the CRS, this is a book-keeping adjustment only
2086
        return static::create($this->latitude, $this->longitude, $this->height, $to, $this->epoch);
2087
    }
2088
2089
    /**
2090
     * Ordnance Survey National Transformation
2091
     * Geodetic transformation between ETRS89 (or WGS 84) and OSGB36 / National Grid.  Uses ETRS89 / National Grid as
2092
     * an intermediate coordinate system for bi-linear interpolation of gridded grid coordinate differences.
2093
     */
2094 2
    public function OSTN15(
2095
        Projected $to,
0 ignored issues
show
Unused Code introduced by
The parameter $to is not used and could be removed. ( Ignorable by Annotation )

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

2095
        /** @scrutinizer ignore-unused */ Projected $to,

This check looks for parameters that have been defined for a function or method, but which are not used in the method body.

Loading history...
2096
        OSTNOSGM15Grid $eastingAndNorthingDifferenceFile
2097
    ): ProjectedPoint {
2098 2
        $osgb36NationalGrid = Projected::fromSRID(Projected::EPSG_OSGB36_BRITISH_NATIONAL_GRID);
2099 2
        $etrs89NationalGrid = new Projected(
2100 2
            'ETRS89 / National Grid',
2101 2
            Cartesian::fromSRID(Cartesian::EPSG_2D_AXES_EASTING_NORTHING_E_N_ORIENTATIONS_EAST_NORTH_UOM_M),
2102 2
            Datum::fromSRID(Datum::EPSG_EUROPEAN_TERRESTRIAL_REFERENCE_SYSTEM_1989_ENSEMBLE),
2103 2
            $osgb36NationalGrid->getBoundingArea()
2104
        );
2105
2106 2
        $projected = $this->transverseMercator($etrs89NationalGrid, new Degree(49), new Degree(-2), new Unity(0.9996012717), new Metre(400000), new Metre(-100000));
2107
2108 2
        return $eastingAndNorthingDifferenceFile->applyForwardHorizontalAdjustment($projected);
2109
    }
2110
2111
    /**
2112
     * Geog3D to Geog2D+GravityRelatedHeight (OSGM-GB).
2113
     * Uses ETRS89 / National Grid as an intermediate coordinate system for bi-linear interpolation of gridded grid
2114
     * coordinate differences.
2115
     */
2116 1
    public function geographic3DTo2DPlusGravityHeightOSGM15(
2117
        Compound $to,
2118
        OSTNOSGM15Grid $geoidHeightCorrectionModelFile
2119
    ): CompoundPoint {
2120 1
        $osgb36NationalGrid = Projected::fromSRID(Projected::EPSG_OSGB36_BRITISH_NATIONAL_GRID);
2121 1
        $etrs89NationalGrid = new Projected(
2122 1
            'ETRS89 / National Grid',
2123 1
            Cartesian::fromSRID(Cartesian::EPSG_2D_AXES_EASTING_NORTHING_E_N_ORIENTATIONS_EAST_NORTH_UOM_M),
2124 1
            Datum::fromSRID(Datum::EPSG_EUROPEAN_TERRESTRIAL_REFERENCE_SYSTEM_1989_ENSEMBLE),
2125 1
            $osgb36NationalGrid->getBoundingArea()
2126
        );
2127
2128 1
        $projected = $this->transverseMercator($etrs89NationalGrid, new Degree(49), new Degree(-2), new Unity(0.9996012717), new Metre(400000), new Metre(-100000));
2129
2130 1
        $horizontalPoint = self::create(
2131 1
            $this->latitude,
2132 1
            $this->longitude,
2133 1
            null,
2134 1
            $to->getHorizontal(),
2135 1
            $this->getCoordinateEpoch()
2136
        );
2137
2138 1
        $verticalPoint = VerticalPoint::create(
2139 1
            $this->height->subtract($geoidHeightCorrectionModelFile->getHeightAdjustment($projected)),
0 ignored issues
show
Bug introduced by
It seems like $this->height->subtract(...Adjustment($projected)) can also be of type null; however, parameter $height of PHPCoord\VerticalPoint::create() does only seem to accept PHPCoord\UnitOfMeasure\Length\Length, maybe add an additional type check? ( Ignorable by Annotation )

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

2139
            /** @scrutinizer ignore-type */ $this->height->subtract($geoidHeightCorrectionModelFile->getHeightAdjustment($projected)),
Loading history...
2140 1
            $to->getVertical(),
2141 1
            $this->getCoordinateEpoch()
2142
        );
2143
2144 1
        return CompoundPoint::create(
2145 1
            $horizontalPoint,
2146 1
            $verticalPoint,
2147 1
            $to,
2148 1
            $this->getCoordinateEpoch()
2149
        );
2150
    }
2151
2152
    /**
2153
     * Geographic3D to GravityRelatedHeight (OSGM-GB).
2154
     * Uses ETRS89 / National Grid as an intermediate coordinate system for bi-linear interpolation of gridded grid
2155
     * coordinate differences.
2156
     */
2157 1
    public function geographic3DToGravityHeightOSGM15(
2158
        Vertical $to,
2159
        OSTNOSGM15Grid $geoidHeightCorrectionModelFile
2160
    ): VerticalPoint {
2161 1
        $osgb36NationalGrid = Projected::fromSRID(Projected::EPSG_OSGB36_BRITISH_NATIONAL_GRID);
2162 1
        $etrs89NationalGrid = new Projected(
2163 1
            'ETRS89 / National Grid',
2164 1
            Cartesian::fromSRID(Cartesian::EPSG_2D_AXES_EASTING_NORTHING_E_N_ORIENTATIONS_EAST_NORTH_UOM_M),
2165 1
            Datum::fromSRID(Datum::EPSG_EUROPEAN_TERRESTRIAL_REFERENCE_SYSTEM_1989_ENSEMBLE),
2166 1
            $osgb36NationalGrid->getBoundingArea()
2167
        );
2168
2169 1
        $projected = $this->transverseMercator($etrs89NationalGrid, new Degree(49), new Degree(-2), new Unity(0.9996012717), new Metre(400000), new Metre(-100000));
2170
2171 1
        return VerticalPoint::create(
2172 1
            $this->height->subtract($geoidHeightCorrectionModelFile->getHeightAdjustment($projected)),
0 ignored issues
show
Bug introduced by
It seems like $this->height->subtract(...Adjustment($projected)) can also be of type null; however, parameter $height of PHPCoord\VerticalPoint::create() does only seem to accept PHPCoord\UnitOfMeasure\Length\Length, maybe add an additional type check? ( Ignorable by Annotation )

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

2172
            /** @scrutinizer ignore-type */ $this->height->subtract($geoidHeightCorrectionModelFile->getHeightAdjustment($projected)),
Loading history...
2173 1
            $to,
2174 1
            $this->getCoordinateEpoch()
2175
        );
2176
    }
2177
2178
    /**
2179
     * Geog3D to Geog2D+GravityRelatedHeight.
2180
     */
2181 2
    public function geographic3DTo2DPlusGravityHeightFromGrid(
2182
        Compound $to,
2183
        GeographicGeoidHeightGrid $geoidHeightCorrectionModelFile
2184
    ): CompoundPoint {
2185 2
        $horizontalPoint = self::create(
2186 2
            $this->latitude,
2187 2
            $this->longitude,
2188 2
            null,
2189 2
            $to->getHorizontal(),
2190 2
            $this->getCoordinateEpoch()
2191
        );
2192
2193 2
        $verticalPoint = VerticalPoint::create(
2194 2
            $this->height->subtract($geoidHeightCorrectionModelFile->getHeightAdjustment($this)),
0 ignored issues
show
Bug introduced by
It seems like $this->height->subtract(...eightAdjustment($this)) can also be of type null; however, parameter $height of PHPCoord\VerticalPoint::create() does only seem to accept PHPCoord\UnitOfMeasure\Length\Length, maybe add an additional type check? ( Ignorable by Annotation )

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

2194
            /** @scrutinizer ignore-type */ $this->height->subtract($geoidHeightCorrectionModelFile->getHeightAdjustment($this)),
Loading history...
2195 2
            $to->getVertical(),
2196 2
            $this->getCoordinateEpoch()
2197
        );
2198
2199 2
        return CompoundPoint::create(
2200 2
            $horizontalPoint,
2201 2
            $verticalPoint,
2202 2
            $to,
2203 2
            $this->getCoordinateEpoch()
2204
        );
2205
    }
2206
2207
    /**
2208
     * Geographic3D to GravityRelatedHeight.
2209
     */
2210 2
    public function geographic3DToGravityHeightFromGrid(
2211
        Vertical $to,
2212
        GeographicGeoidHeightGrid $geoidHeightCorrectionModelFile
2213
    ): VerticalPoint {
2214 2
        return VerticalPoint::create(
2215 2
            $this->height->subtract($geoidHeightCorrectionModelFile->getHeightAdjustment($this)),
0 ignored issues
show
Bug introduced by
It seems like $this->height->subtract(...eightAdjustment($this)) can also be of type null; however, parameter $height of PHPCoord\VerticalPoint::create() does only seem to accept PHPCoord\UnitOfMeasure\Length\Length, maybe add an additional type check? ( Ignorable by Annotation )

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

2215
            /** @scrutinizer ignore-type */ $this->height->subtract($geoidHeightCorrectionModelFile->getHeightAdjustment($this)),
Loading history...
2216 2
            $to,
2217 2
            $this->getCoordinateEpoch()
2218
        );
2219
    }
2220
2221
    /**
2222
     * NADCON5.
2223
     * @internal just a wrapper
2224
     */
2225 7
    public function offsetsFromGridNADCON5(
2226
        Geographic $to,
2227
        NADCON5Grid $latitudeDifferenceFile,
2228
        NADCON5Grid $longitudeDifferenceFile,
2229
        ?NADCON5Grid $ellipsoidalHeightDifferenceFile,
2230
        bool $inReverse
2231
    ): self {
2232 7
        $aggregation = new NADCON5Grids($longitudeDifferenceFile, $latitudeDifferenceFile, $ellipsoidalHeightDifferenceFile);
2233
2234 7
        return $this->offsetsFromGrid($to, $aggregation, $inReverse);
2235
    }
2236
2237
    /**
2238
     * Geographic offsets from grid.
2239
     */
2240 12
    public function offsetsFromGrid(
2241
        Geographic $to,
2242
        GeographicGrid $offsetsFile,
2243
        bool $inReverse
2244
    ): self {
2245 12
        if (!$inReverse) {
2246 7
            return $offsetsFile->applyForwardAdjustment($this, $to);
2247
        }
2248
2249 5
        return $offsetsFile->applyReverseAdjustment($this, $to);
2250
    }
2251
2252 353
    public function asGeographicValue(): GeographicValue
2253
    {
2254 353
        return new GeographicValue($this->latitude, $this->longitude, $this->height, $this->crs->getDatum());
2255
    }
2256
2257 18
    public function asUTMPoint(): UTMPoint
2258
    {
2259 18
        $hemisphere = $this->getLatitude()->asDegrees()->getValue() >= 0 ? UTMPoint::HEMISPHERE_NORTH : UTMPoint::HEMISPHERE_SOUTH;
2260 18
        $latitudeOfNaturalOrigin = new Degree(0);
2261 18
        $initialLongitude = new Degree(-180);
2262 18
        $scaleFactorAtNaturalOrigin = new Unity(0.9996);
2263 18
        $falseEasting = new Metre(500000);
2264 18
        $falseNorthing = $hemisphere === UTMPoint::HEMISPHERE_NORTH ? new Metre(0) : new Metre(10000000);
2265 18
        $Z = (int) ($this->longitude->subtract($initialLongitude)->asDegrees()->getValue() / 6) % (360 / 6) + 1;
2266 18
        $longitudeOrigin = $initialLongitude->add(new Degree($Z * 6 - 3));
2267
2268 18
        $projectedCRS = new Projected(
2269 18
            'UTM/' . $this->crs->getSRID(),
2270 18
            Cartesian::fromSRID(Cartesian::EPSG_2D_AXES_EASTING_NORTHING_E_N_ORIENTATIONS_EAST_NORTH_UOM_M),
2271 18
            $this->crs->getDatum(),
2272 18
            BoundingArea::createWorld() // this is a dummy CRS for the transform only, details don't matter
2273
        );
2274
2275 18
        $asProjected = $this->transverseMercator($projectedCRS, $latitudeOfNaturalOrigin, $longitudeOrigin, $scaleFactorAtNaturalOrigin, $falseEasting, $falseNorthing);
2276
2277 18
        return new UTMPoint($asProjected->getEasting(), $asProjected->getNorthing(), $Z, $hemisphere, $this->crs, $this->epoch);
2278
    }
2279
}
2280