Point::acos()   A
last analyzed

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 5
cp 0
crap 12
rs 10
1
<?php
2
3
/**
4
 * PHPCoord.
5
 *
6
 * @author Doug Wright
7
 */
8
declare(strict_types=1);
9
10
namespace PHPCoord\Point;
11
12
use DateTimeImmutable;
13
use PHPCoord\CoordinateOperation\CoordinateOperationMethods;
14
use PHPCoord\CoordinateOperation\CoordinateOperations;
0 ignored issues
show
Bug introduced by
The type PHPCoord\CoordinateOperation\CoordinateOperations 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...
15
use PHPCoord\CoordinateOperation\Grid;
16
use PHPCoord\CoordinateOperation\GridProvider;
17
use PHPCoord\CoordinateReferenceSystem\Compound;
18
use PHPCoord\CoordinateReferenceSystem\CoordinateReferenceSystem;
19
use PHPCoord\CoordinateReferenceSystem\Geocentric;
20
use PHPCoord\CoordinateReferenceSystem\Geographic2D;
21
use PHPCoord\CoordinateReferenceSystem\Geographic3D;
22
use PHPCoord\CoordinateReferenceSystem\Projected;
23
use PHPCoord\CoordinateReferenceSystem\Vertical;
24
use PHPCoord\UnitOfMeasure\Angle\Angle;
25
use PHPCoord\UnitOfMeasure\Length\Length;
26
use PHPCoord\UnitOfMeasure\Scale\Coefficient;
27
use PHPCoord\UnitOfMeasure\Scale\Scale;
28
use PHPCoord\UnitOfMeasure\UnitOfMeasure;
29
use Stringable;
30
31
use function acos;
32
use function asin;
33
use function assert;
34
use function atan;
35
use function atan2;
36
use function class_exists;
37
use function is_string;
38
use function property_exists;
39
use function sscanf;
40
use function str_ends_with;
41
use function str_starts_with;
42
43
abstract class Point implements Stringable
44
{
45
    protected const ITERATION_CONVERGENCE_FORMULA = 1e-10;
46
    protected const ITERATION_CONVERGENCE_GRID = 0.0001;
47
    protected const METHODS_REQUIRING_HORIZONTAL_POINT = [
48
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_AND_SLOPE => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_AND_SLOPE,
49
        CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019 => CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019,
50
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX,
51
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_PL_TXT => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_PL_TXT,
52
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_BEV_AT => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_BEV_AT,
53
        CoordinateOperationMethods::EPSG_VERTICAL_CHANGE_BY_GEOID_GRID_DIFFERENCE_NRCAN => CoordinateOperationMethods::EPSG_VERTICAL_CHANGE_BY_GEOID_GRID_DIFFERENCE_NRCAN,
54
    ];
55
    protected const METHODS_THAT_REQUIRE_DIRECTION = [
56
        CoordinateOperationMethods::EPSG_SIMILARITY_TRANSFORMATION => CoordinateOperationMethods::EPSG_SIMILARITY_TRANSFORMATION,
57
        CoordinateOperationMethods::EPSG_AFFINE_PARAMETRIC_TRANSFORMATION => CoordinateOperationMethods::EPSG_AFFINE_PARAMETRIC_TRANSFORMATION,
58
        CoordinateOperationMethods::EPSG_NADCON5_2D => CoordinateOperationMethods::EPSG_NADCON5_2D,
59
        CoordinateOperationMethods::EPSG_NADCON5_3D => CoordinateOperationMethods::EPSG_NADCON5_3D,
60
        CoordinateOperationMethods::EPSG_NTV2 => CoordinateOperationMethods::EPSG_NTV2,
61
        CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019 => CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019,
62
        CoordinateOperationMethods::EPSG_GEOCENTRIC_TRANSLATION_BY_GRID_INTERPOLATION_IGN => CoordinateOperationMethods::EPSG_GEOCENTRIC_TRANSLATION_BY_GRID_INTERPOLATION_IGN,
63
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX,
64
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_PL_TXT => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_PL_TXT,
65
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_BEV_AT => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_BEV_AT,
66
        CoordinateOperationMethods::EPSG_VERTICAL_CHANGE_BY_GEOID_GRID_DIFFERENCE_NRCAN => CoordinateOperationMethods::EPSG_VERTICAL_CHANGE_BY_GEOID_GRID_DIFFERENCE_NRCAN,
67
        CoordinateOperationMethods::EPSG_GENERAL_POLYNOMIAL_OF_DEGREE_2 => CoordinateOperationMethods::EPSG_GENERAL_POLYNOMIAL_OF_DEGREE_2,
68
        CoordinateOperationMethods::EPSG_GENERAL_POLYNOMIAL_OF_DEGREE_6 => CoordinateOperationMethods::EPSG_GENERAL_POLYNOMIAL_OF_DEGREE_6,
69
    ];
70
71
    /**
72
     * @var array<class-string, Grid>
0 ignored issues
show
Documentation Bug introduced by
The doc comment array<class-string, Grid> at position 2 could not be parsed: Unknown type name 'class-string' at position 2 in array<class-string, Grid>.
Loading history...
73
     */
74
    protected static array $gridCache = [];
75
76
    /**
77
     * @internal
78
     * @param array{horizontalPoint?: Point} $additionalParams
79 373
     */
80
    public function performOperation(string $srid, Compound|Geocentric|Geographic2D|Geographic3D|Projected|Vertical $to, bool $inReverse, array $additionalParams = []): self
81 373
    {
82
        $operation = CoordinateOperations::getOperationData($srid);
83 373
84
        if ($operation['method'] === CoordinateOperationMethods::EPSG_ALIAS) {
85
            $point = clone $this;
86
            assert(property_exists($point, 'crs'));
87
            $point->crs = $to;
0 ignored issues
show
Bug Best Practice introduced by
The property crs does not exist. Although not strictly required by PHP, it is generally a best practice to declare properties explicitly.
Loading history...
88
89
            return $point;
90 373
        } else {
91 373
            $method = CoordinateOperationMethods::getFunctionName($operation['method']);
92
            $params = self::resolveParamsByOperation($srid, $operation['method'], $inReverse);
93 373
94 9
            if (isset(self::METHODS_REQUIRING_HORIZONTAL_POINT[$operation['method']]) && isset($additionalParams['horizontalPoint'])) {
95
                $params['horizontalPoint'] = $additionalParams['horizontalPoint'];
96
            }
97 373
98
            return $this->$method($to, ...$params);
99
        }
100
    }
101
102
    /**
103
     * @return array<string, mixed>
104 373
     */
105
    protected static function resolveParamsByOperation(string $operationSrid, string $methodSrid, bool $inReverse): array
106 373
    {
107
        $methodData = CoordinateOperationMethods::getMethodData($methodSrid);
108 373
109 373
        $params = [];
110 373
        $powerCoefficients = [];
111 373
        foreach (CoordinateOperations::getParamData($operationSrid) as $paramName => $paramValue) {
112 6
            if (str_ends_with($paramName, 'File') && is_string($paramValue) && class_exists($paramValue) && new $paramValue() instanceof GridProvider) {
113
                $params[$paramName] = static::$gridCache[$paramValue] ??= (new $paramValue())->provideGrid();
114 373
            } else {
115 162
                if ($inReverse && isset($methodData['paramData'][$paramName]) && $methodData['paramData'][$paramName]['reverses']) {
116
                    $paramValue = $paramValue->multiply(-1);
117 373
                }
118 54
                if (str_starts_with($paramName, 'Au') || str_starts_with($paramName, 'Bu')) {
119
                    $powerCoefficients[$paramName] = $paramValue;
120 373
                } else {
121
                    $params[$paramName] = $paramValue;
122
                }
123
            }
124 373
        }
125 54
        if ($powerCoefficients) {
126
            $params['powerCoefficients'] = $powerCoefficients;
127 373
        }
128 39
        if (isset(self::METHODS_THAT_REQUIRE_DIRECTION[$methodSrid])) {
129
            if ($methodSrid === CoordinateOperationMethods::EPSG_VERTICAL_CHANGE_BY_GEOID_GRID_DIFFERENCE_NRCAN) {
130
                $inReverse = !$inReverse; // has a reversed grid
131 39
            }
132
            $params['inReverse'] = $inReverse;
133
        }
134 373
135
        return $params;
136
    }
137 100
138
    protected static function sign(float $number): int
139 100
    {
140
        if ($number < 0) {
141
            return -1;
142
        }
143 100
144
        return 1;
145
    }
146
147
    /**
148
     * General polynomial.
149
     * @param  array<string, Coefficient>  $powerCoefficients
150
     * @return array{xt: float, yt: float}
151 45
     */
152
    protected function generalPolynomialUnitless(
153
        float $xs,
154
        float $ys,
155
        UnitOfMeasure $ordinate1OfEvaluationPointInSourceCRS,
156
        UnitOfMeasure $ordinate2OfEvaluationPointInSourceCRS,
157
        UnitOfMeasure $ordinate1OfEvaluationPointInTargetCRS,
158
        UnitOfMeasure $ordinate2OfEvaluationPointInTargetCRS,
159
        Scale $scalingFactorForSourceCRSCoordDifferences,
160
        Scale $scalingFactorForTargetCRSCoordDifferences,
161
        Scale $A0,
162
        Scale $B0,
163
        array $powerCoefficients,
164
        bool $inReverse
165 45
    ): array {
166 18
        if (!$inReverse) {
167 18
            return $this->generalPolynomialUnitlessForward(
168 18
                $xs,
169 18
                $ys,
170 18
                $ordinate1OfEvaluationPointInSourceCRS,
171 18
                $ordinate2OfEvaluationPointInSourceCRS,
172 18
                $ordinate1OfEvaluationPointInTargetCRS,
173 18
                $ordinate2OfEvaluationPointInTargetCRS,
174 18
                $scalingFactorForSourceCRSCoordDifferences,
175 18
                $scalingFactorForTargetCRSCoordDifferences,
176 18
                $A0,
177 18
                $B0,
178 18
                $powerCoefficients,
179
            );
180 27
        } else {
181 27
            return $this->generalPolynomialUnitlessReverse(
182 27
                $xs,
183 27
                $ys,
184 27
                $ordinate1OfEvaluationPointInSourceCRS,
185 27
                $ordinate2OfEvaluationPointInSourceCRS,
186 27
                $ordinate1OfEvaluationPointInTargetCRS,
187 27
                $ordinate2OfEvaluationPointInTargetCRS,
188 27
                $scalingFactorForSourceCRSCoordDifferences,
189 27
                $scalingFactorForTargetCRSCoordDifferences,
190 27
                $A0,
191 27
                $B0,
192 27
                $powerCoefficients,
193
            );
194
        }
195
    }
196 45
197
    protected function generalPolynomialUnitlessForward(
198
        float $xs,
199
        float $ys,
200
        UnitOfMeasure $ordinate1OfEvaluationPointInSourceCRS,
201
        UnitOfMeasure $ordinate2OfEvaluationPointInSourceCRS,
202
        UnitOfMeasure $ordinate1OfEvaluationPointInTargetCRS,
203
        UnitOfMeasure $ordinate2OfEvaluationPointInTargetCRS,
204
        Scale $scalingFactorForSourceCRSCoordDifferences,
205
        Scale $scalingFactorForTargetCRSCoordDifferences,
206
        Scale $A0,
207
        Scale $B0,
208
        array $powerCoefficients
209 45
    ): array {
210 45
        $xso = $ordinate1OfEvaluationPointInSourceCRS->getValue();
211 45
        $yso = $ordinate2OfEvaluationPointInSourceCRS->getValue();
212 45
        $xto = $ordinate1OfEvaluationPointInTargetCRS->getValue();
213
        $yto = $ordinate2OfEvaluationPointInTargetCRS->getValue();
214 45
215 45
        $U = $scalingFactorForSourceCRSCoordDifferences->asUnity()->getValue() * ($xs - $xso);
216
        $V = $scalingFactorForSourceCRSCoordDifferences->asUnity()->getValue() * ($ys - $yso);
217 45
218 45
        $mTdX = $A0->getValue();
219 45
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
220 45
            if ($coefficientName[0] === 'A') {
221 45
                sscanf($coefficientName, 'Au%dv%d', $uPower, $vPower);
222
                $mTdX += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
223
            }
224
        }
225 45
226 45
        $mTdY = $B0->getValue();
227 45
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
228 45
            if ($coefficientName[0] === 'B') {
229 45
                sscanf($coefficientName, 'Bu%dv%d', $uPower, $vPower);
230
                $mTdY += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
231
            }
232
        }
233 45
234 45
        $xt = $xs - $xso + $xto + $mTdX / $scalingFactorForTargetCRSCoordDifferences->asUnity()->getValue();
235
        $yt = $ys - $yso + $yto + $mTdY / $scalingFactorForTargetCRSCoordDifferences->asUnity()->getValue();
236 45
237
        return ['xt' => $xt, 'yt' => $yt];
238
    }
239 27
240
    protected function generalPolynomialUnitlessReverse(
241
        float $xs,
242
        float $ys,
243
        UnitOfMeasure $ordinate1OfEvaluationPointInSourceCRS,
244
        UnitOfMeasure $ordinate2OfEvaluationPointInSourceCRS,
245
        UnitOfMeasure $ordinate1OfEvaluationPointInTargetCRS,
246
        UnitOfMeasure $ordinate2OfEvaluationPointInTargetCRS,
247
        Scale $scalingFactorForSourceCRSCoordDifferences,
248
        Scale $scalingFactorForTargetCRSCoordDifferences,
249
        Scale $A0,
250
        Scale $B0,
251
        array $powerCoefficients
252 27
    ): array {
253 27
        $result = ['xt' => $xs, 'yt' => $ys];
254 27
        for ($i = 0; $i <= 15; ++$i) {
255 27
            $forwardShiftedCoordinates = $this->generalPolynomialUnitlessForward(
256 27
                $result['xt'],
257 27
                $result['yt'],
258 27
                $ordinate1OfEvaluationPointInSourceCRS,
259 27
                $ordinate2OfEvaluationPointInSourceCRS,
260 27
                $ordinate1OfEvaluationPointInTargetCRS,
261 27
                $ordinate2OfEvaluationPointInTargetCRS,
262 27
                $scalingFactorForSourceCRSCoordDifferences,
263 27
                $scalingFactorForTargetCRSCoordDifferences,
264 27
                $A0,
265 27
                $B0,
266 27
                $powerCoefficients
267 27
            );
268 27
            $deltaError = [
269 27
                'xt' => $forwardShiftedCoordinates['xt'] - $xs,
270 27
                'yt' => $forwardShiftedCoordinates['yt'] - $ys,
271 27
            ];
272 27
            $result['xt'] -= $deltaError['xt'];
273
            $result['yt'] -= $deltaError['yt'];
274
        }
275 27
276
        return $result;
277
    }
278
279
    /**
280
     * Reversible polynomial.
281
     * @param  array<string, Coefficient>  $powerCoefficients
282
     * @return array{xt: float, yt: float}
283 36
     */
284
    protected function reversiblePolynomialUnitless(
285
        float $xs,
286
        float $ys,
287
        Angle $ordinate1OfEvaluationPoint,
288
        Angle $ordinate2OfEvaluationPoint,
289
        Scale $scalingFactorForCoordDifferences,
290
        Scale $A0,
291
        Scale $B0,
292
        array $powerCoefficients
293 36
    ): array {
294 36
        $xo = $ordinate1OfEvaluationPoint->getValue();
295
        $yo = $ordinate2OfEvaluationPoint->getValue();
296 36
297 36
        $U = $scalingFactorForCoordDifferences->asUnity()->getValue() * ($xs - $xo);
298
        $V = $scalingFactorForCoordDifferences->asUnity()->getValue() * ($ys - $yo);
299 36
300 36
        $mTdX = $A0->getValue();
301 36
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
302 36
            if ($coefficientName[0] === 'A') {
303 36
                sscanf($coefficientName, 'Au%dv%d', $uPower, $vPower);
304
                $mTdX += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
305
            }
306
        }
307 36
308 36
        $mTdY = $B0->getValue();
309 36
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
310 36
            if ($coefficientName[0] === 'B') {
311 36
                sscanf($coefficientName, 'Bu%dv%d', $uPower, $vPower);
312
                $mTdY += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
313
            }
314
        }
315 36
316 36
        $xt = $xs + $mTdX * $scalingFactorForCoordDifferences->asUnity()->getValue();
317
        $yt = $ys + $mTdY * $scalingFactorForCoordDifferences->asUnity()->getValue();
318 36
319
        return ['xt' => $xt, 'yt' => $yt];
320
    }
321
322
    /**
323
     * Floating point vagaries mean that it's possible for inputs to be e.g. 1.00000000000001 which makes PHP give a
324
     * silent NaN as output so inputs need to be capped. atan/atan2 are not affected, they seem to cap internally.
325
     */
326
    protected static function acos(float $num): float
327
    {
328
        if ($num > 1.0) {
329
            $num = 1.0;
330
        } elseif ($num < -1) {
331
            $num = -1.0;
332
        }
333
334
        return acos($num);
335
    }
336
337
    /**
338
     * Floating point vagaries mean that it's possible for inputs to be e.g. 1.00000000000001 which makes PHP give a
339
     * silent NaN as output so inputs need to be capped. atan/atan2 are not affected, they seem to cap internally.
340 540
     */
341
    protected static function asin(float $num): float
342 540
    {
343
        if ($num > 1.0) {
344 540
            $num = 1.0;
345
        } elseif ($num < -1.0) {
346
            $num = -1.0;
347
        }
348 540
349
        return asin($num);
350
    }
351
352
    abstract public function getCRS(): CoordinateReferenceSystem;
353
354
    abstract public function getCoordinateEpoch(): ?DateTimeImmutable;
355
356
    abstract public function calculateDistance(self $to): Length;
357
}
358