Passed
Push — extents ( 6d8774...01b6e3 )
by Doug
61:12
created

ProjectedPoint::similarityTransformation()   A

Complexity

Conditions 2
Paths 2

Size

Total Lines 24
Code Lines 13

Duplication

Lines 0
Ratio 0 %

Code Coverage

Tests 15
CRAP Score 2

Importance

Changes 0
Metric Value
cc 2
eloc 13
nc 2
nop 6
dl 0
loc 24
ccs 15
cts 15
cp 1
crap 2
rs 9.8333
c 0
b 0
f 0
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 function count;
19
use DateTime;
20
use DateTimeImmutable;
21
use DateTimeInterface;
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 const M_PI_2;
29
use function max;
30
use PHPCoord\CoordinateOperation\AutoConversion;
31
use PHPCoord\CoordinateOperation\ComplexNumber;
32
use PHPCoord\CoordinateOperation\ConvertiblePoint;
33
use PHPCoord\CoordinateOperation\OSTNOSGM15Grid;
34
use PHPCoord\CoordinateReferenceSystem\CoordinateReferenceSystem;
35
use PHPCoord\CoordinateReferenceSystem\Geographic2D;
36
use PHPCoord\CoordinateReferenceSystem\Geographic3D;
37
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...
38
use PHPCoord\CoordinateSystem\Axis;
39
use PHPCoord\CoordinateSystem\Cartesian;
40
use PHPCoord\Exception\InvalidAxesException;
41
use PHPCoord\Exception\InvalidCoordinateReferenceSystemException;
42
use PHPCoord\Exception\UnknownAxisException;
43
use PHPCoord\UnitOfMeasure\Angle\Angle;
44
use PHPCoord\UnitOfMeasure\Angle\ArcSecond;
45
use PHPCoord\UnitOfMeasure\Angle\Degree;
46
use PHPCoord\UnitOfMeasure\Angle\Radian;
47
use PHPCoord\UnitOfMeasure\Length\Length;
48
use PHPCoord\UnitOfMeasure\Length\Metre;
49
use PHPCoord\UnitOfMeasure\Scale\Coefficient;
50
use PHPCoord\UnitOfMeasure\Scale\Scale;
51
use PHPCoord\UnitOfMeasure\Scale\Unity;
52
use function sin;
53
use function sinh;
54
use function sqrt;
55
use function substr;
56
use function tan;
57
use function tanh;
58
59
/**
60
 * Coordinate representing a point on a map projection.
61
 */
62
class ProjectedPoint extends Point implements ConvertiblePoint
63
{
64
    use AutoConversion {
65
        convert as protected autoConvert;
66
    }
67
68
    /**
69
     * Easting.
70
     */
71
    protected Length $easting;
72
73
    /**
74
     * Northing.
75
     */
76
    protected Length $northing;
77
78
    /**
79
     * Westing.
80
     */
81
    protected Length $westing;
82
83
    /**
84
     * Southing.
85
     */
86
    protected Length $southing;
87
88
    /**
89
     * Height.
90
     */
91
    protected ?Length $height;
92
93
    /**
94 46487
     * Coordinate reference system.
95
     */
96 46487
    protected Projected $crs;
97
98 46487
    /**
99 46487
     * Coordinate epoch (date for which the specified coordinates represented this point).
100 46487
     */
101 46487
    protected ?DateTimeImmutable $epoch;
102
103 46487
    protected function __construct(Projected $crs, ?Length $easting, ?Length $northing, ?Length $westing, ?Length $southing, ?DateTimeInterface $epoch, ?Length $height)
104 45893
    {
105 45893
        if (count($crs->getCoordinateSystem()->getAxes()) === 2 && $height !== null) {
106 612
            throw new InvalidCoordinateReferenceSystemException('A 2D projected point must not include a height');
107 603
        }
108 603
109
        if (count($crs->getCoordinateSystem()->getAxes()) === 3 && $height === null) {
110 9
            throw new InvalidCoordinateReferenceSystemException('A 3D projected point must include a height, none given');
111
        }
112
113 46478
        $this->crs = $crs;
114 46055
115 46055
        $eastingAxis = $this->crs->getCoordinateSystem()->getAxisByName(Axis::EASTING);
116 432
        $westingAxis = $this->crs->getCoordinateSystem()->getAxisByName(Axis::WESTING);
117 423
        $northingAxis = $this->crs->getCoordinateSystem()->getAxisByName(Axis::NORTHING);
118 423
        $southingAxis = $this->crs->getCoordinateSystem()->getAxisByName(Axis::SOUTHING);
119
120 9
        if ($easting && $eastingAxis) {
121
            $this->easting = $easting::convert($easting, $eastingAxis->getUnitOfMeasureId());
122
            $this->westing = $this->easting->multiply(-1);
123 46469
        } elseif ($westing && $westingAxis) {
124 738
            $this->westing = $westing::convert($westing, $westingAxis->getUnitOfMeasureId());
125
            $this->easting = $this->westing->multiply(-1);
126 46469
        } else {
127 46469
            throw new InvalidAxesException($crs->getCoordinateSystem()->getAxes());
128
        }
129 46262
130
        if ($northing && $northingAxis) {
131 46262
            $this->northing = $northing::convert($northing, $northingAxis->getUnitOfMeasureId());
132 102
            $this->southing = $this->northing->multiply(-1);
133
        } elseif ($southing && $southingAxis) {
134
            $this->southing = $southing::convert($southing, $southingAxis->getUnitOfMeasureId());
135 46201
            $this->northing = $this->southing->multiply(-1);
136 18
        } else {
137
            throw new InvalidAxesException($crs->getCoordinateSystem()->getAxes());
138
        }
139 46183
140 27
        if ($epoch instanceof DateTime) {
141
            $epoch = DateTimeImmutable::createFromMutable($epoch);
142
        }
143 46157
        $this->epoch = $epoch;
144
145
        $this->height = $height;
146 580
    }
147
148 580
    public static function create(Projected $crs, ?Length $easting, ?Length $northing, ?Length $westing, ?Length $southing, ?DateTimeInterface $epoch = null, ?Length $height = null): self
149
    {
150
        return match ($crs->getSRID()) {
151 36
            Projected::EPSG_OSGB36_BRITISH_NATIONAL_GRID => new BritishNationalGridPoint($easting, $northing, $epoch),
0 ignored issues
show
Bug introduced by
It seems like $easting can also be of type null; however, parameter $easting of PHPCoord\BritishNationalGridPoint::__construct() 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

151
            Projected::EPSG_OSGB36_BRITISH_NATIONAL_GRID => new BritishNationalGridPoint(/** @scrutinizer ignore-type */ $easting, $northing, $epoch),
Loading history...
Bug introduced by
It seems like $northing can also be of type null; however, parameter $northing of PHPCoord\BritishNationalGridPoint::__construct() 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

151
            Projected::EPSG_OSGB36_BRITISH_NATIONAL_GRID => new BritishNationalGridPoint($easting, /** @scrutinizer ignore-type */ $northing, $epoch),
Loading history...
152
            Projected::EPSG_TM75_IRISH_GRID => new IrishGridPoint($easting, $northing, $epoch),
0 ignored issues
show
Bug introduced by
It seems like $northing can also be of type null; however, parameter $northing of PHPCoord\IrishGridPoint::__construct() 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

152
            Projected::EPSG_TM75_IRISH_GRID => new IrishGridPoint($easting, /** @scrutinizer ignore-type */ $northing, $epoch),
Loading history...
Bug introduced by
It seems like $easting can also be of type null; however, parameter $easting of PHPCoord\IrishGridPoint::__construct() 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

152
            Projected::EPSG_TM75_IRISH_GRID => new IrishGridPoint(/** @scrutinizer ignore-type */ $easting, $northing, $epoch),
Loading history...
153 36
            Projected::EPSG_IRENET95_IRISH_TRANSVERSE_MERCATOR => new IrishTransverseMercatorPoint($easting, $northing, $epoch),
0 ignored issues
show
Bug introduced by
It seems like $easting can also be of type null; however, parameter $easting of PHPCoord\IrishTransverse...torPoint::__construct() 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

153
            Projected::EPSG_IRENET95_IRISH_TRANSVERSE_MERCATOR => new IrishTransverseMercatorPoint(/** @scrutinizer ignore-type */ $easting, $northing, $epoch),
Loading history...
Bug introduced by
It seems like $northing can also be of type null; however, parameter $northing of PHPCoord\IrishTransverse...torPoint::__construct() 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

153
            Projected::EPSG_IRENET95_IRISH_TRANSVERSE_MERCATOR => new IrishTransverseMercatorPoint($easting, /** @scrutinizer ignore-type */ $northing, $epoch),
Loading history...
154
            default => new static($crs, $easting, $northing, $westing, $southing, $epoch, $height),
155
        };
156 63
    }
157
158 63
    public static function createFromEastingNorthing(Projected $crs, Length $easting, Length $northing, ?DateTimeInterface $epoch = null, ?Length $height = null): self
159
    {
160
        return static::create($crs, $easting, $northing, null, null, $epoch, $height);
161 1239
    }
162
163 1239
    public static function createFromWestingNorthing(Projected $crs, Length $westing, Length $northing, ?DateTimeInterface $epoch = null, ?Length $height = null): self
164
    {
165
        return static::create($crs, null, $northing, $westing, null, $epoch, $height);
166 1257
    }
167
168 1257
    public static function createFromWestingSouthing(Projected $crs, Length $westing, Length $southing, ?DateTimeInterface $epoch = null, ?Length $height = null): self
169
    {
170
        return static::create($crs, null, null, $westing, $southing, $epoch, $height);
171 126
    }
172
173 126
    public function getEasting(): Length
174
    {
175
        return $this->easting;
176 108
    }
177
178 108
    public function getNorthing(): Length
179
    {
180
        return $this->northing;
181 45496
    }
182
183 45496
    public function getWesting(): Length
184
    {
185
        return $this->westing;
186 36
    }
187
188 36
    public function getSouthing(): Length
189
    {
190
        return $this->southing;
191
    }
192
193
    public function getHeight(): ?Length
194
    {
195 81
        return $this->height;
196
    }
197
198 81
    public function getCRS(): Projected
199 81
    {
200
        return $this->crs;
201
    }
202 81
203 27
    public function getCoordinateEpoch(): ?DateTimeImmutable
204
    {
205
        return $this->epoch;
206
    }
207 54
208 54
    /**
209 54
     * Calculate distance between two points.
210 54
     * Because this is a simple grid, we can use Pythagoras.
211
     */
212
    public function calculateDistance(Point $to): Length
213
    {
214
        try {
215
            if ($to instanceof ConvertiblePoint) {
216 162
                $to = $to->convert($this->crs);
217
            }
218 162
        } finally {
219 162
            if ($to->getCRS()->getSRID() !== $this->crs->getSRID()) {
220 162
                throw new InvalidCoordinateReferenceSystemException('Can only calculate distances between two points in the same CRS');
221 144
            }
222 162
223 153
            /* @var ProjectedPoint $to */
224 18
            return new Metre(
225 18
                sqrt(
226 9
                    ($to->getEasting()->getValue() - $this->getEasting()->getValue()) ** 2 +
227 9
                    ($to->getNorthing()->getValue() - $this->getNorthing()->getValue()) ** 2
228
                )
229
            );
230
        }
231
    }
232
233 162
    public function asGeographicPoint(): GeographicPoint
234
    {
235
        return $this->performOperation($this->crs->getDerivingConversion(), $this->crs->getBaseCRS(), true);
0 ignored issues
show
Bug Best Practice introduced by
The expression return $this->performOpe...rs->getBaseCRS(), true) returns the type PHPCoord\Point which includes types incompatible with the type-hinted return PHPCoord\GeographicPoint.
Loading history...
236
    }
237
238
    public function convert(CoordinateReferenceSystem $to, bool $ignoreBoundaryRestrictions = false): Point
239 72
    {
240
        if ($to->getSRID() === $this->crs->getBaseCRS()->getSRID()) {
241
            return $this->performOperation($this->crs->getDerivingConversion(), $this->crs->getBaseCRS(), true);
242
        }
243
244
        return $this->autoConvert($to, $ignoreBoundaryRestrictions);
245
    }
246
247
    public function __toString(): string
248
    {
249 72
        $values = [];
250 72
        foreach ($this->getCRS()->getCoordinateSystem()->getAxes() as $axis) {
251
            if ($axis->getName() === Axis::EASTING) {
252 72
                $values[] = $this->easting;
253 63
            } elseif ($axis->getName() === Axis::NORTHING) {
254 63
                $values[] = $this->northing;
255 63
            } elseif ($axis->getName() === Axis::WESTING) {
256 63
                $values[] = $this->westing;
257 63
            } elseif ($axis->getName() === Axis::SOUTHING) {
258 63
                $values[] = $this->southing;
259 63
            } elseif ($axis->getName() === Axis::ELLIPSOIDAL_HEIGHT) {
260
                $values[] = $this->height;
261 72
            } else {
262 72
                throw new UnknownAxisException(); // @codeCoverageIgnore
263 72
            }
264 72
        }
265 72
266 72
        return '(' . implode(', ', $values) . ')';
267
    }
268
269 72
    /**
270 72
     * Affine parametric transformation.
271
     */
272 72
    public function affineParametricTransform(
273
        Projected $to,
274
        Length $A0,
275
        Coefficient $A1,
276
        Coefficient $A2,
277
        Length $B0,
278 342
        Coefficient $B1,
279
        Coefficient $B2,
280
        bool $inReverse
281
    ): self {
282
        $xs = $this->easting->getValue(); // native unit to metre conversion already embedded in the scale factor
283
        $ys = $this->northing->getValue(); // native unit to metre conversion already embedded in the scale factor
284
285
        if ($inReverse) {
286
            $D = ($A1->getValue() * $B2->getValue()) - ($A2->getValue() * $B1->getValue());
287 342
            $a0 = (($A2->getValue() * $B0->asMetres()->getValue()) - ($B2->getValue() * $A0->asMetres()->getValue())) / $D;
288 342
            $b0 = (($B1->getValue() * $A0->asMetres()->getValue()) - ($A1->getValue() * $B0->asMetres()->getValue())) / $D;
289 342
            $a1 = $B2->getValue() / $D;
290 342
            $a2 = -$A2->getValue() / $D;
291 342
            $b1 = -$B1->getValue() / $D;
292 342
            $b2 = $A1->getValue() / $D;
293 342
        } else {
294 342
            $a0 = $A0->asMetres()->getValue();
295 342
            $a1 = $A1->getValue();
296 342
            $a2 = $A2->getValue();
297 342
            $b0 = $B0->asMetres()->getValue();
298
            $b1 = $B1->getValue();
299 342
            $b2 = $B2->getValue();
300 342
        }
301
302 342
        $xt = $a0 + ($a1 * $xs) + ($a2 * $ys);
303 342
        $yt = $b0 + ($b1 * $xs) + ($b2 * $ys);
304 342
305
        return static::create($to, new Metre($xt), new Metre($yt), new Metre(-$xt), new Metre(-$yt), $this->epoch);
306 342
    }
307 342
308 342
    /**
309 342
     * Albers Equal Area.
310 342
     */
311 342
    public function albersEqualArea(
312 342
        Geographic2D|Geographic3D $to,
313 288
        Angle $latitudeOfFalseOrigin,
314
        Angle $longitudeOfFalseOrigin,
315 54
        Angle $latitudeOf1stStandardParallel,
316
        Angle $latitudeOf2ndStandardParallel,
317
        Length $eastingAtFalseOrigin,
318 342
        Length $northingAtFalseOrigin
319 342
    ): GeographicPoint {
320
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
321 342
        $easting = $this->easting->asMetres()->getValue() - $eastingAtFalseOrigin->asMetres()->getValue();
322
        $northing = $this->northing->asMetres()->getValue() - $northingAtFalseOrigin->asMetres()->getValue();
323
        $phiOrigin = $latitudeOfFalseOrigin->asRadians()->getValue();
324
        $phi1 = $latitudeOf1stStandardParallel->asRadians()->getValue();
325
        $phi2 = $latitudeOf2ndStandardParallel->asRadians()->getValue();
326
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
327 45
        $e = $ellipsoid->getEccentricity();
328
        $e2 = $ellipsoid->getEccentricitySquared();
329
        $e4 = $e ** 4;
330
        $e6 = $e ** 6;
331
332
        $centralMeridianFirstParallel = cos($phi1) / sqrt(1 - $e2 * sin($phi1) ** 2);
333
        $centralMeridianSecondParallel = cos($phi2) / sqrt(1 - $e2 * sin($phi2) ** 2);
334 45
335 45
        $alphaOrigin = (1 - $e2) * (sin($phiOrigin) / (1 - $e2 * sin($phiOrigin) ** 2) - (1 / 2 / $e) * log((1 - $e * sin($phiOrigin)) / (1 + $e * sin($phiOrigin))));
336 45
        $alphaFirstParallel = (1 - $e2) * (sin($phi1) / (1 - $e2 * sin($phi1) ** 2) - (1 / 2 / $e) * log((1 - $e * sin($phi1)) / (1 + $e * sin($phi1))));
337 45
        $alphaSecondParallel = (1 - $e2) * (sin($phi2) / (1 - $e2 * sin($phi2) ** 2) - (1 / 2 / $e) * log((1 - $e * sin($phi2)) / (1 + $e * sin($phi2))));
338 45
339 45
        $n = ($centralMeridianFirstParallel ** 2 - $centralMeridianSecondParallel ** 2) / ($alphaSecondParallel - $alphaFirstParallel);
340 45
        $C = $centralMeridianFirstParallel ** 2 + $n * $alphaFirstParallel;
341 45
        $rhoOrigin = $a * sqrt($C - $n * $alphaOrigin) / $n;
342 45
        $rhoPrime = hypot($easting, $rhoOrigin - $northing);
343 45
        $alphaPrime = ($C - $rhoPrime ** 2 * $n ** 2 / $a ** 2) / $n;
344
        $betaPrime = self::asin($alphaPrime / (1 - (1 - $e2) / 2 / $e * log((1 - $e) / (1 + $e))));
345 45
        if ($n > 0) {
346 45
            $theta = atan2($easting, $rhoOrigin - $northing);
347 45
        } else {
348 45
            $theta = atan2(-$easting, $northing - $rhoOrigin);
349
        }
350 45
351
        $latitude = $betaPrime + (($e2 / 3 + 31 * $e4 / 180 + 517 * $e6 / 5040) * sin(2 * $betaPrime)) + ((23 * $e4 / 360 + 251 * $e6 / 3780) * sin(4 * $betaPrime)) + ((761 * $e6 / 45360) * sin(6 * $betaPrime));
352 45
        $longitude = $longitudeOfFalseOrigin->asRadians()->getValue() + ($theta / $n);
353
354
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
355
    }
356 45
357 45
    /**
358
     * American Polyconic.
359 45
     */
360 45
    public function americanPolyconic(
361
        Geographic2D|Geographic3D $to,
362 45
        Angle $latitudeOfNaturalOrigin,
363 45
        Angle $longitudeOfNaturalOrigin,
364 45
        Length $falseEasting,
365 45
        Length $falseNorthing
366 45
    ): GeographicPoint {
367 45
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
368
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
369 45
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
370
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
371
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
372 45
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
373
        $e = $ellipsoid->getEccentricity();
374
        $e2 = $ellipsoid->getEccentricitySquared();
375
        $e4 = $e ** 4;
376
        $e6 = $e ** 6;
377
378 9
        $i = (1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256);
379
        $ii = (3 * $e2 / 8 + 3 * $e4 / 32 + 45 * $e6 / 1024);
380
        $iii = (15 * $e4 / 256 + 45 * $e6 / 1024);
381
        $iv = (35 * $e6 / 3072);
382
383
        $MO = $a * ($i * $latitudeOrigin - $ii * sin(2 * $latitudeOrigin) + $iii * sin(4 * $latitudeOrigin) - $iv * sin(6 * $latitudeOrigin));
384
385 9
        if ($MO === $northing) {
386 9
            $latitude = 0;
387 9
            $longitude = $longitudeOrigin + $easting / $a;
388 9
        } else {
389 9
            $A = ($MO + $northing) / $a;
390 9
            $B = $A ** 2 + $easting ** 2 / $a ** 2;
391 9
392 9
            $latitude = $A;
393 9
            $C = sqrt(1 - $e2 * sin($latitude) ** 2) * tan($latitude);
394 9
            do {
395
                $latitudeN = $latitude;
396 9
                $Ma = $i * $latitude - $ii * sin(2 * $latitude) + $iii * sin(4 * $latitude) - $iv * sin(6 * $latitude);
397 9
                $MnPrime = $i - 2 * $ii * cos(2 * $latitude) + 4 * $iii * cos(4 * $latitude) - 6 * $iv * cos(6 * $latitude);
398 9
                $latitude = $latitude - ($A * ($C * $Ma + 1) - $Ma - $C * ($Ma ** 2 + $B) / 2) / ($e2 * sin(2 * $latitude) * ($Ma ** 2 + $B - 2 * $A * $Ma) / 4 * $C + ($A - $Ma) * ($C * $MnPrime - (2 / sin(2 * $latitude))) - $MnPrime);
399
                $C = sqrt(1 - $e2 * sin($latitude) ** 2) * tan($latitude);
400 9
            } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE_FORMULA);
401 9
402 9
            $longitude = $longitudeOrigin + (self::asin($easting * $C / $a)) / sin($latitude);
403
        }
404 9
405 9
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
406
    }
407 9
408
    /**
409 9
     * Bonne.
410 9
     */
411
    public function bonne(
412
        Geographic2D|Geographic3D $to,
413
        Angle $latitudeOfNaturalOrigin,
414
        Angle $longitudeOfNaturalOrigin,
415 9
        Length $falseEasting,
416
        Length $falseNorthing
417
    ): GeographicPoint {
418
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
419
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
420
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
421 27
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
422
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
423
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
424
        $e = $ellipsoid->getEccentricity();
425
        $e2 = $ellipsoid->getEccentricitySquared();
426
        $e4 = $e ** 4;
427
        $e6 = $e ** 6;
428 27
429 27
        $mO = cos($latitudeOrigin) / sqrt(1 - $e2 * sin($latitudeOrigin) ** 2);
430 27
        $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));
431 27
        $rho = hypot($easting, $a * $mO / sin($latitudeOrigin) - $northing) * static::sign($latitudeOrigin);
432 27
433 27
        $M = $a * $mO / sin($latitudeOrigin) + $MO - $rho;
434 27
        $mu = $M / ($a * (1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256));
435 27
        $e1 = (1 - sqrt(1 - $e2)) / (1 + sqrt(1 - $e2));
436 27
437 27
        $latitude = $mu + ((3 * $e1 / 2) - (27 * $e1 ** 3 / 32)) * sin(2 * $mu) + ((21 * $e1 ** 2 / 16) - (55 * $e1 ** 4 / 32)) * sin(4 * $mu) + ((151 * $e1 ** 3 / 96)) * sin(6 * $mu) + ((1097 * $e1 ** 4 / 512)) * sin(8 * $mu);
438
        $m = cos($latitude) / sqrt(1 - $e2 * sin($latitude) ** 2);
439 27
440 27
        if ($m === 0.0) {
441 27
            $longitude = $longitudeOrigin; // pole
442
        } elseif ($latitudeOrigin >= 0) {
443 27
            $longitude = $longitudeOrigin + $rho * atan2($easting, $a * $mO / sin($latitudeOrigin) - $northing) / $a / $m;
444 27
        } else {
445 27
            $longitude = $longitudeOrigin + $rho * atan2(-$easting, -($a * $mO / sin($latitudeOrigin) - $northing)) / $a / $m;
446
        }
447 27
448 27
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
449
    }
450 27
451
    /**
452 27
     * Bonne South Orientated.
453 27
     */
454
    public function bonneSouthOrientated(
455
        Geographic2D|Geographic3D $to,
456
        Angle $latitudeOfNaturalOrigin,
457
        Angle $longitudeOfNaturalOrigin,
458 27
        Length $falseEasting,
459
        Length $falseNorthing
460
    ): GeographicPoint {
461
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
462
        $westing = $falseEasting->asMetres()->getValue() - $this->westing->asMetres()->getValue();
463
        $southing = $falseNorthing->asMetres()->getValue() - $this->southing->asMetres()->getValue();
464
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
465
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
466 360
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
467
        $e = $ellipsoid->getEccentricity();
468
        $e2 = $ellipsoid->getEccentricitySquared();
469
        $e4 = $e ** 4;
470
        $e6 = $e ** 6;
471 360
472 360
        $mO = cos($latitudeOrigin) / sqrt(1 - $e2 * sin($latitudeOrigin) ** 2);
473
        $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));
474 360
        $rho = hypot($westing, $a * $mO / sin($latitudeOrigin) - $southing) * static::sign($latitudeOrigin);
475
476
        $M = $a * $mO / sin($latitudeOrigin) + $MO - $rho;
477
        $mu = $M / ($a * (1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256));
478
        $e1 = (1 - sqrt(1 - $e2)) / (1 + sqrt(1 - $e2));
479
480 279
        $latitude = $mu + ((3 * $e1 / 2) - (27 * $e1 ** 3 / 32)) * sin(2 * $mu) + ((21 * $e1 ** 2 / 16) - (55 * $e1 ** 4 / 32)) * sin(4 * $mu) + ((151 * $e1 ** 3 / 96)) * sin(6 * $mu) + ((1097 * $e1 ** 4 / 512)) * sin(8 * $mu);
481
        $m = cos($latitude) / sqrt(1 - $e2 * sin($latitude) ** 2);
482
483
        if ($m === 0.0) {
484
            $longitude = $longitudeOrigin; // pole
485
        } elseif ($latitudeOrigin >= 0) {
486
            $longitude = $longitudeOrigin + $rho * atan2($westing, $a * $mO / sin($latitudeOrigin) - $southing) / $a / $m;
487 279
        } else {
488 279
            $longitude = $longitudeOrigin + $rho * atan2(-$westing, -($a * $mO / sin($latitudeOrigin) - $southing)) / $a / $m;
489 279
        }
490 279
491 279
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
492 279
    }
493 279
494 279
    /**
495 279
     * Cartesian Grid Offsets
496 279
     * This transformation allows calculation of coordinates in the target system by adding the parameter value to the
497
     * coordinate values of the point in the source system.
498 279
     */
499
    public function offsets(
500 279
        Projected $to,
501 279
        Length $eastingOffset,
502 279
        Length $northingOffset
503 279
    ): self {
504
        $easting = $this->easting->asMetres()->getValue() + $eastingOffset->asMetres()->getValue();
505 279
        $northing = $this->northing->asMetres()->getValue() + $northingOffset->asMetres()->getValue();
506 279
507
        return static::create($to, new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $this->epoch);
508 279
    }
509 279
510
    /**
511 279
     * Cassini-Soldner.
512 279
     */
513
    public function cassiniSoldner(
514 279
        Geographic2D|Geographic3D $to,
515
        Angle $latitudeOfNaturalOrigin,
516
        Angle $longitudeOfNaturalOrigin,
517
        Length $falseEasting,
518
        Length $falseNorthing
519
    ): GeographicPoint {
520 18
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
521
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
522
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
523
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
524
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
525
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
526
        $e = $ellipsoid->getEccentricity();
527 18
        $e2 = $ellipsoid->getEccentricitySquared();
528 18
        $e4 = $e ** 4;
529 18
        $e6 = $e ** 6;
530 18
531 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));
532 18
533 18
        $e1 = (1 - sqrt(1 - $e2)) / (1 + sqrt(1 - $e2));
534 18
        $M = $MO + $northing;
535 18
        $mu = $M / ($a * (1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256));
536 18
        $latitudeCentralMeridian = $mu + (3 * $e1 / 2 - 27 * $e1 ** 3 / 32) * sin(2 * $mu) + (21 * $e1 ** 2 / 16 - 55 * $e1 ** 4 / 32) * sin(4 * $mu) + (151 * $e1 ** 3 / 96) * sin(6 * $mu) + (1097 * $e1 ** 4 / 512) * sin(8 * $mu);
537
538 18
        $nu = $a / sqrt((1 - $e2 * sin($latitudeCentralMeridian) ** 2));
539 18
        $rho = $a * (1 - $e2) / (1 - $e2 * sin($latitudeCentralMeridian) ** 2) ** 1.5;
540
541 18
        $T = tan($latitudeCentralMeridian) ** 2;
542
        $D = $easting / $nu;
543 18
544 18
        $latitude = $latitudeCentralMeridian - ($nu * tan($latitudeCentralMeridian) / $rho) * ($D ** 2 / 2 - (1 + 3 * $T) * $D ** 4 / 24);
545
        $longitude = $longitudeOrigin + ($D - $T * $D ** 3 / 3 + (1 + 3 * $T) * $T * $D ** 5 / 15) / cos($latitudeCentralMeridian);
546 18
547 18
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
548 18
    }
549
550 18
    /**
551 18
     * Hyperbolic Cassini-Soldner.
552
     */
553 18
    public function hyperbolicCassiniSoldner(
554 18
        Geographic2D|Geographic3D $to,
555
        Angle $latitudeOfNaturalOrigin,
556 18
        Angle $longitudeOfNaturalOrigin,
557 18
        Length $falseEasting,
558
        Length $falseNorthing
559 18
    ): GeographicPoint {
560
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
561
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
562
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
563
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
564
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
565 297
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
566
        $e = $ellipsoid->getEccentricity();
567
        $e2 = $ellipsoid->getEccentricitySquared();
568
        $e4 = $e ** 4;
569
        $e6 = $e ** 6;
570
571
        $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));
572
        $e1 = (1 - sqrt(1 - $e2)) / (1 + sqrt(1 - $e2));
573 297
574 297
        $latitude1 = $latitudeOrigin + $northing / 1567446;
575 297
576 297
        $nu = $a / sqrt((1 - $e2 * sin($latitude1) ** 2));
577 297
        $rho = $a * (1 - $e2) / (1 - $e2 * sin($latitude1) ** 2) ** 1.5;
578 297
579 297
        $qPrime = $northing ** 3 / (6 * $rho * $nu);
580 297
        $q = ($northing + $qPrime) ** 3 / (6 * $rho * $nu);
581
        $M = $MO + $northing + $q;
582 297
583
        $mu = $M / ($a * (1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256));
584 297
        $latitudeCentralMeridian = $mu + (3 * $e1 / 2 - 27 * $e1 ** 3 / 32) * sin(2 * $mu) + (21 * $e1 ** 2 / 16 - 55 * $e1 ** 4 / 32) * sin(4 * $mu) + (151 * $e1 ** 3 / 96) * sin(6 * $mu) + (1097 * $e1 ** 4 / 512) * sin(8 * $mu);
585
586 297
        $T = tan($latitudeCentralMeridian) ** 2;
587 297
        $D = $easting / $nu;
588 297
589
        $latitude = $latitudeCentralMeridian - ($nu * tan($latitudeCentralMeridian) / $rho) * ($D ** 2 / 2 - (1 + 3 * $T) * $D ** 4 / 24);
590 297
        $longitude = $longitudeOrigin + ($D - $T * $D ** 3 / 3 + (1 + 3 * $T) * $T * $D ** 5 / 15) / cos($latitudeCentralMeridian);
591 297
592 297
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
593
    }
594 297
595
    /**
596
     * Colombia Urban.
597
     */
598
    public function columbiaUrban(
599
        Geographic2D|Geographic3D $to,
600 36
        Angle $latitudeOfNaturalOrigin,
601
        Angle $longitudeOfNaturalOrigin,
602
        Length $falseEasting,
603
        Length $falseNorthing,
604
        Length $projectionPlaneOriginHeight
605
    ): GeographicPoint {
606 36
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
607 36
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
608 36
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
609 36
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
610 36
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
611 36
        $heightOrigin = $projectionPlaneOriginHeight->asMetres()->getValue();
612 36
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
613 36
        $e2 = $ellipsoid->getEccentricitySquared();
614 36
615
        $rhoOrigin = $a * (1 - $e2) / (1 - $e2 * sin($latitudeOrigin) ** 2) ** 1.5;
616 36
617 36
        $nuOrigin = $a / sqrt(1 - $e2 * (sin($latitudeOrigin) ** 2));
618
619 36
        $B = tan($latitudeOrigin) / (2 * $rhoOrigin * $nuOrigin);
620
        $C = 1 + $heightOrigin / $a;
621 36
        $D = $rhoOrigin * (1 + $heightOrigin / ($a * (1 - $e2)));
622 36
623 36
        $latitude = $latitudeOrigin + ($northing / $D) - $B * ($easting / $C) ** 2;
624 36
        $nu = $a / sqrt(1 - $e2 * (sin($latitude) ** 2));
625
        $longitude = $longitudeOrigin + $easting / ($C * $nu * cos($latitude));
626 36
627
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
628 36
    }
629 36
630
    /**
631 36
     * Equal Earth.
632
     */
633
    public function equalEarth(
634
        Geographic2D|Geographic3D $to,
635
        Angle $longitudeOfNaturalOrigin,
636
        Length $falseEasting,
637
        Length $falseNorthing
638 18
    ): GeographicPoint {
639
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
640
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
641
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
642
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
643
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
644
        $e = $ellipsoid->getEccentricity();
645 18
        $e2 = $ellipsoid->getEccentricitySquared();
646 18
        $e4 = $e ** 4;
647 18
        $e6 = $e ** 6;
648 18
649 18
        $qP = (1 - $e2) * ((1 / (1 - $e2)) - (1 / (2 * $e) * log((1 - $e) / (1 + $e))));
650 18
        $Rq = $a * sqrt($qP / 2);
651 18
652 18
        $theta = $northing / $Rq;
653 18
        do {
654 18
            $thetaN = $theta;
655 18
            $correctionFactor = ($theta * (1.340264 - 0.081106 * $theta ** 2 + $theta ** 6 * (0.000893 + 0.003796 * $theta ** 2)) - $northing / $Rq) / (1.340264 - 0.243318 * $theta ** 2 + $theta ** 6 * (0.006251 + 0.034164 * $theta ** 2));
656 18
            $theta = $theta - $correctionFactor;
657 18
        } while (abs($theta - $thetaN) >= static::ITERATION_CONVERGENCE_FORMULA);
658 18
659
        $beta = self::asin(2 * sin($theta) / sqrt(3));
660 18
661 18
        $latitude = $beta + (($e2 / 3 + 31 * $e4 / 180 + 517 * $e6 / 5040) * sin(2 * $beta)) + ((23 * $e4 / 360 + 251 * $e6 / 3780) * sin(4 * $beta)) + ((761 * $e6 / 45360) * sin(6 * $beta));
662 18
        $longitude = $longitudeOrigin + sqrt(3) * $easting * (1.340264 - 0.243318 * $theta ** 2 + $theta ** 6 * (0.006251 + 0.034164 * $theta ** 2)) / (2 * $Rq * cos($theta));
663 18
664 18
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
665 18
    }
666 18
667 18
    /**
668
     * Equidistant Cylindrical
669 18
     * See method code 1029 for spherical development. See also Pseudo Plate Carree, method code 9825.
670 18
     */
671 18
    public function equidistantCylindrical(
672 18
        Geographic2D|Geographic3D $to,
673 18
        Angle $latitudeOf1stStandardParallel,
674 18
        Angle $longitudeOfNaturalOrigin,
675 18
        Length $falseEasting,
676
        Length $falseNorthing
677 18
    ): GeographicPoint {
678
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
679 18
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
680
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
681
        $latitudeFirstParallel = $latitudeOf1stStandardParallel->asRadians()->getValue();
682
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
683
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
684
        $e = $ellipsoid->getEccentricity();
685
        $e2 = $ellipsoid->getEccentricitySquared();
686 18
        $e4 = $e ** 4;
687
        $e6 = $e ** 6;
688
        $e8 = $e ** 8;
689
        $e10 = $e ** 10;
690
        $e12 = $e ** 12;
691
        $e14 = $e ** 14;
692
693 18
        $n = (1 - sqrt(1 - $e2)) / (1 + sqrt(1 - $e2));
694 18
        $n2 = $n ** 2;
695 18
        $n3 = $n ** 3;
696 18
        $n4 = $n ** 4;
697 18
        $n5 = $n ** 5;
698 18
        $n6 = $n ** 6;
699 18
        $n7 = $n ** 7;
700 18
        $mu = $northing / ($a * (1 - 1 / 4 * $e2 - 3 / 64 * $e4 - 5 / 256 * $e6 - 175 / 16384 * $e8 - 441 / 65536 * $e10 - 4851 / 1048576 * $e12 - 14157 / 4194304 * $e14));
701 18
702 18
        $latitude = $mu + (3 / 2 * $n - 27 / 32 * $n3 + 269 / 512 * $n5 - 6607 / 24576 * $n7) * sin(2 * $mu)
703
            + (21 / 16 * $n2 - 55 / 32 * $n4 + 6759 / 4096 * $n6) * sin(4 * $mu)
704 18
            + (151 / 96 * $n3 - 417 / 128 * $n5 + 87963 / 20480 * $n7) * sin(6 * $mu)
705 18
            + (1097 / 512 * $n4 - 15543 / 2560 * $n6) * sin(8 * $mu)
706 18
            + (8011 / 2560 * $n5 - 69119 / 6144 * $n7) * sin(10 * $mu)
707
            + (293393 / 61440 * $n6) * sin(12 * $mu)
708 18
            + (6845701 / 860160 * $n7) * sin(14 * $mu);
709
710 18
        $longitude = $longitudeOrigin + $easting * sqrt(1 - $e2 * sin($latitudeFirstParallel) ** 2) / ($a * cos($latitudeFirstParallel));
711 18
712 18
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
713 18
    }
714 18
715
    /**
716 18
     * Guam Projection
717
     * Simplified form of Oblique Azimuthal Equidistant projection method.
718 18
     */
719
    public function guamProjection(
720
        Geographic2D|Geographic3D $to,
721
        Angle $latitudeOfNaturalOrigin,
722
        Angle $longitudeOfNaturalOrigin,
723
        Length $falseEasting,
724 126
        Length $falseNorthing
725
    ): GeographicPoint {
726
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
727
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
728
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
729
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
730
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
731
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
732
        $e = $ellipsoid->getEccentricity();
733
        $e2 = $ellipsoid->getEccentricitySquared();
734 126
        $e4 = $e ** 4;
735 126
        $e6 = $e ** 6;
736 126
737 126
        $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));
738 126
        $e1 = (1 - sqrt(1 - $e2)) / (1 + sqrt(1 - $e2));
739 126
        $i = (1 - $e2 / 4 - 3 * $e4 / 64 - 5 * $e6 / 256);
740 126
741 126
        $latitude = $latitudeOrigin;
742 126
        do {
743 126
            $latitudeN = $latitude;
744 126
            $M = $MO + $northing - ($easting ** 2 * tan($latitude) * sqrt(1 - $e2 * sin($latitude) ** 2) / (2 * $a));
745 126
            $mu = $M / ($a * $i);
746
            $latitude = $mu + (3 * $e1 / 2 - 27 * $e1 ** 3 / 32) * sin(2 * $mu) + (21 * $e1 ** 2 / 16 - 55 * $e1 ** 4 / 32) * sin(4 * $mu) + (151 * $e1 ** 3 / 96) * sin(6 * $mu) + (1097 * $e1 ** 4 / 512) * sin(8 * $mu);
747 126
        } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE_FORMULA);
748 126
749 126
        $longitude = $longitudeOrigin + $easting * sqrt(1 - $e2 * sin($latitude) ** 2) / ($a * cos($latitude));
750 126
751 126
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
752 126
    }
753
754 126
    /**
755 126
     * Krovak.
756 126
     */
757 126
    public function krovak(
758 126
        Geographic2D|Geographic3D $to,
759 126
        Angle $latitudeOfProjectionCentre,
760
        Angle $longitudeOfOrigin,
761 126
        Angle $coLatitudeOfConeAxis,
762
        Angle $latitudeOfPseudoStandardParallel,
763 126
        Scale $scaleFactorOnPseudoStandardParallel,
764 126
        Length $falseEasting,
765 126
        Length $falseNorthing
766
    ): GeographicPoint {
767 126
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
768
        $longitudeOffset = $this->crs->getDatum()->getPrimeMeridian()->getGreenwichLongitude()->asRadians()->getValue() - $to->getDatum()->getPrimeMeridian()->getGreenwichLongitude()->asRadians()->getValue();
769 126
        $westing = $this->westing->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
770
        $southing = $this->southing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
771
        $latitudeC = $latitudeOfProjectionCentre->asRadians()->getValue();
772
        $longitudeO = $longitudeOfOrigin->asRadians()->getValue();
773
        $alphaC = $coLatitudeOfConeAxis->asRadians()->getValue();
774
        $latitudeP = $latitudeOfPseudoStandardParallel->asRadians()->getValue();
775
        $kP = $scaleFactorOnPseudoStandardParallel->asUnity()->getValue();
776
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
777 54
        $e = $ellipsoid->getEccentricity();
778
        $e2 = $ellipsoid->getEccentricitySquared();
779
780
        $A = $a * sqrt(1 - $e2) / (1 - $e2 * sin($latitudeC) ** 2);
781
        $B = sqrt(1 + $e2 * cos($latitudeC) ** 4 / (1 - $e2));
782
        $upsilonO = self::asin(sin($latitudeC) / $B);
783
        $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);
784
        $n = sin($latitudeP);
785
        $rO = $kP * $A / tan($latitudeP);
786
787
        $r = hypot($southing, $westing) ?: 1;
788
        $theta = atan2($westing, $southing);
789
        $D = $theta / sin($latitudeP);
790
        $T = 2 * (atan(($rO / $r) ** (1 / $n) * tan(M_PI / 4 + $latitudeP / 2)) - M_PI / 4);
791
        $U = self::asin(cos($alphaC) * sin($T) - sin($alphaC) * cos($T) * cos($D));
792
        $V = self::asin(cos($T) * sin($D) / cos($U));
793
794
        $latitude = $U;
795
        do {
796
            $latitudeN = $latitude;
797
            $latitude = 2 * (atan($tO ** (-1 / $B) * tan($U / 2 + M_PI / 4) ** (1 / $B) * ((1 + $e * sin($latitude)) / (1 - $e * sin($latitude))) ** ($e / 2)) - M_PI / 4);
798
        } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE_FORMULA);
799 54
800 54
        $longitude = $longitudeO + $longitudeOffset - $V / $B;
801 54
802 54
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
803 54
    }
804 54
805 54
    /**
806 54
     * Krovak Modified
807 54
     * Incorporates a polynomial transformation which is defined to be exact and for practical purposes is considered
808 54
     * to be a map projection.
809 54
     */
810 54
    public function krovakModified(
811
        Geographic2D|Geographic3D $to,
812 54
        Angle $latitudeOfProjectionCentre,
813 54
        Angle $longitudeOfOrigin,
814
        Angle $coLatitudeOfConeAxis,
815 54
        Angle $latitudeOfPseudoStandardParallel,
816 54
        Scale $scaleFactorOnPseudoStandardParallel,
817
        Length $falseEasting,
818 54
        Length $falseNorthing,
819
        Length $ordinate1OfEvaluationPoint,
820 54
        Length $ordinate2OfEvaluationPoint,
821
        Coefficient $C1,
822
        Coefficient $C2,
823
        Coefficient $C3,
824
        Coefficient $C4,
825
        Coefficient $C5,
826
        Coefficient $C6,
827 135
        Coefficient $C7,
828
        Coefficient $C8,
829
        Coefficient $C9,
830
        Coefficient $C10
831
    ): GeographicPoint {
832
        $Xr = $this->getSouthing()->asMetres()->getValue() - $falseNorthing->asMetres()->getValue() - $ordinate1OfEvaluationPoint->asMetres()->getValue();
833
        $Yr = $this->getWesting()->asMetres()->getValue() - $falseEasting->asMetres()->getValue() - $ordinate2OfEvaluationPoint->asMetres()->getValue();
834 135
        $c1 = $C1->asUnity()->getValue();
835 135
        $c2 = $C2->asUnity()->getValue();
836 135
        $c3 = $C3->asUnity()->getValue();
837 135
        $c4 = $C4->asUnity()->getValue();
838 135
        $c5 = $C5->asUnity()->getValue();
839 135
        $c6 = $C6->asUnity()->getValue();
840 135
        $c7 = $C7->asUnity()->getValue();
841 135
        $c8 = $C8->asUnity()->getValue();
842 135
        $c9 = $C9->asUnity()->getValue();
843 135
        $c10 = $C10->asUnity()->getValue();
844
845 135
        $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);
846 135
        $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);
847 135
848 135
        $Xp = $this->getSouthing()->asMetres()->getValue() - $falseNorthing->asMetres()->getValue() + $dX;
849 135
        $Yp = $this->getWesting()->asMetres()->getValue() - $falseEasting->asMetres()->getValue() + $dY;
850 135
851 135
        $asKrovak = self::create($this->crs, new Metre(-$Yp), new Metre(-$Xp), new Metre($Yp), new Metre($Xp), $this->epoch);
852 135
853
        return $asKrovak->krovak($to, $latitudeOfProjectionCentre, $longitudeOfOrigin, $coLatitudeOfConeAxis, $latitudeOfPseudoStandardParallel, $scaleFactorOnPseudoStandardParallel, new Metre(0), new Metre(0));
854 135
    }
855 135
856
    /**
857 135
     * Lambert Azimuthal Equal Area
858
     * This is the ellipsoidal form of the projection.
859
     */
860
    public function lambertAzimuthalEqualArea(
861
        Geographic2D|Geographic3D $to,
862
        Angle $latitudeOfNaturalOrigin,
863
        Angle $longitudeOfNaturalOrigin,
864
        Length $falseEasting,
865
        Length $falseNorthing
866 18
    ): GeographicPoint {
867
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
868
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
869
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
870
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
871
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
872
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
873 18
        $e = $ellipsoid->getEccentricity();
874 18
        $e2 = $ellipsoid->getEccentricitySquared();
875 18
        $e4 = $e ** 4;
876 18
        $e6 = $e ** 6;
877 18
878 18
        $qO = (1 - $e2) * ((sin($latitudeOrigin) / (1 - $e2 * sin($latitudeOrigin) ** 2)) - ((1 / (2 * $e)) * log((1 - $e * sin($latitudeOrigin)) / (1 + $e * sin($latitudeOrigin)))));
879
        $qP = (1 - $e2) * ((1 / (1 - $e2)) - ((1 / (2 * $e)) * log((1 - $e) / (1 + $e))));
880 18
        $betaO = self::asin($qO / $qP);
881 18
        $Rq = $a * sqrt($qP / 2);
882
        $D = $a * (cos($latitudeOrigin) / sqrt(1 - $e2 * sin($latitudeOrigin) ** 2)) / ($Rq * cos($betaO));
883 18
        $rho = hypot($easting / $D, $D * $northing) ?: 1;
884
        $C = 2 * self::asin($rho / (2 * $Rq));
885 18
        $beta = self::asin(cos($C) * sin($betaO) + ($D * $northing * sin($C) * cos($betaO)) / $rho);
886
887 18
        $latitude = $beta + (($e2 / 3 + 31 * $e4 / 180 + 517 * $e6 / 5040) * sin(2 * $beta)) + ((23 * $e4 / 360 + 251 * $e6 / 3780) * sin(4 * $beta)) + ((761 * $e6 / 45360) * sin(6 * $beta));
888
        $longitude = $longitudeOrigin + atan2($easting * sin($C), $D * $rho * cos($betaO) * cos($C) - $D ** 2 * $northing * sin($betaO) * sin($C));
889
890 18
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
891 18
    }
892 18
893 9
    /**
894
     * Lambert Azimuthal Equal Area (Spherical)
895
     * This is the spherical form of the projection.  See coordinate operation method Lambert Azimuthal Equal Area
896
     * (code 9820) for ellipsoidal form.  Differences of several tens of metres result from comparison of the two
897 18
     * methods.
898
     */
899
    public function lambertAzimuthalEqualAreaSpherical(
900
        Geographic2D|Geographic3D $to,
901
        Angle $latitudeOfNaturalOrigin,
902
        Angle $longitudeOfNaturalOrigin,
903 1953
        Length $falseEasting,
904
        Length $falseNorthing
905
    ): GeographicPoint {
906
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
907
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
908
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
909
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
910
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
911 1953
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
912 1953
913 1953
        $rho = hypot($easting, $northing) ?: 1;
914 1953
        $c = 2 * self::asin($rho / (2 * $a));
915 1953
916 1953
        $latitude = self::asin(cos($c) * sin($latitudeOrigin) + ($northing * sin($c) * cos($latitudeOrigin) / $rho));
917 1953
918 1953
        if ($latitudeOrigin === 90) {
0 ignored issues
show
introduced by
The condition $latitudeOrigin === 90 is always false.
Loading history...
919 1953
            $longitude = $longitudeOrigin + atan($easting / -$northing);
920
        } elseif ($latitudeOrigin === -90) {
0 ignored issues
show
introduced by
The condition $latitudeOrigin === -90 is always false.
Loading history...
921 1953
            $longitude = $longitudeOrigin + atan($easting / $northing);
922 1953
        } else {
923 1953
            $longitudeDenominator = ($rho * cos($latitudeOrigin) * cos($c) - $northing * sin($latitudeOrigin) * sin($c));
924 1953
            $longitude = $longitudeOrigin + atan($easting * sin($c) / $longitudeDenominator);
925 1953
            if ($longitudeDenominator < 0) {
926 1953
                $longitude += M_PI;
927 1953
            }
928 1935
        }
929
930 18
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
931 18
    }
932
933
    /**
934 1953
     * Lambert Conic Conformal (1SP).
935
     */
936 1953
    public function lambertConicConformal1SP(
937
        Geographic2D|Geographic3D $to,
938 1953
        Angle $latitudeOfNaturalOrigin,
939 1953
        Angle $longitudeOfNaturalOrigin,
940 1953
        Scale $scaleFactorAtNaturalOrigin,
941
        Length $falseEasting,
942 1953
        Length $falseNorthing
943
    ): GeographicPoint {
944 1953
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
945
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
946
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
947
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
948
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
949
        $scaleFactorOrigin = $scaleFactorAtNaturalOrigin->asUnity()->getValue();
950 135
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
951
        $e = $ellipsoid->getEccentricity();
952
        $e2 = $ellipsoid->getEccentricitySquared();
953
954
        $mO = cos($latitudeOrigin) / sqrt(1 - $e2 * sin($latitudeOrigin) ** 2);
955
        $tO = tan(M_PI / 4 - $latitudeOrigin / 2) / ((1 - $e * sin($latitudeOrigin)) / (1 + $e * sin($latitudeOrigin))) ** ($e / 2);
956
        $n = sin($latitudeOrigin);
957
        $F = $mO / ($n * $tO ** $n);
958 135
        $rO = $a * $F * $tO ** $n * $scaleFactorOrigin;
959 135
        $r = hypot($easting, $rO - $northing);
960 135
        if ($n >= 0) {
961 135
            $theta = atan2($easting, $rO - $northing);
962 135
        } else {
963 135
            $r = -$r;
964 135
            $theta = atan2(-$easting, -($rO - $northing));
965 135
        }
966 135
967
        $t = ($r / ($a * $scaleFactorOrigin * $F)) ** (1 / $n);
968 135
969 135
        $latitude = M_PI / (2 - 2 * atan($t));
970 135
        do {
971 135
            $latitudeN = $latitude;
972 135
            $latitude = M_PI / 2 - 2 * atan($t * ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2));
973 135
        } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE_FORMULA);
974 135
975 135
        $longitude = $theta / $n + $longitudeOrigin;
976
977
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
978
    }
979
980
    /**
981 135
     * Lambert Conic Conformal (west orientated).
982
     */
983 135
    public function lambertConicConformalWestOrientated(
984
        Geographic2D|Geographic3D $to,
985 135
        Angle $latitudeOfNaturalOrigin,
986 135
        Angle $longitudeOfNaturalOrigin,
987 135
        Scale $scaleFactorAtNaturalOrigin,
988
        Length $falseEasting,
989 135
        Length $falseNorthing
990
    ): GeographicPoint {
991 135
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
992
        $westing = $falseEasting->asMetres()->getValue() - $this->westing->asMetres()->getValue();
993
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
994
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
995
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
996
        $scaleFactorOrigin = $scaleFactorAtNaturalOrigin->asUnity()->getValue();
997 9
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
998
        $e = $ellipsoid->getEccentricity();
999
        $e2 = $ellipsoid->getEccentricitySquared();
1000
1001
        $mO = cos($latitudeOrigin) / sqrt(1 - $e2 * sin($latitudeOrigin) ** 2);
1002
        $tO = tan(M_PI / 4 - $latitudeOrigin / 2) / ((1 - $e * sin($latitudeOrigin)) / (1 + $e * sin($latitudeOrigin))) ** ($e / 2);
1003
        $n = sin($latitudeOrigin);
1004
        $F = $mO / ($n * $tO ** $n);
1005
        $rO = $a * $F * $tO ** $n ** $scaleFactorOrigin;
1006 9
        $r = hypot($westing, $rO - $northing);
1007 9
        if ($n >= 0) {
1008 9
            $theta = atan2($westing, $rO - $northing);
1009 9
        } else {
1010 9
            $r = -$r;
1011 9
            $theta = atan2(-$westing, -($rO - $northing));
1012 9
        }
1013 9
1014 9
        $t = ($r / ($a * $scaleFactorOrigin * $F)) ** (1 / $n);
1015 9
1016
        $latitude = M_PI / (2 - 2 * atan($t));
1017 9
        do {
1018 9
            $latitudeN = $latitude;
1019 9
            $latitude = M_PI / 2 - 2 * atan($t * ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2));
1020 9
        } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE_FORMULA);
1021 9
1022 9
        $longitude = $theta / $n + $longitudeOrigin;
1023 9
1024 9
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
1025 9
    }
1026
1027
    /**
1028
     * Lambert Conic Conformal (1SP) Variant B.
1029
     */
1030
    public function lambertConicConformal1SPVariantB(
1031 9
        Geographic2D|Geographic3D $to,
1032
        Angle $latitudeOfNaturalOrigin,
1033 9
        Scale $scaleFactorAtNaturalOrigin,
1034
        Angle $latitudeOfFalseOrigin,
1035 9
        Angle $longitudeOfFalseOrigin,
1036 9
        Length $eastingAtFalseOrigin,
1037 9
        Length $northingAtFalseOrigin
1038
    ): GeographicPoint {
1039 9
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1040
        $easting = $this->easting->asMetres()->getValue() - $eastingAtFalseOrigin->asMetres()->getValue();
1041 9
        $northing = $this->northing->asMetres()->getValue() - $northingAtFalseOrigin->asMetres()->getValue();
1042
        $latitudeNaturalOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
1043
        $latitudeFalseOrigin = $latitudeOfFalseOrigin->asRadians()->getValue();
1044
        $longitudeFalseOrigin = $longitudeOfFalseOrigin->asRadians()->getValue();
1045
        $scaleFactorOrigin = $scaleFactorAtNaturalOrigin->asUnity()->getValue();
1046
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1047 8244
        $e = $ellipsoid->getEccentricity();
1048
        $e2 = $ellipsoid->getEccentricitySquared();
1049
1050
        $mO = cos($latitudeNaturalOrigin) / sqrt(1 - $e2 * sin($latitudeNaturalOrigin) ** 2);
1051
        $tO = tan(M_PI / 4 - $latitudeNaturalOrigin / 2) / ((1 - $e * sin($latitudeNaturalOrigin)) / (1 + $e * sin($latitudeNaturalOrigin))) ** ($e / 2);
1052
        $tF = tan(M_PI / 4 - $latitudeFalseOrigin / 2) / ((1 - $e * sin($latitudeFalseOrigin)) / (1 + $e * sin($latitudeFalseOrigin))) ** ($e / 2);
1053
        $n = sin($latitudeNaturalOrigin);
1054
        $F = $mO / ($n * $tO ** $n);
1055
        $rF = $a * $F * $tF ** $n * $scaleFactorOrigin;
1056 8244
        $r = hypot($easting, $rF - $northing);
1057 8244
        if ($n >= 0) {
1058 8244
            $theta = atan2($easting, $rF - $northing);
1059 8244
        } else {
1060 8244
            $r = -$r;
1061 8244
            $theta = atan2(-$easting, -($rF - $northing));
1062 8244
        }
1063 8244
1064 8244
        $t = ($r / ($a * $scaleFactorOrigin * $F)) ** (1 / $n);
1065 8244
1066
        $latitude = M_PI / (2 - 2 * atan($t));
1067 8244
        do {
1068 8244
            $latitudeN = $latitude;
1069 8244
            $latitude = M_PI / 2 - 2 * atan($t * ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2));
1070 8244
        } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE_FORMULA);
1071 8244
1072 8244
        $longitude = $theta / $n + $longitudeFalseOrigin;
1073 8244
1074 8244
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
1075 8244
    }
1076 8244
1077 8244
    /**
1078 7407
     * Lambert Conic Conformal (2SP).
1079
     */
1080 837
    public function lambertConicConformal2SP(
1081
        Geographic2D|Geographic3D $to,
1082
        Angle $latitudeOfFalseOrigin,
1083 8244
        Angle $longitudeOfFalseOrigin,
1084
        Angle $latitudeOf1stStandardParallel,
1085 8244
        Angle $latitudeOf2ndStandardParallel,
1086 8244
        Length $eastingAtFalseOrigin,
1087 8244
        Length $northingAtFalseOrigin
1088
    ): GeographicPoint {
1089 8244
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1090
        $easting = $this->easting->asMetres()->getValue() - $eastingAtFalseOrigin->asMetres()->getValue();
1091 8244
        $northing = $this->northing->asMetres()->getValue() - $northingAtFalseOrigin->asMetres()->getValue();
1092
        $lambdaOrigin = $longitudeOfFalseOrigin->asRadians()->getValue();
1093
        $phiF = $latitudeOfFalseOrigin->asRadians()->getValue();
1094
        $phi1 = $latitudeOf1stStandardParallel->asRadians()->getValue();
1095
        $phi2 = $latitudeOf2ndStandardParallel->asRadians()->getValue();
1096
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1097 36
        $e = $ellipsoid->getEccentricity();
1098
        $e2 = $ellipsoid->getEccentricitySquared();
1099
1100
        $m1 = cos($phi1) / sqrt(1 - $e2 * sin($phi1) ** 2);
1101
        $m2 = cos($phi2) / sqrt(1 - $e2 * sin($phi2) ** 2);
1102
        $t1 = tan(M_PI / 4 - $phi1 / 2) / ((1 - $e * sin($phi1)) / (1 + $e * sin($phi1))) ** ($e / 2);
1103
        $t2 = tan(M_PI / 4 - $phi2 / 2) / ((1 - $e * sin($phi2)) / (1 + $e * sin($phi2))) ** ($e / 2);
1104
        $tF = tan(M_PI / 4 - $phiF / 2) / ((1 - $e * sin($phiF)) / (1 + $e * sin($phiF))) ** ($e / 2);
1105
        $n = (log($m1) - log($m2)) / (log($t1) - log($t2));
1106
        $F = $m1 / ($n * $t1 ** $n);
1107 36
        $rF = $a * $F * $tF ** $n;
1108 36
        $r = hypot($easting, $rF - $northing) * static::sign($n);
1109 36
        $t = ($r / ($a * $F)) ** (1 / $n);
1110 36
        if ($n >= 0) {
1111 36
            $theta = atan2($easting, $rF - $northing);
1112 36
        } else {
1113 36
            $theta = atan2(-$easting, -($rF - $northing));
1114 36
        }
1115 36
1116 36
        $latitude = M_PI / 2 - 2 * atan($t);
1117 36
        do {
1118
            $latitudeN = $latitude;
1119 36
            $latitude = M_PI / 2 - 2 * atan($t * ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2));
1120 36
        } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE_FORMULA);
1121 36
1122 36
        $longitude = $theta / $n + $lambdaOrigin;
1123 36
1124 36
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
1125 36
    }
1126 36
1127 36
    /**
1128 36
     * Lambert Conic Conformal (2SP Michigan).
1129 36
     */
1130 36
    public function lambertConicConformal2SPMichigan(
1131
        Geographic2D|Geographic3D $to,
1132
        Angle $latitudeOfFalseOrigin,
1133
        Angle $longitudeOfFalseOrigin,
1134
        Angle $latitudeOf1stStandardParallel,
1135 36
        Angle $latitudeOf2ndStandardParallel,
1136
        Length $eastingAtFalseOrigin,
1137 36
        Length $northingAtFalseOrigin,
1138 36
        Scale $ellipsoidScalingFactor
1139 36
    ): GeographicPoint {
1140
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1141 36
        $easting = $this->easting->asMetres()->getValue() - $eastingAtFalseOrigin->asMetres()->getValue();
1142
        $northing = $this->northing->asMetres()->getValue() - $northingAtFalseOrigin->asMetres()->getValue();
1143 36
        $lambdaOrigin = $longitudeOfFalseOrigin->asRadians()->getValue();
1144
        $phiF = $latitudeOfFalseOrigin->asRadians()->getValue();
1145
        $phi1 = $latitudeOf1stStandardParallel->asRadians()->getValue();
1146
        $phi2 = $latitudeOf2ndStandardParallel->asRadians()->getValue();
1147
        $K = $ellipsoidScalingFactor->asUnity()->getValue();
1148
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1149
        $e = $ellipsoid->getEccentricity();
1150
        $e2 = $ellipsoid->getEccentricitySquared();
1151 18
1152
        $m1 = cos($phi1) / sqrt(1 - $e2 * sin($phi1) ** 2);
1153
        $m2 = cos($phi2) / sqrt(1 - $e2 * sin($phi2) ** 2);
1154
        $t1 = tan(M_PI / 4 - $phi1 / 2) / ((1 - $e * sin($phi1)) / (1 + $e * sin($phi1))) ** ($e / 2);
1155
        $t2 = tan(M_PI / 4 - $phi2 / 2) / ((1 - $e * sin($phi2)) / (1 + $e * sin($phi2))) ** ($e / 2);
1156
        $tF = tan(M_PI / 4 - $phiF / 2) / ((1 - $e * sin($phiF)) / (1 + $e * sin($phiF))) ** ($e / 2);
1157
        $n = (log($m1) - log($m2)) / (log($t1) - log($t2));
1158
        $F = $m1 / ($n * $t1 ** $n);
1159
        $rF = $a * $K * $F * $tF ** $n;
1160 18
        $r = sqrt($easting ** 2 + ($rF - $northing) ** 2) * static::sign($n);
1161 18
        $t = ($r / ($a * $K * $F)) ** (1 / $n);
1162 18
        if ($n >= 0) {
1163 18
            $theta = atan2($easting, $rF - $northing);
1164 18
        } else {
1165 18
            $theta = atan2(-$easting, -($rF - $northing));
1166 18
        }
1167 18
1168 18
        $latitude = M_PI / 2 - 2 * atan($t);
1169 18
        do {
1170
            $latitudeN = $latitude;
1171 18
            $latitude = M_PI / 2 - 2 * atan($t * ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2));
1172 18
        } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE_FORMULA);
1173 18
1174 18
        $longitude = $theta / $n + $lambdaOrigin;
1175 18
1176 18
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
1177 18
    }
1178 18
1179 18
    /**
1180 9
     * Lambert Conic Conformal (2SP Belgium)
1181
     * In 2000 this modification was replaced through use of the regular Lambert Conic Conformal (2SP) method [9802]
1182 18
     * with appropriately modified parameter values.
1183 18
     */
1184 18
    public function lambertConicConformal2SPBelgium(
1185 18
        Geographic2D|Geographic3D $to,
1186
        Angle $latitudeOfFalseOrigin,
1187
        Angle $longitudeOfFalseOrigin,
1188
        Angle $latitudeOf1stStandardParallel,
1189
        Angle $latitudeOf2ndStandardParallel,
1190 18
        Length $eastingAtFalseOrigin,
1191
        Length $northingAtFalseOrigin
1192 18
    ): GeographicPoint {
1193 18
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1194 18
        $easting = $this->easting->asMetres()->getValue() - $eastingAtFalseOrigin->asMetres()->getValue();
1195
        $northing = $this->northing->asMetres()->getValue() - $northingAtFalseOrigin->asMetres()->getValue();
1196 18
        $lambdaOrigin = $longitudeOfFalseOrigin->asRadians()->getValue();
1197
        $phiF = $latitudeOfFalseOrigin->asRadians()->getValue();
1198 18
        $phi1 = $latitudeOf1stStandardParallel->asRadians()->getValue();
1199
        $phi2 = $latitudeOf2ndStandardParallel->asRadians()->getValue();
1200
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1201
        $e = $ellipsoid->getEccentricity();
1202
        $e2 = $ellipsoid->getEccentricitySquared();
1203
1204
        $m1 = cos($phi1) / sqrt(1 - $e2 * sin($phi1) ** 2);
1205
        $m2 = cos($phi2) / sqrt(1 - $e2 * sin($phi2) ** 2);
1206 18
        $t1 = tan(M_PI / 4 - $phi1 / 2) / ((1 - $e * sin($phi1)) / (1 + $e * sin($phi1))) ** ($e / 2);
1207
        $t2 = tan(M_PI / 4 - $phi2 / 2) / ((1 - $e * sin($phi2)) / (1 + $e * sin($phi2))) ** ($e / 2);
1208
        $tF = tan(M_PI / 4 - $phiF / 2) / ((1 - $e * sin($phiF)) / (1 + $e * sin($phiF))) ** ($e / 2);
1209
        $n = (log($m1) - log($m2)) / (log($t1) - log($t2));
1210
        $F = $m1 / ($n * $t1 ** $n);
1211
        $rF = $a * $F * $tF ** $n;
1212
        if (is_nan($rF)) {
1213
            $rF = 0;
1214 18
        }
1215 18
        $r = hypot($easting, $rF - $northing) * static::sign($n);
1216 18
        $t = ($r / ($a * $F)) ** (1 / $n);
1217 18
        if ($n >= 0) {
1218 18
            $theta = atan2($easting, $rF - $northing);
1219 18
        } else {
1220 18
            $theta = atan2(-$easting, -($rF - $northing));
1221 18
        }
1222 18
1223
        $latitude = M_PI / 2 - 2 * atan($t);
1224 18
        do {
1225 18
            $latitudeN = $latitude;
1226 18
            $latitude = M_PI / 2 - 2 * atan($t * ((1 - $e * sin($latitude)) / (1 + $e * sin($latitude))) ** ($e / 2));
1227 18
        } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE_FORMULA);
1228 18
1229 18
        $longitude = ($theta + (new ArcSecond(29.2985))->asRadians()->getValue()) / $n + $lambdaOrigin;
1230 18
1231 18
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
1232 18
    }
1233 18
1234 18
    /**
1235
     * Lambert Conic Near-Conformal
1236 18
     * The Lambert Near-Conformal projection is derived from the Lambert Conformal Conic projection by truncating the
1237 18
     * series expansion of the projection formulae.
1238 18
     */
1239
    public function lambertConicNearConformal(
1240 18
        Geographic2D|Geographic3D $to,
1241
        Angle $latitudeOfNaturalOrigin,
1242 18
        Angle $longitudeOfNaturalOrigin,
1243 18
        Scale $scaleFactorAtNaturalOrigin,
1244 18
        Length $falseEasting,
1245
        Length $falseNorthing
1246 18
    ): GeographicPoint {
1247
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1248 18
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
1249 18
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
1250 18
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
1251
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
1252 18
        $scaleFactorOrigin = $scaleFactorAtNaturalOrigin->asUnity()->getValue();
1253
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1254 18
        $e2 = $ellipsoid->getEccentricitySquared();
1255
        $f = $ellipsoid->getFlattening();
1256
1257
        $n = $f / (2 - $f);
1258
        $rhoO = $a * (1 - $e2) / (1 - $e2 * sin($latitudeOrigin) ** 2) ** (3 / 2);
1259
        $nuO = $a / sqrt(1 - $e2 * (sin($latitudeOrigin) ** 2));
1260
        $A = 1 / (6 * $rhoO * $nuO);
1261 18
        $APrime = $a * (1 - $n + 5 * ($n ** 2 - $n ** 3) / 4 + 81 * ($n ** 4 - $n ** 5) / 64);
1262
        $BPrime = 3 * $a * ($n - $n ** 2 + 7 * ($n ** 3 - $n ** 4) / 8 + 55 * $n ** 5 / 64) / 2;
1263
        $CPrime = 15 * $a * ($n ** 2 - $n ** 3 + 3 * ($n ** 4 - $n ** 5) / 4) / 16;
1264
        $DPrime = 35 * $a * ($n ** 3 - $n ** 4 + 11 * $n ** 5 / 16) / 48;
1265
        $EPrime = 315 * $a * ($n ** 4 - $n ** 5) / 512;
1266
        $rO = $scaleFactorOrigin * $nuO / tan($latitudeOrigin);
1267
        $sO = $APrime * $latitudeOrigin - $BPrime * sin(2 * $latitudeOrigin) + $CPrime * sin(4 * $latitudeOrigin) - $DPrime * sin(6 * $latitudeOrigin) + $EPrime * sin(8 * $latitudeOrigin);
1268 18
1269 18
        $theta = atan2($easting, $rO - $northing);
1270 18
        $r = hypot($easting, $rO - $northing) * static::sign($latitudeOrigin);
1271 18
        $M = $rO - $r;
1272 18
1273 18
        $m = $M;
1274 18
        do {
1275 18
            $mN = $m;
1276 18
            $m = $m - ($M - $scaleFactorOrigin * $m - $scaleFactorOrigin * $A * $m ** 3) / (-$scaleFactorOrigin - 3 * $scaleFactorOrigin * $A * $m ** 2);
1277 18
        } while (abs($m - $mN) >= static::ITERATION_CONVERGENCE_FORMULA);
1278
1279 18
        $latitude = $latitudeOrigin + $m / $A;
1280 18
        do {
1281 18
            $latitudeN = $latitude;
1282
            $latitude = $latitude + ($m + $sO - ($APrime * $latitude - $BPrime * sin(2 * $latitude) + $CPrime * sin(4 * $latitude) - $DPrime * sin(6 * $latitude) + $EPrime * sin(8 * $latitude))) / $APrime;
1283 18
        } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE_FORMULA);
1284 18
1285
        $longitude = $longitudeOrigin + $theta / sin($latitudeOrigin);
1286 18
1287
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
1288
    }
1289
1290
    /**
1291
     * Lambert Cylindrical Equal Area
1292
     * This is the ellipsoidal form of the projection.
1293
     */
1294 18
    public function lambertCylindricalEqualArea(
1295
        Geographic2D|Geographic3D $to,
1296
        Angle $latitudeOf1stStandardParallel,
1297
        Angle $longitudeOfNaturalOrigin,
1298
        Length $falseEasting,
1299
        Length $falseNorthing
1300
    ): GeographicPoint {
1301 18
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1302 18
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
1303 18
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
1304 18
        $latitudeFirstParallel = $latitudeOf1stStandardParallel->asRadians()->getValue();
1305 18
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
1306 18
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1307 18
        $e = $ellipsoid->getEccentricity();
1308
        $e2 = $ellipsoid->getEccentricitySquared();
1309 18
        $e4 = $e ** 4;
1310 18
        $e6 = $e ** 6;
1311 18
1312 18
        $k = cos($latitudeFirstParallel) / sqrt(1 - $e2 * sin($latitudeFirstParallel) ** 2);
1313 18
        $qP = (1 - $e2) * ((sin(M_PI_2) / (1 - $e2 * sin(M_PI_2) ** 2)) - (1 / (2 * $e)) * log((1 - $e * sin(M_PI_2)) / (1 + $e * sin(M_PI_2))));
1314 18
        $beta = self::asin(2 * $northing * $k / ($a * $qP));
1315 18
1316 18
        $latitude = $beta + (($e2 / 3 + 31 * $e4 / 180 + 517 * $e6 / 5040) * sin(2 * $beta)) + ((23 * $e4 / 360 + 251 * $e6 / 3780) * sin(4 * $beta)) + ((761 * $e6 / 45360) * sin(6 * $beta));
1317 18
        $longitude = $longitudeOrigin + $easting / ($a * $k);
1318
1319 18
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
1320 18
    }
1321
1322 18
    /**
1323
     * Modified Azimuthal Equidistant
1324
     * Modified form of Oblique Azimuthal Equidistant projection method developed for Polynesian islands. For the
1325
     * distances over which these projections are used (under 800km) this modification introduces no significant error.
1326
     */
1327
    public function modifiedAzimuthalEquidistant(
1328
        Geographic2D|Geographic3D $to,
1329
        Angle $latitudeOfNaturalOrigin,
1330 144
        Angle $longitudeOfNaturalOrigin,
1331
        Length $falseEasting,
1332
        Length $falseNorthing
1333
    ): GeographicPoint {
1334
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1335
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
1336
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
1337
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
1338 144
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
1339 144
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1340 144
        $e2 = $ellipsoid->getEccentricitySquared();
1341 144
1342 144
        $nuO = $a / sqrt(1 - $e2 * sin($latitudeOrigin) ** 2);
1343 144
        $c = hypot($easting, $northing);
1344 144
        $alpha = atan2($easting, $northing);
1345 144
        $A = -$e2 * cos($latitudeOrigin) ** 2 * cos($alpha) ** 2 / (1 - $e2);
1346 144
        $B = 3 * $e2 * (1 - $A) * sin($latitudeOrigin) * cos($latitudeOrigin) * cos($alpha) / (1 - $e2);
1347
        $D = $c / $nuO;
1348 144
        $J = $D - ($A * (1 + $A) * $D ** 3 / 6) - ($B * (1 + 3 * $A) * $D ** 4 / 24);
1349 144
        $K = 1 - ($A * $J ** 2 / 2) - ($B * $J ** 3 / 6);
1350 144
        $psi = self::asin(sin($latitudeOrigin) * cos($J) + cos($latitudeOrigin) * sin($J) * cos($alpha));
1351
1352 144
        $latitude = atan((1 - $e2 * $K * sin($latitudeOrigin) / sin($psi)) * tan($psi) / (1 - $e2));
1353 144
        $longitude = $longitudeOrigin + self::asin(sin($alpha) * sin($J) / cos($psi));
1354 144
1355 144
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
1356 144
    }
1357 144
1358 144
    /**
1359
     * Oblique Stereographic
1360 144
     * This is not the same as the projection method of the same name in USGS Professional Paper no. 1395, "Map
1361 144
     * Projections - A Working Manual" by John P. Snyder.
1362 144
     */
1363 144
    public function obliqueStereographic(
1364 144
        Geographic2D|Geographic3D $to,
1365 144
        Angle $latitudeOfNaturalOrigin,
1366
        Angle $longitudeOfNaturalOrigin,
1367 144
        Scale $scaleFactorAtNaturalOrigin,
1368
        Length $falseEasting,
1369 144
        Length $falseNorthing
1370
    ): GeographicPoint {
1371 144
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1372
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
1373 144
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
1374 144
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
1375 144
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
1376 144
        $scaleFactorOrigin = $scaleFactorAtNaturalOrigin->asUnity()->getValue();
1377
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1378 144
        $e = $ellipsoid->getEccentricity();
1379
        $e2 = $ellipsoid->getEccentricitySquared();
1380
1381
        $rhoOrigin = $a * (1 - $e2) / (1 - $e2 * sin($latitudeOrigin) ** 2) ** (3 / 2);
1382
        $nuOrigin = $a / sqrt(1 - $e2 * (sin($latitudeOrigin) ** 2));
1383
        $R = sqrt($rhoOrigin * $nuOrigin);
1384
1385 99
        $n = sqrt(1 + ($e2 * cos($latitudeOrigin) ** 4 / (1 - $e2)));
1386
        $S1 = (1 + sin($latitudeOrigin)) / (1 - sin($latitudeOrigin));
1387
        $S2 = (1 - $e * sin($latitudeOrigin)) / (1 + $e * sin($latitudeOrigin));
1388
        $w1 = ($S1 * ($S2 ** $e)) ** $n;
1389
        $c = (($n + sin($latitudeOrigin)) * (1 - ($w1 - 1) / ($w1 + 1))) / (($n - sin($latitudeOrigin)) * (1 + ($w1 - 1) / ($w1 + 1)));
1390
        $w2 = $c * $w1;
1391
        $chiOrigin = self::asin(($w2 - 1) / ($w2 + 1));
1392
1393 99
        $g = 2 * $R * $scaleFactorOrigin * tan(M_PI / 4 - $chiOrigin / 2);
1394 99
        $h = 4 * $R * $scaleFactorOrigin * tan($chiOrigin) + $g;
1395 99
        $i = atan2($easting, ($h + $northing));
1396 99
        $j = atan2($easting, ($g - $northing)) - $i;
1397 99
        $chi = $chiOrigin + 2 * atan(($northing - $easting * tan($j / 2)) / (2 * $R * $scaleFactorOrigin));
1398 99
        $lambda = $j + 2 * $i + $longitudeOrigin;
1399 99
1400 99
        $longitude = ($lambda - $longitudeOrigin) / $n + $longitudeOrigin;
1401 99
1402 99
        $psi = 0.5 * log((1 + sin($chi)) / ($c * (1 - sin($chi)))) / $n;
1403 99
1404 99
        $latitude = 2 * atan(M_E ** $psi) - M_PI / 2;
1405
        do {
1406 99
            $latitudeN = $latitude;
1407 99
            $psiN = log((tan($latitudeN / 2 + M_PI / 4)) * ((1 - $e * sin($latitudeN)) / (1 + $e * sin($latitudeN))) ** ($e / 2));
1408
            $latitude = $latitudeN - ($psiN - $psi) * cos($latitudeN) * (1 - $e2 * sin($latitudeN) ** 2) / (1 - $e2);
1409 99
        } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE_FORMULA);
1410 27
1411
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
1412 72
    }
1413
1414
    /**
1415 99
     * Polar Stereographic (variant A)
1416
     * Latitude of natural origin must be either 90 degrees or -90 degrees (or equivalent in alternative angle unit).
1417 99
     */
1418 8
    public function polarStereographicVariantA(
1419 91
        Geographic2D|Geographic3D $to,
1420 19
        Angle $latitudeOfNaturalOrigin,
1421
        Angle $longitudeOfNaturalOrigin,
1422 72
        Scale $scaleFactorAtNaturalOrigin,
1423
        Length $falseEasting,
1424
        Length $falseNorthing
1425 99
    ): GeographicPoint {
1426
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1427
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
1428
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
1429
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
1430
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
1431 243
        $scaleFactorOrigin = $scaleFactorAtNaturalOrigin->asUnity()->getValue();
1432
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1433
        $e = $ellipsoid->getEccentricity();
1434
        $e2 = $ellipsoid->getEccentricitySquared();
1435
        $e4 = $e ** 4;
1436
        $e6 = $e ** 6;
1437
        $e8 = $e ** 8;
1438 243
1439 243
        $rho = hypot($easting, $northing);
1440 243
        $t = $rho * sqrt((1 + $e) ** (1 + $e) * (1 - $e) ** (1 - $e)) / (2 * $a * $scaleFactorOrigin);
1441 243
1442 243
        if ($latitudeOrigin < 0) {
1443 243
            $chi = 2 * atan($t) - M_PI / 2;
1444 243
        } else {
1445 243
            $chi = M_PI / 2 - 2 * atan($t);
1446 243
        }
1447 243
1448 243
        $latitude = $chi + ($e2 / 2 + 5 * $e4 / 24 + $e6 / 12 + 13 * $e8 / 360) * sin(2 * $chi) + (7 * $e4 / 48 + 29 * $e6 / 240 + 811 * $e8 / 11520) * sin(4 * $chi) + (7 * $e6 / 120 + 81 * $e8 / 1120) * sin(6 * $chi) + (4279 * $e8 / 161280) * sin(8 * $chi);
1449
1450 243
        if ($easting === 0.0) {
0 ignored issues
show
introduced by
The condition $easting === 0.0 is always false.
Loading history...
1451 243
            $longitude = $longitudeOrigin;
1452 216
        } elseif ($latitudeOrigin < 0) {
1453
            $longitude = $longitudeOrigin + atan2($easting, $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue());
1454 27
        } else {
1455
            $longitude = $longitudeOrigin + atan2($easting, $falseNorthing->asMetres()->getValue() - $this->northing->asMetres()->getValue());
1456 243
        }
1457 243
1458 243
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
1459
    }
1460 243
1461 216
    /**
1462
     * Polar Stereographic (variant B).
1463 27
     */
1464
    public function polarStereographicVariantB(
1465
        Geographic2D|Geographic3D $to,
1466 243
        Angle $latitudeOfStandardParallel,
1467
        Angle $longitudeOfOrigin,
1468 243
        Length $falseEasting,
1469 161
        Length $falseNorthing
1470 82
    ): GeographicPoint {
1471 55
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1472
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
1473 27
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
1474
        $standardParallel = $latitudeOfStandardParallel->asRadians()->getValue();
1475
        $longitudeOrigin = $longitudeOfOrigin->asRadians()->getValue();
1476 243
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1477
        $e = $ellipsoid->getEccentricity();
1478
        $e2 = $ellipsoid->getEccentricitySquared();
1479
        $e4 = $e ** 4;
1480
        $e6 = $e ** 6;
1481
        $e8 = $e ** 8;
1482 27
1483
        $rho = hypot($easting, $northing);
1484
        if ($standardParallel < 0) {
1485
            $tF = tan(M_PI / 4 + $standardParallel / 2) / (((1 + $e * sin($standardParallel)) / (1 - $e * sin($standardParallel))) ** ($e / 2));
1486
        } else {
1487
            $tF = tan(M_PI / 4 - $standardParallel / 2) * (((1 + $e * sin($standardParallel)) / (1 - $e * sin($standardParallel))) ** ($e / 2));
1488
        }
1489 27
        $mF = cos($standardParallel) / sqrt(1 - $e2 * sin($standardParallel) ** 2);
1490 27
        $kO = $mF * sqrt((1 + $e) ** (1 + $e) * (1 - $e) ** (1 - $e)) / (2 * $tF);
1491 27
        $t = $rho * sqrt((1 + $e) ** (1 + $e) * (1 - $e) ** (1 - $e)) / (2 * $a * $kO);
1492 27
1493 27
        if ($standardParallel < 0) {
1494 27
            $chi = 2 * atan($t) - M_PI / 2;
1495 27
        } else {
1496 27
            $chi = M_PI / 2 - 2 * atan($t);
1497 27
        }
1498 27
1499 27
        $latitude = $chi + ($e2 / 2 + 5 * $e4 / 24 + $e6 / 12 + 13 * $e8 / 360) * sin(2 * $chi) + (7 * $e4 / 48 + 29 * $e6 / 240 + 811 * $e8 / 11520) * sin(4 * $chi) + (7 * $e6 / 120 + 81 * $e8 / 1120) * sin(6 * $chi) + (4279 * $e8 / 161280) * sin(8 * $chi);
1500
1501 27
        if ($easting === 0.0) {
0 ignored issues
show
introduced by
The condition $easting === 0.0 is always false.
Loading history...
1502 27
            $longitude = $longitudeOrigin;
1503
        } elseif ($standardParallel < 0) {
1504
            $longitude = $longitudeOrigin + atan2($easting, $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue());
1505
        } else {
1506 27
            $longitude = $longitudeOrigin + atan2($easting, $falseNorthing->asMetres()->getValue() - $this->northing->asMetres()->getValue());
1507 27
        }
1508 27
1509 27
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
1510 27
    }
1511 27
1512
    /**
1513
     * Polar Stereographic (variant C).
1514
     */
1515
    public function polarStereographicVariantC(
1516
        Geographic2D|Geographic3D $to,
1517
        Angle $latitudeOfStandardParallel,
1518 27
        Angle $longitudeOfOrigin,
1519
        Length $eastingAtFalseOrigin,
1520 27
        Length $northingAtFalseOrigin
1521
    ): GeographicPoint {
1522 27
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1523 27
        $easting = $this->easting->asMetres()->getValue() - $eastingAtFalseOrigin->asMetres()->getValue();
1524
        $northing = $this->northing->asMetres()->getValue() - $northingAtFalseOrigin->asMetres()->getValue();
1525
        $standardParallel = $latitudeOfStandardParallel->asRadians()->getValue();
1526
        $longitudeOrigin = $longitudeOfOrigin->asRadians()->getValue();
1527
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1528 27
        $e = $ellipsoid->getEccentricity();
1529
        $e2 = $ellipsoid->getEccentricitySquared();
1530
        $e4 = $e ** 4;
1531
        $e6 = $e ** 6;
1532
        $e8 = $e ** 8;
1533
1534
        if ($standardParallel < 0) {
1535 27
            $tF = tan(M_PI / 4 + $standardParallel / 2) / (((1 + $e * sin($standardParallel)) / (1 - $e * sin($standardParallel))) ** ($e / 2));
1536
        } else {
1537
            $tF = tan(M_PI / 4 - $standardParallel / 2) * (((1 + $e * sin($standardParallel)) / (1 - $e * sin($standardParallel))) ** ($e / 2));
1538
        }
1539
        $mF = cos($standardParallel) / sqrt(1 - $e2 * sin($standardParallel) ** 2);
1540
        $rhoF = $a * $mF;
1541
        if ($standardParallel < 0) {
1542 27
            $rho = hypot($easting, $northing + $rhoF);
1543 27
            $t = $rho * $tF / $rhoF;
1544 27
            $chi = 2 * atan($t) - M_PI / 2;
1545 27
        } else {
1546 27
            $rho = hypot($easting, $northing - $rhoF);
1547 27
            $t = $rho * $tF / $rhoF;
1548
            $chi = M_PI / 2 - 2 * atan($t);
1549 27
        }
1550 27
1551 27
        $latitude = $chi + ($e2 / 2 + 5 * $e4 / 24 + $e6 / 12 + 13 * $e8 / 360) * sin(2 * $chi) + (7 * $e4 / 48 + 29 * $e6 / 240 + 811 * $e8 / 11520) * sin(4 * $chi) + (7 * $e6 / 120 + 81 * $e8 / 1120) * sin(6 * $chi) + (4279 * $e8 / 161280) * sin(8 * $chi);
1552
1553 27
        if ($easting === 0.0) {
0 ignored issues
show
introduced by
The condition $easting === 0.0 is always false.
Loading history...
1554
            $longitude = $longitudeOrigin;
1555
        } elseif ($standardParallel < 0) {
1556
            $longitude = $longitudeOrigin + atan2($easting, $this->northing->asMetres()->getValue() - $northingAtFalseOrigin->asMetres()->getValue() + $rhoF);
1557
        } else {
1558
            $longitude = $longitudeOrigin + atan2($easting, $northingAtFalseOrigin->asMetres()->getValue() - $this->northing->asMetres()->getValue() + $rhoF);
1559
        }
1560 234
1561
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
1562
    }
1563
1564
    /**
1565
     * Popular Visualisation Pseudo Mercator
1566
     * Applies spherical formulas to the ellipsoid. As such does not have the properties of a true Mercator projection.
1567
     */
1568 234
    public function popularVisualisationPseudoMercator(
1569 234
        Geographic2D|Geographic3D $to,
1570 234
        Angle $latitudeOfNaturalOrigin,
1571 234
        Angle $longitudeOfNaturalOrigin,
1572 234
        Length $falseEasting,
1573 234
        Length $falseNorthing
1574
    ): GeographicPoint {
1575 234
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1576 225
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
1577 225
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
1578
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
0 ignored issues
show
Unused Code introduced by
The assignment to $latitudeOrigin is dead and can be removed.
Loading history...
1579 234
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
1580 234
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1581
1582
        $D = -$northing / $a;
1583 234
        $latitude = M_PI / 2 - 2 * atan(M_E ** $D);
1584
        $longitude = $easting / $a + $longitudeOrigin;
1585
1586
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
1587
    }
1588
1589
    /**
1590
     * Similarity transformation
1591
     * Defined for two-dimensional coordinate systems.
1592 90
     */
1593
    public function similarityTransformation(
1594
        Projected $to,
1595
        Length $ordinate1OfEvaluationPointInTargetCRS,
1596
        Length $ordinate2OfEvaluationPointInTargetCRS,
1597
        Scale $scaleFactorForSourceCRSAxes,
1598
        Angle $rotationAngleOfSourceCRSAxes,
1599
        bool $inReverse
1600 90
    ): self {
1601 90
        $xs = $this->easting->asMetres()->getValue();
1602 90
        $ys = $this->northing->asMetres()->getValue();
1603 90
        $xo = $ordinate1OfEvaluationPointInTargetCRS->asMetres()->getValue();
1604 90
        $yo = $ordinate2OfEvaluationPointInTargetCRS->asMetres()->getValue();
1605 90
        $M = $scaleFactorForSourceCRSAxes->asUnity()->getValue();
1606 90
        $theta = $rotationAngleOfSourceCRSAxes->asRadians()->getValue();
1607 90
1608 90
        if ($inReverse) {
1609 90
            $easting = (($xs - $xo) * cos($theta) - ($ys - $yo) * sin($theta)) / $M;
1610 90
            $northing = (($xs - $xo) * sin($theta) + ($ys - $yo) * cos($theta)) / $M;
1611 90
        } else {
1612
            $easting = $xo + $xs * $M * cos($theta) + $ys * $M * sin($theta);
1613 90
            $northing = $yo - $xs * $M * sin($theta) + $ys * $M * cos($theta);
1614 90
        }
1615
1616 90
        return self::create($to, new Metre($easting), new Metre($northing), new Metre(-$easting), new Metre(-$northing), $this->epoch);
1617 90
    }
1618
1619 90
    /**
1620
     * Mercator (variant A)
1621
     * Note that in these formulas the parameter latitude of natural origin (latO) is not used. However for this
1622
     * Mercator (variant A) method the EPSG dataset includes this parameter, which must have a value of zero, for
1623
     * completeness in CRS labelling.
1624
     */
1625
    public function mercatorVariantA(
1626 36
        Geographic2D|Geographic3D $to,
1627
        Angle $latitudeOfNaturalOrigin,
1628
        Angle $longitudeOfNaturalOrigin,
1629
        Scale $scaleFactorAtNaturalOrigin,
1630
        Length $falseEasting,
1631
        Length $falseNorthing
1632
    ): GeographicPoint {
1633 36
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1634 36
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
1635 36
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
0 ignored issues
show
Unused Code introduced by
The assignment to $northing is dead and can be removed.
Loading history...
1636 36
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
0 ignored issues
show
Unused Code introduced by
The assignment to $latitudeOrigin is dead and can be removed.
Loading history...
1637 36
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
1638 36
        $scaleFactorOrigin = $scaleFactorAtNaturalOrigin->asUnity()->getValue();
1639 36
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1640 36
        $e = $ellipsoid->getEccentricity();
1641 36
        $e2 = $ellipsoid->getEccentricitySquared();
1642 36
        $e4 = $e ** 4;
1643 36
        $e6 = $e ** 6;
1644
        $e8 = $e ** 8;
1645 36
1646
        $t = M_E ** (($falseNorthing->asMetres()->getValue() - $this->northing->asMetres()->getValue()) / ($a * $scaleFactorOrigin));
1647 36
        $chi = M_PI / 2 - 2 * atan($t);
1648 36
1649
        $latitude = $chi + ($e2 / 2 + 5 * $e4 / 24 + $e6 / 12 + 13 * $e8 / 360) * sin(2 * $chi) + (7 * $e4 / 48 + 29 * $e6 / 240 + 811 * $e8 / 11520) * sin(4 * $chi) + (7 * $e6 / 120 + 81 * $e8 / 1120) * sin(6 * $chi) + (4279 * $e8 / 161280) * sin(8 * $chi);
1650 36
        $longitude = $easting / ($a * $scaleFactorOrigin) + $longitudeOrigin;
1651 36
1652
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
1653 36
    }
1654
1655
    /**
1656
     * Mercator (variant B)
1657
     * Used for most nautical charts.
1658
     */
1659 207
    public function mercatorVariantB(
1660
        Geographic2D|Geographic3D $to,
1661
        Angle $latitudeOf1stStandardParallel,
1662
        Angle $longitudeOfNaturalOrigin,
1663
        Length $falseEasting,
1664
        Length $falseNorthing
1665
    ): GeographicPoint {
1666
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1667
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
1668
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
0 ignored issues
show
Unused Code introduced by
The assignment to $northing is dead and can be removed.
Loading history...
1669 207
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
1670 207
        $firstStandardParallel = $latitudeOf1stStandardParallel->asRadians()->getValue();
1671 207
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1672 207
        $e = $ellipsoid->getEccentricity();
1673 207
        $e2 = $ellipsoid->getEccentricitySquared();
1674 207
        $e4 = $e ** 4;
1675 207
        $e6 = $e ** 6;
1676 207
        $e8 = $e ** 8;
1677 207
1678 207
        $scaleFactorOrigin = cos($firstStandardParallel) / sqrt(1 - $e2 * sin($firstStandardParallel) ** 2);
1679 207
1680 207
        $t = M_E ** (($falseNorthing->asMetres()->getValue() - $this->northing->asMetres()->getValue()) / ($a * $scaleFactorOrigin));
1681 207
        $chi = M_PI / 2 - 2 * atan($t);
1682 207
1683
        $latitude = $chi + ($e2 / 2 + 5 * $e4 / 24 + $e6 / 12 + 13 * $e8 / 360) * sin(2 * $chi) + (7 * $e4 / 48 + 29 * $e6 / 240 + 811 * $e8 / 11520) * sin(4 * $chi) + (7 * $e6 / 120 + 81 * $e8 / 1120) * sin(6 * $chi) + (4279 * $e8 / 161280) * sin(8 * $chi);
1684 207
        $longitude = $easting / ($a * $scaleFactorOrigin) + $longitudeOrigin;
1685 207
1686 207
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
1687 207
    }
1688 207
1689 207
    /**
1690 207
     * Hotine Oblique Mercator (variant A).
1691 207
     */
1692 207
    public function obliqueMercatorHotineVariantA(
1693 207
        Geographic2D|Geographic3D $to,
1694
        Angle $latitudeOfProjectionCentre,
1695 207
        Angle $longitudeOfProjectionCentre,
1696 207
        Angle $azimuthOfInitialLine,
1697
        Angle $angleFromRectifiedToSkewGrid,
1698 207
        Scale $scaleFactorOnInitialLine,
1699 207
        Length $falseEasting,
1700 207
        Length $falseNorthing
1701 207
    ): GeographicPoint {
1702 207
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1703 207
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
1704
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
1705 207
        $latC = $latitudeOfProjectionCentre->asRadians()->getValue();
1706
        $lonC = $longitudeOfProjectionCentre->asRadians()->getValue();
1707 207
        $alphaC = $azimuthOfInitialLine->asRadians()->getValue();
1708 207
        $kC = $scaleFactorOnInitialLine->asUnity()->getValue();
1709
        $gammaC = $angleFromRectifiedToSkewGrid->asRadians()->getValue();
1710 207
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1711
        $e = $ellipsoid->getEccentricity();
1712
        $e2 = $ellipsoid->getEccentricitySquared();
1713
        $e4 = $e ** 4;
1714
        $e6 = $e ** 6;
1715
        $e8 = $e ** 8;
1716 108
1717
        $B = sqrt(1 + ($e2 * cos($latC) ** 4 / (1 - $e2)));
1718
        $A = $a * $B * $kC * sqrt(1 - $e2) / (1 - $e2 * sin($latC) ** 2);
1719
        $tO = tan(M_PI / 4 - $latC / 2) / ((1 - $e * sin($latC)) / (1 + $e * sin($latC))) ** ($e / 2);
1720
        $D = $B * sqrt((1 - $e2)) / (cos($latC) * sqrt(1 - $e2 * sin($latC) ** 2));
1721
        $DD = max(1, $D ** 2);
1722
        $F = $D + sqrt($DD - 1) * static::sign($latC);
1723
        $H = $F * ($tO) ** $B;
1724
        $G = ($F - 1 / $F) / 2;
1725
        $gammaO = self::asin(sin($alphaC) / $D);
1726 108
        $lonO = $lonC - (self::asin($G * tan($gammaO))) / $B;
1727 108
1728 108
        $v = $easting * cos($gammaC) - $northing * sin($gammaC);
1729 108
        $u = $northing * cos($gammaC) + $easting * sin($gammaC);
1730 108
1731 108
        $Q = M_E ** -($B * $v / $A);
1732 108
        $S = ($Q - 1 / $Q) / 2;
1733 108
        $T = ($Q + 1 / $Q) / 2;
1734 108
        $V = sin($B * $u / $A);
1735 108
        $U = ($V * cos($gammaO) + $S * sin($gammaO)) / $T;
1736 108
        $t = ($H / sqrt((1 + $U) / (1 - $U))) ** (1 / $B);
1737 108
1738 108
        $chi = M_PI / 2 - 2 * atan($t);
1739 108
1740
        $latitude = $chi + sin(2 * $chi) * ($e2 / 2 + 5 * $e4 / 24 + $e6 / 12 + 13 * $e8 / 360) + sin(4 * $chi) * (7 * $e4 / 48 + 29 * $e6 / 240 + 811 * $e8 / 11520) + sin(6 * $chi) * (7 * $e6 / 120 + 81 * $e8 / 1120) + sin(8 * $chi) * (4279 * $e8 / 161280);
1741 108
        $longitude = $lonO - atan2(($S * cos($gammaO) - $V * sin($gammaO)), cos($B * $u / $A)) / $B;
1742 108
1743 108
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
1744 108
    }
1745 108
1746 108
    /**
1747 108
     * Hotine Oblique Mercator (variant B).
1748 108
     */
1749 108
    public function obliqueMercatorHotineVariantB(
1750 108
        Geographic2D|Geographic3D $to,
1751 108
        Angle $latitudeOfProjectionCentre,
1752 108
        Angle $longitudeOfProjectionCentre,
1753 45
        Angle $azimuthOfInitialLine,
1754
        Angle $angleFromRectifiedToSkewGrid,
1755 63
        Scale $scaleFactorOnInitialLine,
1756
        Length $eastingAtProjectionCentre,
1757
        Length $northingAtProjectionCentre
1758 108
    ): GeographicPoint {
1759 108
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1760
        $easting = $this->easting->asMetres()->getValue() - $eastingAtProjectionCentre->asMetres()->getValue();
1761 108
        $northing = $this->northing->asMetres()->getValue() - $northingAtProjectionCentre->asMetres()->getValue();
1762 108
        $latC = $latitudeOfProjectionCentre->asRadians()->getValue();
1763 108
        $lonC = $longitudeOfProjectionCentre->asRadians()->getValue();
1764 108
        $alphaC = $azimuthOfInitialLine->asRadians()->getValue();
1765 108
        $kC = $scaleFactorOnInitialLine->asUnity()->getValue();
1766 108
        $gammaC = $angleFromRectifiedToSkewGrid->asRadians()->getValue();
1767
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1768 108
        $e = $ellipsoid->getEccentricity();
1769
        $e2 = $ellipsoid->getEccentricitySquared();
1770 108
        $e4 = $e ** 4;
1771 108
        $e6 = $e ** 6;
1772
        $e8 = $e ** 8;
1773 108
1774
        $B = sqrt(1 + ($e2 * cos($latC) ** 4 / (1 - $e2)));
1775
        $A = $a * $B * $kC * sqrt(1 - $e2) / (1 - $e2 * sin($latC) ** 2);
1776
        $tO = tan(M_PI / 4 - $latC / 2) / ((1 - $e * sin($latC)) / (1 + $e * sin($latC))) ** ($e / 2);
1777
        $D = $B * sqrt((1 - $e2)) / (cos($latC) * sqrt(1 - $e2 * sin($latC) ** 2));
1778
        $DD = max(1, $D ** 2);
1779 27
        $F = $D + sqrt($DD - 1) * static::sign($latC);
1780
        $H = $F * ($tO) ** $B;
1781
        $G = ($F - 1 / $F) / 2;
1782
        $gammaO = self::asin(sin($alphaC) / $D);
1783
        $lonO = $lonC - (self::asin($G * tan($gammaO))) / $B;
1784
        $vC = 0;
0 ignored issues
show
Unused Code introduced by
The assignment to $vC is dead and can be removed.
Loading history...
1785
        if ($alphaC === M_PI / 2) {
1786
            $uC = $A * ($lonC - $lonO);
1787
        } else {
1788 27
            $uC = ($A / $B) * atan2(sqrt($DD - 1), cos($alphaC)) * static::sign($latC);
1789 27
        }
1790 27
1791 27
        $v = $easting * cos($gammaC) - $northing * sin($gammaC);
1792 27
        $u = $northing * cos($gammaC) + $easting * sin($gammaC) + (abs($uC) * static::sign($latC));
1793 27
1794 27
        $Q = M_E ** -($B * $v / $A);
1795 27
        $S = ($Q - 1 / $Q) / 2;
1796 27
        $T = ($Q + 1 / $Q) / 2;
1797 27
        $V = sin($B * $u / $A);
1798
        $U = ($V * cos($gammaO) + $S * sin($gammaO)) / $T;
1799 27
        $t = ($H / sqrt((1 + $U) / (1 - $U))) ** (1 / $B);
1800 27
1801 27
        $chi = M_PI / 2 - 2 * atan($t);
1802 27
1803
        $latitude = $chi + sin(2 * $chi) * ($e2 / 2 + 5 * $e4 / 24 + $e6 / 12 + 13 * $e8 / 360) + sin(4 * $chi) * (7 * $e4 / 48 + 29 * $e6 / 240 + 811 * $e8 / 11520) + sin(6 * $chi) * (7 * $e6 / 120 + 81 * $e8 / 1120) + sin(8 * $chi) * (4279 * $e8 / 161280);
1804 27
        $longitude = $lonO - atan2(($S * cos($gammaO) - $V * sin($gammaO)), cos($B * $u / $A)) / $B;
1805
1806 27
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
1807 27
    }
1808
1809 27
    /**
1810 27
     * Laborde Oblique Mercator.
1811 27
     */
1812
    public function obliqueMercatorLaborde(
1813 27
        Geographic2D|Geographic3D $to,
1814 27
        Angle $latitudeOfProjectionCentre,
1815 27
        Angle $longitudeOfProjectionCentre,
1816 27
        Angle $azimuthOfInitialLine,
1817 27
        Scale $scaleFactorOnInitialLine,
1818
        Length $falseEasting,
1819 27
        Length $falseNorthing
1820 27
    ): GeographicPoint {
1821
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1822
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
1823
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
1824 27
        $latC = $latitudeOfProjectionCentre->asRadians()->getValue();
1825 27
        $lonC = $longitudeOfProjectionCentre->asRadians()->getValue();
1826
        $alphaC = $azimuthOfInitialLine->asRadians()->getValue();
1827
        $kC = $scaleFactorOnInitialLine->asUnity()->getValue();
1828 27
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1829
        $e = $ellipsoid->getEccentricity();
1830 27
        $e2 = $ellipsoid->getEccentricitySquared();
1831
1832 27
        $B = sqrt(1 + ($e2 * cos($latC) ** 4 / (1 - $e2)));
1833
        $latS = self::asin(sin($latC) / $B);
1834 27
        $R = $a * $kC * (sqrt(1 - $e2) / (1 - $e2 * sin($latC) ** 2));
1835 27
        $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));
1836 27
1837
        $G = (new ComplexNumber(1 - cos(2 * $alphaC), sin(2 * $alphaC)))->divide(new ComplexNumber(12, 0));
1838 27
1839
        $H0 = new ComplexNumber($northing / $R, $easting / $R);
1840
        $H = $H0->divide($H0->pow(3)->multiply($G)->add($H0));
1841
        do {
1842
            $HN = $H;
1843
            $H = ($HN->pow(3)->multiply($G)->multiply(new ComplexNumber(2, 0))->add($H0))->divide($HN->pow(2)->multiply($G)->multiply(new ComplexNumber(3, 0))->add(new ComplexNumber(1, 0)));
1844 32051
        } while (abs($H0->subtract($H)->subtract($H->pow(3)->multiply($G))->getReal()) >= static::ITERATION_CONVERGENCE_FORMULA);
1845
1846
        $LPrime = -1 * $H->getReal();
1847
        $PPrime = 2 * atan(M_E ** $H->getImaginary()) - M_PI / 2;
1848
        $U = cos($PPrime) * cos($LPrime) * cos($latS) + cos($PPrime) * sin($LPrime) * sin($latS);
1849
        $V = sin($PPrime);
1850
        $W = cos($PPrime) * cos($LPrime) * sin($latS) - cos($PPrime) * sin($LPrime) * cos($latS);
1851
1852 32051
        $d = hypot($U, $V);
1853 32051
        if ($d === 0) {
1854 32051
            $L = 0;
1855 32051
            $P = static::sign($W) * M_PI / 2;
1856 32051
        } else {
1857 32051
            $L = 2 * atan($V / ($U + $d));
1858 32051
            $P = atan($W / $d);
1859 32051
        }
1860 32051
1861
        $longitude = $lonC + ($L / $B);
1862 32051
1863 32051
        $q = (log(tan(M_PI / 4 + $P / 2)) - $C) / $B;
1864
1865 32051
        $latitude = 2 * atan(M_E ** $q) - M_PI / 2;
1866 32051
        do {
1867 32051
            $latitudeN = $latitude;
1868 32051
            $latitude = 2 * atan(((1 + $e * sin($latitude)) / (1 - $e * sin($latitude))) ** ($e / 2) * M_E ** $q) - M_PI / 2;
1869
        } while (abs($latitude - $latitudeN) >= static::ITERATION_CONVERGENCE_FORMULA);
1870 32051
1871 21906
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), null, $this->epoch);
1872 10145
    }
1873
1874 10145
    /**
1875 333
     * Transverse Mercator.
1876
     */
1877 9812
    public function transverseMercator(
1878 9812
        Geographic2D|Geographic3D $to,
1879 9812
        Angle $latitudeOfNaturalOrigin,
1880 9812
        Angle $longitudeOfNaturalOrigin,
1881 9812
        Scale $scaleFactorAtNaturalOrigin,
1882 9812
        Length $falseEasting,
1883 9812
        Length $falseNorthing
1884 9812
    ): GeographicPoint {
1885 9812
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
1886
        $easting = $this->easting->asMetres()->getValue() - $falseEasting->asMetres()->getValue();
1887
        $northing = $this->northing->asMetres()->getValue() - $falseNorthing->asMetres()->getValue();
1888 32051
        $latitudeOrigin = $latitudeOfNaturalOrigin->asRadians()->getValue();
1889 32051
        $longitudeOrigin = $longitudeOfNaturalOrigin->asRadians()->getValue();
1890 32051
        $kO = $scaleFactorAtNaturalOrigin->asUnity()->getValue();
1891 32051
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
1892 32051
        $e = $ellipsoid->getEccentricity();
1893 32051
        $f = $ellipsoid->getFlattening();
1894 32051
1895 32051
        $n = $f / (2 - $f);
1896 32051
        $B = ($a / (1 + $n)) * (1 + $n ** 2 / 4 + $n ** 4 / 64);
1897 32051
1898 32051
        $h1 = $n / 2 - (2 / 3) * $n ** 2 + (37 / 96) * $n ** 3 - (1 / 360) * $n ** 4;
1899 32051
        $h2 = (1 / 48) * $n ** 2 + (1 / 15) * $n ** 3 - (437 / 1440) * $n ** 4;
1900
        $h3 = (17 / 480) * $n ** 3 - (37 / 840) * $n ** 4;
1901 32051
        $h4 = (4397 / 161280) * $n ** 4;
1902
1903 32051
        if ($latitudeOrigin === 0.0) {
0 ignored issues
show
introduced by
The condition $latitudeOrigin === 0.0 is always false.
Loading history...
1904 32051
            $mO = 0;
1905
        } elseif ($latitudeOrigin === M_PI / 2) {
1906 32051
            $mO = $B * M_PI / 2;
1907 32051
        } elseif ($latitudeOrigin === -M_PI / 2) {
1908 32051
            $mO = $B * -M_PI / 2;
1909
        } else {
1910 32051
            $qO = asinh(tan($latitudeOrigin)) - ($e * atanh($e * sin($latitudeOrigin)));
1911 32051
            $betaO = atan(sinh($qO));
1912
            $xiO0 = self::asin(sin($betaO));
1913 32051
            $xiO1 = $h1 * sin(2 * $xiO0);
1914
            $xiO2 = $h2 * sin(4 * $xiO0);
1915
            $xiO3 = $h3 * sin(6 * $xiO0);
1916
            $xiO4 = $h4 * sin(8 * $xiO0);
1917
            $xiO = $xiO0 + $xiO1 + $xiO2 + $xiO3 + $xiO4;
1918
            $mO = $B * $xiO;
1919
        }
1920
1921 36
        $eta = $easting / ($B * $kO);
1922
        $xi = ($northing + $kO * $mO) / ($B * $kO);
1923
        $xi1 = $h1 * sin(2 * $xi) * cosh(2 * $eta);
1924
        $eta1 = $h1 * cos(2 * $xi) * sinh(2 * $eta);
1925
        $xi2 = $h2 * sin(4 * $xi) * cosh(4 * $eta);
1926
        $eta2 = $h2 * cos(4 * $xi) * sinh(4 * $eta);
1927
        $xi3 = $h3 * sin(6 * $xi) * cosh(6 * $eta);
1928
        $eta3 = $h3 * cos(6 * $xi) * sinh(6 * $eta);
1929
        $xi4 = $h4 * sin(8 * $xi) * cosh(8 * $eta);
1930 36
        $eta4 = $h4 * cos(8 * $xi) * sinh(8 * $eta);
1931 36
        $xi0 = $xi - ($xi1 + $xi2 + $xi3 + $xi4);
1932
        $eta0 = $eta - ($eta1 + $eta2 + $eta3 + $eta4);
1933 36
1934 36
        $beta = self::asin(sin($xi0) / cosh($eta0));
1935
1936 36
        $QPrime = asinh(tan($beta));
1937
        $Q = asinh(tan($beta));
1938
        do {
1939
            $QN = $Q;
1940
            $Q = $QPrime + ($e * atanh($e * tanh($Q)));
1941
        } while (abs($Q - $QN) >= static::ITERATION_CONVERGENCE_FORMULA);
1942
1943 18
        $latitude = atan(sinh($Q));
1944
        $longitude = $longitudeOrigin + self::asin(tanh($eta0) / cos($beta));
1945
1946
        $height = $this->height && $to instanceof Geographic3D ? $this->height : null;
1947
1948
        return GeographicPoint::create($to, new Radian($latitude), new Radian($longitude), $height, $this->epoch);
1949
    }
1950
1951
    /**
1952
     * Transverse Mercator Zoned Grid System
1953
     * If locations fall outwith the fixed zones the general Transverse Mercator method (code 9807) must be used for
1954
     * each zone.
1955 18
     */
1956 18
    public function transverseMercatorZonedGrid(
1957
        Geographic2D|Geographic3D $to,
1958 18
        Angle $latitudeOfNaturalOrigin,
1959 18
        Angle $initialLongitude,
1960
        Angle $zoneWidth,
1961
        Scale $scaleFactorAtNaturalOrigin,
1962
        Length $falseEasting,
1963
        Length $falseNorthing
1964
    ): GeographicPoint {
1965
        $Z = substr((string) $this->easting->asMetres()->getValue(), 0, 2);
1966
        $falseEasting = $falseEasting->add(new Metre($Z * 1000000));
1967
1968
        $W = $zoneWidth->asDegrees()->getValue();
1969
        $longitudeOrigin = $initialLongitude->add(new Degree($Z * $W - $W / 2));
1970
1971
        return $this->transverseMercator($to, $latitudeOfNaturalOrigin, $longitudeOrigin, $scaleFactorAtNaturalOrigin, $falseEasting, $falseNorthing);
1972 18
    }
1973 18
1974
    /**
1975 18
     * General polynomial.
1976 18
     * @param Coefficient[] $powerCoefficients
1977 18
     */
1978
    public function generalPolynomial(
1979 18
        Projected $to,
1980
        Length $ordinate1OfEvaluationPointInSourceCRS,
1981
        Length $ordinate2OfEvaluationPointInSourceCRS,
1982
        Length $ordinate1OfEvaluationPointInTargetCRS,
1983
        Length $ordinate2OfEvaluationPointInTargetCRS,
1984
        Scale $scalingFactorForSourceCRSCoordDifferences,
1985
        Scale $scalingFactorForTargetCRSCoordDifferences,
1986 36
        Scale $A0,
1987
        Scale $B0,
1988
        array $powerCoefficients
1989
    ): self {
1990
        $xs = $this->easting->getValue();
1991
        $ys = $this->northing->getValue();
1992
1993 36
        $t = $this->generalPolynomialUnitless(
1994 36
            $xs,
1995
            $ys,
1996 36
            $ordinate1OfEvaluationPointInSourceCRS,
1997 36
            $ordinate2OfEvaluationPointInSourceCRS,
1998 36
            $ordinate1OfEvaluationPointInTargetCRS,
1999
            $ordinate2OfEvaluationPointInTargetCRS,
2000
            $scalingFactorForSourceCRSCoordDifferences,
2001 36
            $scalingFactorForTargetCRSCoordDifferences,
2002 36
            $A0,
2003 36
            $B0,
2004 36
            $powerCoefficients
2005 36
        );
2006 36
2007 36
        $xtUnit = $to->getCoordinateSystem()->getAxes()[0]->getUnitOfMeasureId();
2008 36
        $ytUnit = $to->getCoordinateSystem()->getAxes()[1]->getUnitOfMeasureId();
2009 36
2010 36
        return static::createFromEastingNorthing(
2011 36
            $to,
2012 36
            Length::makeUnit($t['xt'], $xtUnit),
2013
            Length::makeUnit($t['yt'], $ytUnit),
2014 36
            $this->epoch
2015 36
        );
2016 36
    }
2017 36
2018 36
    /**
2019 36
     * New Zealand Map Grid.
2020 36
     */
2021
    public function newZealandMapGrid(
2022 36
        Geographic2D|Geographic3D $to,
2023 36
        Angle $latitudeOfNaturalOrigin,
2024 36
        Angle $longitudeOfNaturalOrigin,
2025 36
        Length $falseEasting,
2026 36
        Length $falseNorthing
2027 36
    ): GeographicPoint {
2028 36
        $ellipsoid = $this->crs->getDatum()->getEllipsoid();
2029
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
2030 36
2031 36
        $z = new ComplexNumber(
2032 36
            $this->northing->subtract($falseNorthing)->divide($a)->asMetres()->getValue(),
2033 36
            $this->easting->subtract($falseEasting)->divide($a)->asMetres()->getValue(),
2034 36
        );
2035 36
2036
        $B1 = new ComplexNumber(0.7557853228, 0.0);
2037 36
        $B2 = new ComplexNumber(0.249204646, 0.003371507);
2038
        $B3 = new ComplexNumber(-0.001541739, 0.041058560);
2039
        $B4 = new ComplexNumber(-0.10162907, 0.01727609);
2040 36
        $B5 = new ComplexNumber(-0.26623489, -0.36249218);
2041 36
        $B6 = new ComplexNumber(-0.6870983, -1.1651967);
2042 36
        $b1 = new ComplexNumber(1.3231270439, 0.0);
2043 36
        $b2 = new ComplexNumber(-0.577245789, -0.007809598);
2044 36
        $b3 = new ComplexNumber(0.508307513, -0.112208952);
2045 36
        $b4 = new ComplexNumber(-0.15094762, 0.18200602);
2046 36
        $b5 = new ComplexNumber(1.01418179, 1.64497696);
2047 36
        $b6 = new ComplexNumber(1.9660549, 2.5127645);
2048 36
2049 36
        $zeta = new ComplexNumber(0, 0);
2050 36
        $zeta = $zeta->add($b1->multiply($z->pow(1)));
2051
        $zeta = $zeta->add($b2->multiply($z->pow(2)));
2052 36
        $zeta = $zeta->add($b3->multiply($z->pow(3)));
2053 36
        $zeta = $zeta->add($b4->multiply($z->pow(4)));
2054
        $zeta = $zeta->add($b5->multiply($z->pow(5)));
2055 36
        $zeta = $zeta->add($b6->multiply($z->pow(6)));
2056
2057
        for ($iterations = 0; $iterations < 2; ++$iterations) {
2058
            $numerator = $z;
2059
            $numerator = $numerator->add($B2->multiply($zeta->pow(2))->multiply(new ComplexNumber(1, 0)));
2060
            $numerator = $numerator->add($B3->multiply($zeta->pow(3))->multiply(new ComplexNumber(2, 0)));
2061
            $numerator = $numerator->add($B4->multiply($zeta->pow(4))->multiply(new ComplexNumber(3, 0)));
2062
            $numerator = $numerator->add($B5->multiply($zeta->pow(5))->multiply(new ComplexNumber(4, 0)));
2063 81
            $numerator = $numerator->add($B6->multiply($zeta->pow(6))->multiply(new ComplexNumber(5, 0)));
2064
2065
            $denominator = $B1;
2066
            $denominator = $denominator->add($B2->multiply($zeta->pow(1))->multiply(new ComplexNumber(2, 0)));
2067
            $denominator = $denominator->add($B3->multiply($zeta->pow(2))->multiply(new ComplexNumber(3, 0)));
2068
            $denominator = $denominator->add($B4->multiply($zeta->pow(3))->multiply(new ComplexNumber(4, 0)));
2069
            $denominator = $denominator->add($B5->multiply($zeta->pow(4))->multiply(new ComplexNumber(5, 0)));
2070
            $denominator = $denominator->add($B6->multiply($zeta->pow(5))->multiply(new ComplexNumber(6, 0)));
2071
2072
            $zeta = $numerator->divide($denominator);
2073
        }
2074
2075
        $deltaPsi = $zeta->getReal();
2076
        $deltaLatitudeToOrigin = 0;
2077
        $deltaLatitudeToOrigin += 1.5627014243 * $deltaPsi ** 1;
2078
        $deltaLatitudeToOrigin += 0.5185406398 * $deltaPsi ** 2;
2079
        $deltaLatitudeToOrigin += -0.03333098 * $deltaPsi ** 3;
2080 81
        $deltaLatitudeToOrigin += -0.1052906 * $deltaPsi ** 4;
2081 81
        $deltaLatitudeToOrigin += -0.0368594 * $deltaPsi ** 5;
2082 81
        $deltaLatitudeToOrigin += 0.007317 * $deltaPsi ** 6;
2083 81
        $deltaLatitudeToOrigin += 0.01220 * $deltaPsi ** 7;
2084 81
        $deltaLatitudeToOrigin += 0.00394 * $deltaPsi ** 8;
2085 81
        $deltaLatitudeToOrigin += -0.0013 * $deltaPsi ** 9;
2086
2087 81
        $latitude = $latitudeOfNaturalOrigin->add(new ArcSecond($deltaLatitudeToOrigin / 0.00001));
2088 81
        $longitude = $longitudeOfNaturalOrigin->add(new Radian($zeta->getImaginary()));
2089
2090 81
        return GeographicPoint::create($to, $latitude, $longitude, null, $this->epoch);
2091 81
    }
2092 81
2093 81
    /**
2094 81
     * Complex polynomial.
2095
     * Coordinate pairs treated as complex numbers.  This exploits the correlation between the polynomial coefficients
2096 81
     * and leads to a smaller number of coefficients than the general polynomials.
2097 81
     */
2098
    public function complexPolynomial(
2099 81
        Projected $to,
2100 81
        Length $ordinate1OfEvaluationPointInSourceCRS,
2101
        Length $ordinate2OfEvaluationPointInSourceCRS,
2102 81
        Length $ordinate1OfEvaluationPointInTargetCRS,
2103 81
        Length $ordinate2OfEvaluationPointInTargetCRS,
2104 81
        Scale $scalingFactorForSourceCRSCoordDifferences,
2105
        Scale $scalingFactorForTargetCRSCoordDifferences,
2106 81
        Scale $A1,
2107
        Scale $A2,
2108
        Scale $A3,
2109
        Scale $A4,
2110
        Scale $A5,
2111
        Scale $A6,
2112
        ?Scale $A7 = null,
2113
        ?Scale $A8 = null
2114
    ): self {
2115 7
        $xs = $this->easting->getValue();
2116
        $ys = $this->northing->getValue();
2117
        $xso = $ordinate1OfEvaluationPointInSourceCRS->getValue();
2118
        $yso = $ordinate2OfEvaluationPointInSourceCRS->getValue();
2119 7
        $xto = $ordinate1OfEvaluationPointInTargetCRS->getValue();
2120
        $yto = $ordinate2OfEvaluationPointInTargetCRS->getValue();
2121 7
2122
        $U = $scalingFactorForSourceCRSCoordDifferences->asUnity()->getValue() * ($xs - $xso);
2123
        $V = $scalingFactorForSourceCRSCoordDifferences->asUnity()->getValue() * ($ys - $yso);
2124
2125
        $mTdXdY = new ComplexNumber(0, 0);
2126
        $mTdXdY = $mTdXdY->add((new ComplexNumber($A1->getValue(), $A2->getValue()))->multiply(new ComplexNumber($U, $V))->pow(1));
2127
        $mTdXdY = $mTdXdY->add((new ComplexNumber($A3->getValue(), $A4->getValue()))->multiply((new ComplexNumber($U, $V))->pow(2)));
2128
        $mTdXdY = $mTdXdY->add((new ComplexNumber($A5->getValue(), $A6->getValue()))->multiply((new ComplexNumber($U, $V))->pow(3)));
2129
        $mTdXdY = $mTdXdY->add((new ComplexNumber($A7 ? $A7->getValue() : 0, $A8 ? $A8->getValue() : 0))->multiply((new ComplexNumber($U, $V))->pow(4)));
2130
2131
        $xt = $xs - $xso + $xto + $mTdXdY->getReal() / $scalingFactorForTargetCRSCoordDifferences->asUnity()->getValue();
2132
        $yt = $ys - $yso + $yto + $mTdXdY->getImaginary() / $scalingFactorForTargetCRSCoordDifferences->asUnity()->getValue();
2133
2134
        $xtUnit = $to->getCoordinateSystem()->getAxes()[0]->getUnitOfMeasureId();
2135
        $ytUnit = $to->getCoordinateSystem()->getAxes()[1]->getUnitOfMeasureId();
2136
2137
        return static::createFromEastingNorthing(
2138
            $to,
2139
            Length::makeUnit($xt, $xtUnit),
2140
            Length::makeUnit($yt, $ytUnit),
2141
            $this->epoch
2142
        );
2143
    }
2144
2145
    /**
2146
     * Ordnance Survey National Transformation
2147
     * Geodetic transformation between ETRS89 (or WGS 84) and OSGB36 / National Grid.  Uses ETRS89 / National Grid as
2148
     * an intermediate coordinate system for bi-linear interpolation of gridded grid coordinate differences.
2149
     */
2150
    public function OSTN15(
2151
        Geographic2D $to,
2152
        OSTNOSGM15Grid $eastingAndNorthingDifferenceFile
2153
    ): GeographicPoint {
2154
        $asETRS89 = $eastingAndNorthingDifferenceFile->applyReverseHorizontalAdjustment($this);
2155
2156
        return $asETRS89->transverseMercator($to, new Degree(49), new Degree(-2), new Unity(0.9996012717), new Metre(400000), new Metre(-100000));
2157
    }
2158
}
2159