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

Point::generalPolynomialUnitless()   A

Complexity

Conditions 5
Paths 9

Size

Total Lines 41
Code Lines 19

Duplication

Lines 0
Ratio 0 %

Code Coverage

Tests 20
CRAP Score 5

Importance

Changes 0
Metric Value
eloc 19
c 0
b 0
f 0
dl 0
loc 41
ccs 20
cts 20
cp 1
rs 9.3222
cc 5
nc 9
nop 11
crap 5

How to fix   Many Parameters   

Many Parameters

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

There are several approaches to avoid long parameter lists:

1
<?php
2
/**
3
 * PHPCoord.
4
 *
5
 * @author Doug Wright
6
 */
7
declare(strict_types=1);
8
9
namespace PHPCoord;
10
11
use function abs;
12
use function 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