Passed
Push — master ( 6a435f...dad6c0 )
by Doug
13:08
created

Point::acos()   A

Complexity

Conditions 3
Paths 3

Size

Total Lines 9
Code Lines 5

Duplication

Lines 0
Ratio 0 %

Code Coverage

Tests 0
CRAP Score 12

Importance

Changes 0
Metric Value
cc 3
eloc 5
c 0
b 0
f 0
nc 3
nop 1
dl 0
loc 9
ccs 0
cts 6
cp 0
crap 12
rs 10
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 acos;
13
use function array_reverse;
14
use function array_values;
15
use function asin;
16
use function atan;
17
use function atan2;
18
use function cos;
19
use DateTimeImmutable;
20
use function lcfirst;
21
use const M_PI;
22
use const PHP_MAJOR_VERSION;
23
use PHPCoord\CoordinateOperation\ConvertiblePoint;
24
use PHPCoord\CoordinateOperation\CoordinateOperationMethods;
25
use PHPCoord\CoordinateOperation\CoordinateOperationParams;
0 ignored issues
show
Bug introduced by
The type PHPCoord\CoordinateOpera...ordinateOperationParams 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...
26
use PHPCoord\CoordinateOperation\CoordinateOperations;
27
use PHPCoord\CoordinateOperation\GeographicValue;
28
use PHPCoord\CoordinateReferenceSystem\CoordinateReferenceSystem;
29
use PHPCoord\CoordinateSystem\Axis;
30
use PHPCoord\Datum\Ellipsoid;
31
use PHPCoord\UnitOfMeasure\Angle\Angle;
32
use PHPCoord\UnitOfMeasure\Length\Length;
33
use PHPCoord\UnitOfMeasure\Length\Metre;
34
use PHPCoord\UnitOfMeasure\Scale\Coefficient;
35
use PHPCoord\UnitOfMeasure\Scale\Scale;
36
use PHPCoord\UnitOfMeasure\UnitOfMeasure;
37
use PHPCoord\UnitOfMeasure\UnitOfMeasureFactory;
38
use function preg_match;
39
use function sin;
40
use function sqrt;
41
use function sscanf;
42
use function str_replace;
43
use function str_starts_with;
44
use Stringable;
45
use function tan;
46
use function ucwords;
47
48
abstract class Point implements Stringable
49
{
50
    protected const ITERATION_CONVERGENCE_FORMULA = 1e-10;
51
    protected const ITERATION_CONVERGENCE_GRID = 0.0001;
52
    protected const METHODS_REQUIRING_HORIZONTAL_POINT = [
53
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_AND_SLOPE => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_AND_SLOPE,
54
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX,
55
        CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019 => CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019,
56
    ];
57
    protected const METHODS_THAT_REQUIRE_DIRECTION = [
58
        CoordinateOperationMethods::EPSG_SIMILARITY_TRANSFORMATION => CoordinateOperationMethods::EPSG_SIMILARITY_TRANSFORMATION,
59
        CoordinateOperationMethods::EPSG_AFFINE_PARAMETRIC_TRANSFORMATION => CoordinateOperationMethods::EPSG_AFFINE_PARAMETRIC_TRANSFORMATION,
60
        CoordinateOperationMethods::EPSG_NADCON5_2D => CoordinateOperationMethods::EPSG_NADCON5_2D,
61
        CoordinateOperationMethods::EPSG_NADCON5_3D => CoordinateOperationMethods::EPSG_NADCON5_3D,
62
        CoordinateOperationMethods::EPSG_NTV2 => CoordinateOperationMethods::EPSG_NTV2,
63
        CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019 => CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019,
64
        CoordinateOperationMethods::EPSG_GEOCENTRIC_TRANSLATION_BY_GRID_INTERPOLATION_IGN => CoordinateOperationMethods::EPSG_GEOCENTRIC_TRANSLATION_BY_GRID_INTERPOLATION_IGN,
65
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX,
66
    ];
67
68
    /**
69
     * @internal
70
     */
71 245
    public function performOperation(string $srid, CoordinateReferenceSystem $to, bool $inReverse, array $additionalParams = []): self
72
    {
73 245
        $operations = self::resolveConcatenatedOperations($srid, $inReverse);
74
75 245
        $point = $this;
76 245
        foreach ($operations as $operationSrid => $operation) {
77 245
            $method = CoordinateOperationMethods::getFunctionName($operation['method']);
78 245
            if (isset($operation['source_crs'])) {
79
                $sourceCRS = $inReverse ? $operation['target_crs'] : $operation['source_crs'];
80
                $destCRS = CoordinateReferenceSystem::fromSRID($inReverse ? $operation['source_crs'] : $operation['target_crs']);
81
82
                // occasionally the CRS domain switches part way through the sequence (e.g. NAD86 geog to NAD86 geocen)
83
                if ($point instanceof ConvertiblePoint && $point->getCRS()->getSRID() !== $sourceCRS) {
84
                    $point = $point->convert(CoordinateReferenceSystem::fromSRID($sourceCRS), true);
85
                }
86
            } else {
87 245
                $destCRS = $to;
88
            }
89
90 245
            $params = self::resolveParamsByOperation($operationSrid, $operation['method'], $inReverse);
91
92 245
            if (isset(self::METHODS_REQUIRING_HORIZONTAL_POINT[$operation['method']])) {
93 9
                $params['horizontalPoint'] = $additionalParams['horizontalPoint'];
94
            }
95
96 245
            if (PHP_MAJOR_VERSION >= 8) {
97 245
                $point = $point->$method($destCRS, ...$params);
98
            } else {
99
                $point = $point->$method($destCRS, ...array_values($params));
100
            }
101
        }
102
103
        // some operations are reused across CRSses (e.g. ETRS89 and WGS84), so the $destCRS of the final suboperation might not be the intended target
104 245
        if ($point instanceof ConvertiblePoint && !$point->getCRS() instanceof $to) {// occasionally the CRS domain switches part way through the sequence (e.g. NAD86 geog to NAD86 geocen)
0 ignored issues
show
Bug introduced by
The method getCRS() does not exist on PHPCoord\CoordinateOperation\ConvertiblePoint. Since it exists in all sub-types, consider adding an abstract or default implementation to PHPCoord\CoordinateOperation\ConvertiblePoint. ( Ignorable by Annotation )

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

104
        if ($point instanceof ConvertiblePoint && !$point->/** @scrutinizer ignore-call */ getCRS() instanceof $to) {// occasionally the CRS domain switches part way through the sequence (e.g. NAD86 geog to NAD86 geocen)
Loading history...
105
            $point = $point->convert($to, true);
106
        }
107 245
        $point->crs = $to;
108
109 245
        return $point;
110
    }
111
112
    /**
113
     * @deprecated
114
     */
115
    protected static function camelCase(string $string): string
116
    {
117
        $string = str_replace([' ', '-', '(', ')'], '', ucwords($string, ' -()'));
118
        if (!preg_match('/^(EPSG|[ABC][uv\d])/', $string)) {
119
            $string = lcfirst($string);
120
        }
121
122
        return $string;
123
    }
124
125 254
    protected static function resolveConcatenatedOperations(string $operationSrid, bool $inReverse): array
126
    {
127 254
        $operations = [];
128 254
        $operation = CoordinateOperations::getOperationData($operationSrid);
129 254
        if (isset($operation['operations'])) {
130 16
            foreach ($operation['operations'] as $subOperation) {
131 16
                $subOperationData = CoordinateOperations::getOperationData($subOperation['operation']);
132 16
                $subOperationData['source_crs'] = $subOperation['source_crs'];
133 16
                $subOperationData['target_crs'] = $subOperation['target_crs'];
134 16
                $operations[$subOperation['operation']] = $subOperationData;
135
            }
136
        } else {
137 254
            $operations[$operationSrid] = $operation;
138
        }
139
140 254
        if ($inReverse) {
141 115
            $operations = array_reverse($operations, true);
142
        }
143
144 254
        return $operations;
145
    }
146
147 245
    protected static function resolveParamsByOperation(string $operationSrid, string $methodSrid, bool $inReverse): array
148
    {
149 245
        $params = [];
150 245
        $powerCoefficients = [];
151 245
        foreach (CoordinateOperationParams::getParamData($operationSrid) as $paramName => $paramData) {
152 245
            if (isset($paramData['fileProvider'])) {
153 11
                $params[$paramName] = (new $paramData['fileProvider']())->provideGrid();
154
            } else {
155 245
                if ($inReverse && $paramData['reverses']) {
156 50
                    $paramData['value'] *= -1;
157
                }
158 245
                if ($paramData['uom']) {
159 245
                    $param = UnitOfMeasureFactory::makeUnit($paramData['value'], $paramData['uom']);
160
                } else {
161 10
                    $param = $paramData['value'];
162
                }
163 245
                if (str_starts_with($paramName, 'Au') || str_starts_with($paramName, 'Bu')) {
164 27
                    $powerCoefficients[$paramName] = $param;
165
                } else {
166 245
                    $params[$paramName] = $param;
167
                }
168
            }
169
        }
170 245
        if ($powerCoefficients) {
171 27
            $params['powerCoefficients'] = $powerCoefficients;
172
        }
173 245
        if (isset(self::METHODS_THAT_REQUIRE_DIRECTION[$methodSrid])) {
174 4
            $params['inReverse'] = $inReverse;
175
        }
176
177 245
        return $params;
178
    }
179
180
    /**
181
     * @deprecated
182
     */
183
    protected function getAxisByName(string $name): ?Axis
184
    {
185
        foreach ($this->getCRS()->getCoordinateSystem()->getAxes() as $axis) {
186
            if ($axis->getName() === $name) {
187
                return $axis;
188
            }
189
        }
190
191
        return null;
192
    }
193
194 90
    protected static function sign(float $number): int
195
    {
196 90
        if ($number < 0) {
197
            return -1;
198
        }
199
200 90
        return 1;
201
    }
202
203
    /**
204
     * Calculate surface distance between two points.
205
     */
206 153
    protected static function vincenty(GeographicValue $from, GeographicValue $to, Ellipsoid $ellipsoid): Length
207
    {
208 153
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
209 153
        $b = $ellipsoid->getSemiMinorAxis()->asMetres()->getValue();
210 153
        $f = $ellipsoid->getInverseFlattening();
211 153
        $U1 = atan((1 - $f) * tan($from->getLatitude()->asRadians()->getValue()));
212 153
        $U2 = atan((1 - $f) * tan($to->getLatitude()->asRadians()->getValue()));
213 153
        $L = $to->getLongitude()->subtract($from->getLongitude())->asRadians()->getValue();
214
215 153
        $lambda = $L;
216
        do {
217 153
            $lambdaN = $lambda;
218
219 153
            $sinSigma = sqrt((cos($U2) * sin($lambda)) ** 2 + (cos($U1) * sin($U2) - sin($U1) * cos($U2) * cos($lambda)) ** 2);
220 153
            $cosSigma = sin($U1) * sin($U2) + cos($U1) * cos($U2) * cos($lambda);
221 153
            $sigma = atan2($sinSigma, $cosSigma);
222
223 153
            $sinAlpha = cos($U1) * cos($U2) * sin($lambda) / $sinSigma;
224 153
            $cosSqAlpha = (1 - $sinAlpha ** 2);
225 153
            $cos2SigmaM = $cosSqAlpha ? $cosSigma - (2 * sin($U1) * sin($U2) / $cosSqAlpha) : 0;
226 153
            $C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
227 153
            $lambda = $L + (1 - $C) * $f * $sinAlpha * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (-1 + 2 * $cos2SigmaM ** 2)));
228 153
        } while (abs($lambda - $lambdaN) >= static::ITERATION_CONVERGENCE_FORMULA && abs($lambda) < M_PI);
229
230
        // Antipodal case
231 153
        if (abs($lambda) >= M_PI) {
232 27
            if ($L >= 0) {
233 27
                $LPrime = M_PI - $L;
234
            } else {
235
                $LPrime = -M_PI - $L;
236
            }
237
238 27
            $lambdaPrime = 0;
239 27
            $sigma = M_PI - abs($U1 + $U2);
240 27
            $sinSigma = sin($sigma);
241 27
            $cosSqAlpha = 0.5;
242 27
            $sinAlpha = 0;
243
244
            do {
245 27
                $sinAlphaN = $sinAlpha;
246
247 27
                $C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
248 27
                $cos2SigmaM = cos($sigma) - 2 * sin($U1) * sin($U2) / $cosSqAlpha;
249 27
                $D = (1 - $C) * $f * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * cos($sigma) * (-1 + 2 * $cos2SigmaM ** 2)));
250 27
                $sinAlpha = ($LPrime - $lambdaPrime) / $D;
251 27
                $cosSqAlpha = (1 - $sinAlpha ** 2);
252 27
                $sinLambdaPrime = ($sinAlpha * $sinSigma) / (cos($U1) * cos($U2));
253 27
                $lambdaPrime = self::asin($sinLambdaPrime);
254 27
                $sinSqSigma = (cos($U2) * $sinLambdaPrime) ** 2 + (cos($U1) * sin($U2) + sin($U1) * cos($U2) * cos($lambdaPrime)) ** 2;
255 27
                $sinSigma = sqrt($sinSqSigma);
256 27
            } while (abs($sinAlpha - $sinAlphaN) >= static::ITERATION_CONVERGENCE_FORMULA);
257
        }
258
259 153
        $E = sqrt(1 + (($a ** 2 - $b ** 2) / $b ** 2) * $cosSqAlpha);
260 153
        $F = ($E - 1) / ($E + 1);
261
262 153
        $A = (1 + $F ** 2 / 4) / (1 - $F);
263 153
        $B = $F * (1 - 3 / 8 * $F ** 2);
264
265 153
        $deltaSigma = $B * $sinSigma * ($cos2SigmaM + $B / 4 * ($cosSigma * (-1 + 2 * $cos2SigmaM ** 2) - $B / 6 * $cos2SigmaM * (-3 + 4 * $sinSigma ** 2) * (-3 + 4 * $cos2SigmaM ** 2)));
266
267 153
        return new Metre($b * $A * ($sigma - $deltaSigma));
268
    }
269
270
    /**
271
     * General polynomial.
272
     * @param Coefficient[] $powerCoefficients
273
     */
274 18
    protected function generalPolynomialUnitless(
275
        float $xs,
276
        float $ys,
277
        UnitOfMeasure $ordinate1OfEvaluationPointInSourceCRS,
278
        UnitOfMeasure $ordinate2OfEvaluationPointInSourceCRS,
279
        UnitOfMeasure $ordinate1OfEvaluationPointInTargetCRS,
280
        UnitOfMeasure $ordinate2OfEvaluationPointInTargetCRS,
281
        Scale $scalingFactorForSourceCRSCoordDifferences,
282
        Scale $scalingFactorForTargetCRSCoordDifferences,
283
        Scale $A0,
284
        Scale $B0,
285
        array $powerCoefficients
286
    ): array {
287 18
        $xso = $ordinate1OfEvaluationPointInSourceCRS->getValue();
288 18
        $yso = $ordinate2OfEvaluationPointInSourceCRS->getValue();
289 18
        $xto = $ordinate1OfEvaluationPointInTargetCRS->getValue();
290 18
        $yto = $ordinate2OfEvaluationPointInTargetCRS->getValue();
291
292 18
        $U = $scalingFactorForSourceCRSCoordDifferences->asUnity()->getValue() * ($xs - $xso);
293 18
        $V = $scalingFactorForSourceCRSCoordDifferences->asUnity()->getValue() * ($ys - $yso);
294
295 18
        $mTdX = $A0->getValue();
296 18
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
297 18
            if ($coefficientName[0] === 'A') {
298 18
                sscanf($coefficientName, 'Au%dv%d', $uPower, $vPower);
1 ignored issue
show
Comprehensibility Best Practice introduced by
The variable $vPower seems to be never defined.
Loading history...
299 18
                $mTdX += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
300
            }
301
        }
302
303 18
        $mTdY = $B0->getValue();
304 18
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
305 18
            if ($coefficientName[0] === 'B') {
306 18
                sscanf($coefficientName, 'Bu%dv%d', $uPower, $vPower);
307 18
                $mTdY += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
308
            }
309
        }
310
311 18
        $xt = $xs - $xso + $xto + $mTdX / $scalingFactorForTargetCRSCoordDifferences->asUnity()->getValue();
312 18
        $yt = $ys - $yso + $yto + $mTdY / $scalingFactorForTargetCRSCoordDifferences->asUnity()->getValue();
313
314 18
        return ['xt' => $xt, 'yt' => $yt];
315
    }
316
317
    /**
318
     * Reversible polynomial.
319
     */
320 36
    protected function reversiblePolynomialUnitless(
321
        float $xs,
322
        float $ys,
323
        Angle $ordinate1OfEvaluationPoint,
324
        Angle $ordinate2OfEvaluationPoint,
325
        Scale $scalingFactorForCoordDifferences,
326
        Scale $A0,
327
        Scale $B0,
328
        array $powerCoefficients
329
    ): array {
330 36
        $xo = $ordinate1OfEvaluationPoint->getValue();
331 36
        $yo = $ordinate2OfEvaluationPoint->getValue();
332
333 36
        $U = $scalingFactorForCoordDifferences->asUnity()->getValue() * ($xs - $xo);
334 36
        $V = $scalingFactorForCoordDifferences->asUnity()->getValue() * ($ys - $yo);
335
336 36
        $mTdX = $A0->getValue();
337 36
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
338 36
            if ($coefficientName[0] === 'A') {
339 36
                sscanf($coefficientName, 'Au%dv%d', $uPower, $vPower);
1 ignored issue
show
Comprehensibility Best Practice introduced by
The variable $vPower seems to be never defined.
Loading history...
340 36
                $mTdX += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
341
            }
342
        }
343
344 36
        $mTdY = $B0->getValue();
345 36
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
346 36
            if ($coefficientName[0] === 'B') {
347 36
                sscanf($coefficientName, 'Bu%dv%d', $uPower, $vPower);
348 36
                $mTdY += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
349
            }
350
        }
351
352 36
        $xt = $xs + $mTdX * $scalingFactorForCoordDifferences->asUnity()->getValue();
353 36
        $yt = $ys + $mTdY * $scalingFactorForCoordDifferences->asUnity()->getValue();
354
355 36
        return ['xt' => $xt, 'yt' => $yt];
356
    }
357
358
    /**
359
     * Floating point vagaries mean that it's possible for inputs to be e.g. 1.00000000000001 which makes PHP give a
360
     * silent NaN as output so inputs need to be capped. atan/atan2 are not affected, they seem to cap internally.
361
     */
362
    protected static function acos(float $num): float
363
    {
364
        if ($num > 1.0) {
365
            $num = 1.0;
366
        } elseif ($num < -1) {
367
            $num = -1.0;
368
        }
369
370
        return acos($num);
371
    }
372
373
    /**
374
     * Floating point vagaries mean that it's possible for inputs to be e.g. 1.00000000000001 which makes PHP give a
375
     * silent NaN as output so inputs need to be capped. atan/atan2 are not affected, they seem to cap internally.
376
     */
377 500
    protected static function asin(float $num): float
378
    {
379 500
        if ($num > 1.0) {
380
            $num = 1.0;
381 500
        } elseif ($num < -1.0) {
382
            $num = -1.0;
383
        }
384
385 500
        return asin($num);
386
    }
387
388
    abstract public function getCRS(): CoordinateReferenceSystem;
389
390
    abstract public function getCoordinateEpoch(): ?DateTimeImmutable;
391
392
    abstract public function calculateDistance(self $to): Length;
393
}
394