Passed
Push — master ( aca036...b09582 )
by Doug
83:35 queued 68:07
created

Point::reversiblePolynomialUnitless()   A

Complexity

Conditions 5
Paths 9

Size

Total Lines 36
Code Lines 17

Duplication

Lines 0
Ratio 0 %

Code Coverage

Tests 18
CRAP Score 5

Importance

Changes 0
Metric Value
eloc 17
c 0
b 0
f 0
dl 0
loc 36
ccs 18
cts 18
cp 1
rs 9.3888
cc 5
nc 9
nop 8
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 in_array;
21
use function lcfirst;
22
use const M_PI;
23
use const PHP_MAJOR_VERSION;
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\CoordinateOperation\GridProvider;
29
use PHPCoord\CoordinateReferenceSystem\CoordinateReferenceSystem;
30
use PHPCoord\CoordinateSystem\Axis;
31
use PHPCoord\Datum\Ellipsoid;
32
use PHPCoord\UnitOfMeasure\Angle\Angle;
33
use PHPCoord\UnitOfMeasure\Length\Length;
34
use PHPCoord\UnitOfMeasure\Length\Metre;
35
use PHPCoord\UnitOfMeasure\Scale\Coefficient;
36
use PHPCoord\UnitOfMeasure\Scale\Scale;
37
use PHPCoord\UnitOfMeasure\UnitOfMeasure;
38
use PHPCoord\UnitOfMeasure\UnitOfMeasureFactory;
39
use function preg_match;
40
use function sin;
41
use function sqrt;
42
use function sscanf;
43
use function str_replace;
44
use Stringable;
45
use function tan;
46
use function ucwords;
47
48
abstract class Point implements Stringable
49
{
50
    protected const ITERATION_CONVERGENCE = 1e-12;
51
52
    /**
53
     * @internal
54
     */
55 225
    public function performOperation(string $srid, CoordinateReferenceSystem $to, bool $inReverse, array $additionalParams = []): self
56
    {
57 225
        $operations = self::resolveConcatenatedOperations($srid, $inReverse);
58
59 225
        $point = $this;
60 225
        foreach ($operations as $operationSrid => $operation) {
61 225
            $method = CoordinateOperationMethods::getFunctionName($operation['method']);
62 225
            if (isset($operation['source_crs'])) {
63
                $destCRS = CoordinateReferenceSystem::fromSRID($inReverse ? $operation['source_crs'] : $operation['target_crs']);
64
            } else {
65 225
                $destCRS = $to;
66
            }
67
68 225
            $params = self::resolveParamsByOperation($operationSrid, $operation['method'], $inReverse);
69
70 225
            if ($operation['method'] === CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_AND_SLOPE) {
71 9
                $params['horizontalPoint'] = $additionalParams['horizontalPoint'];
72
            }
73
74 225
            if (PHP_MAJOR_VERSION >= 8) {
75
                $point = $point->$method($destCRS, ...$params);
76
            } else {
77 225
                $point = $point->$method($destCRS, ...array_values($params));
78
            }
79
        }
80
81 225
        $point->crs = $to; //some operations are reused across CRSses (e.g. ETRS89 and WGS84), so the $destCRS of the final suboperation might not be the intended target
82
83 225
        return $point;
84
    }
85
86 225
    protected static function camelCase(string $string): string
87
    {
88 225
        $string = str_replace([' ', '-'], '', ucwords($string, ' -'));
89 225
        if (!preg_match('/^(EPSG|[ABC][uv\d])/', $string)) {
90 225
            $string = lcfirst($string);
91
        }
92
93 225
        return $string;
94
    }
95
96 225
    protected static function resolveConcatenatedOperations(string $operationSrid, bool $inReverse): array
97
    {
98 225
        $operations = [];
99 225
        $operation = CoordinateOperations::getOperationData($operationSrid);
100 225
        if (isset($operation['operations'])) {
101 8
            foreach ($operation['operations'] as $subOperation) {
102 8
                $subOperationData = CoordinateOperations::getOperationData($subOperation['operation']);
103 8
                $subOperationData['source_crs'] = $subOperation['source_crs'];
104 8
                $subOperationData['target_crs'] = $subOperation['target_crs'];
105 8
                $operations[$subOperation['operation']] = $subOperationData;
106
            }
107
        } else {
108 225
            $operations[$operationSrid] = $operation;
109
        }
110
111 225
        if ($inReverse) {
112 108
            $operations = array_reverse($operations, true);
113
        }
114
115 225
        return $operations;
116
    }
117
118 225
    protected static function resolveParamsByOperation(string $operationSrid, string $methodSrid, bool $inReverse): array
119
    {
120 225
        $params = [];
121 225
        $powerCoefficients = [];
122 225
        foreach (CoordinateOperationParams::getParamData($operationSrid) as $paramName => $paramData) {
123 225
            if (isset($paramData['fileProvider'])) {
124
                $paramName = static::camelCase($paramName);
125
                /** @var GridProvider $provider */
126
                $provider = new $paramData['fileProvider']();
127
                $params[$paramName] = $provider->provideGrid();
128
            } else {
129 225
                $value = $paramData['value'];
130 225
                if ($inReverse && $paramData['reverses']) {
131 45
                    $value *= -1;
132
                }
133 225
                if ($paramData['uom']) {
134 225
                    $param = UnitOfMeasureFactory::makeUnit($value, $paramData['uom']);
135
                } else {
136 9
                    $param = $paramData['value'];
137
                }
138 225
                $paramName = static::camelCase($paramName);
139 225
                if (preg_match('/^(Au|Bu)/', $paramName)) {
140 27
                    $powerCoefficients[$paramName] = $param;
141
                } else {
142 225
                    $params[$paramName] = $param;
143
                }
144
            }
145
        }
146 225
        if ($powerCoefficients) {
147 27
            $params['powerCoefficients'] = $powerCoefficients;
148
        }
149 225
        if (in_array($methodSrid, [
150 225
            CoordinateOperationMethods::EPSG_SIMILARITY_TRANSFORMATION,
151 225
            CoordinateOperationMethods::EPSG_AFFINE_PARAMETRIC_TRANSFORMATION,
152 225
        ], true)) {
153
            $params['inReverse'] = $inReverse;
154
        }
155
156 225
        return $params;
157
    }
158
159 2043
    protected function getAxisByName(string $name): ?Axis
160
    {
161 2043
        foreach ($this->getCRS()->getCoordinateSystem()->getAxes() as $axis) {
162 2043
            if ($axis->getName() === $name) {
163 2043
                return $axis;
164
            }
165
        }
166
167 1314
        return null;
168
    }
169
170 90
    protected static function sign(float $number): int
171
    {
172 90
        if ($number < 0) {
173
            return -1;
174
        }
175
176 90
        return 1;
177
    }
178
179
    /**
180
     * Calculate surface distance between two points.
181
     */
182 153
    protected static function vincenty(GeographicValue $from, GeographicValue $to, Ellipsoid $ellipsoid): Length
183
    {
184 153
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
185 153
        $b = $ellipsoid->getSemiMinorAxis()->asMetres()->getValue();
186 153
        $f = $ellipsoid->getInverseFlattening();
187 153
        $U1 = atan((1 - $f) * tan($from->getLatitude()->asRadians()->getValue()));
188 153
        $U2 = atan((1 - $f) * tan($to->getLatitude()->asRadians()->getValue()));
189 153
        $L = $to->getLongitude()->subtract($from->getLongitude())->asRadians()->getValue();
190
191 153
        $lambda = $L;
192
        do {
193 153
            $lambdaN = $lambda;
194
195 153
            $sinSigma = sqrt((cos($U2) * sin($lambda)) ** 2 + (cos($U1) * sin($U2) - sin($U1) * cos($U2) * cos($lambda)) ** 2);
196 153
            $cosSigma = sin($U1) * sin($U2) + cos($U1) * cos($U2) * cos($lambda);
197 153
            $sigma = atan2($sinSigma, $cosSigma);
198
199 153
            $sinAlpha = cos($U1) * cos($U2) * sin($lambda) / $sinSigma;
200 153
            $cosSqAlpha = (1 - $sinAlpha ** 2);
201 153
            $cos2SigmaM = $cosSqAlpha ? $cosSigma - (2 * sin($U1) * sin($U2) / $cosSqAlpha) : 0;
202 153
            $C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
203 153
            $lambda = $L + (1 - $C) * $f * $sinAlpha * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (-1 + 2 * $cos2SigmaM ** 2)));
204 153
        } while (abs($lambda - $lambdaN) >= static::ITERATION_CONVERGENCE && abs($lambda) < M_PI);
205
206
        // Antipodal case
207 153
        if (abs($lambda) >= M_PI) {
208 27
            if ($L >= 0) {
209 27
                $LPrime = M_PI - $L;
210
            } else {
211
                $LPrime = -M_PI - $L;
212
            }
213
214 27
            $lambdaPrime = 0;
215 27
            $sigma = M_PI - abs($U1 + $U2);
216 27
            $sinSigma = sin($sigma);
217 27
            $cosSqAlpha = 0.5;
218 27
            $sinAlpha = 0;
219
220
            do {
221 27
                $sinAlphaN = $sinAlpha;
222
223 27
                $C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
224 27
                $cos2SigmaM = cos($sigma) - 2 * sin($U1) * sin($U2) / $cosSqAlpha;
225 27
                $D = (1 - $C) * $f * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * cos($sigma) * (-1 + 2 * $cos2SigmaM ** 2)));
226 27
                $sinAlpha = ($LPrime - $lambdaPrime) / $D;
227 27
                $cosSqAlpha = (1 - $sinAlpha ** 2);
228 27
                $sinLambdaPrime = ($sinAlpha * $sinSigma) / (cos($U1) * cos($U2));
229 27
                $lambdaPrime = self::asin($sinLambdaPrime);
230 27
                $sinSqSigma = (cos($U2) * $sinLambdaPrime) ** 2 + (cos($U1) * sin($U2) + sin($U1) * cos($U2) * cos($lambdaPrime)) ** 2;
231 27
                $sinSigma = sqrt($sinSqSigma);
232 27
            } while (abs($sinAlpha - $sinAlphaN) >= static::ITERATION_CONVERGENCE);
233
        }
234
235 153
        $E = sqrt(1 + (($a ** 2 - $b ** 2) / $b ** 2) * $cosSqAlpha);
236 153
        $F = ($E - 1) / ($E + 1);
237
238 153
        $A = (1 + $F ** 2 / 4) / (1 - $F);
239 153
        $B = $F * (1 - 3 / 8 * $F ** 2);
240
241 153
        $deltaSigma = $B * $sinSigma * ($cos2SigmaM + $B / 4 * ($cosSigma * (-1 + 2 * $cos2SigmaM ** 2) - $B / 6 * $cos2SigmaM * (-3 + 4 * $sinSigma ** 2) * (-3 + 4 * $cos2SigmaM ** 2)));
242
243 153
        return new Metre($b * $A * ($sigma - $deltaSigma));
244
    }
245
246
    /**
247
     * General polynomial.
248
     * @param Coefficient[] $powerCoefficients
249
     */
250 18
    protected function generalPolynomialUnitless(
251
        float $xs,
252
        float $ys,
253
        UnitOfMeasure $ordinate1OfEvaluationPointInSourceCRS,
254
        UnitOfMeasure $ordinate2OfEvaluationPointInSourceCRS,
255
        UnitOfMeasure $ordinate1OfEvaluationPointInTargetCRS,
256
        UnitOfMeasure $ordinate2OfEvaluationPointInTargetCRS,
257
        Scale $scalingFactorForSourceCRSCoordDifferences,
258
        Scale $scalingFactorForTargetCRSCoordDifferences,
259
        Scale $A0,
260
        Scale $B0,
261
        array $powerCoefficients
262
    ): array {
263 18
        $xso = $ordinate1OfEvaluationPointInSourceCRS->getValue();
264 18
        $yso = $ordinate2OfEvaluationPointInSourceCRS->getValue();
265 18
        $xto = $ordinate1OfEvaluationPointInTargetCRS->getValue();
266 18
        $yto = $ordinate2OfEvaluationPointInTargetCRS->getValue();
267
268 18
        $U = $scalingFactorForSourceCRSCoordDifferences->asUnity()->getValue() * ($xs - $xso);
269 18
        $V = $scalingFactorForSourceCRSCoordDifferences->asUnity()->getValue() * ($ys - $yso);
270
271 18
        $mTdX = $A0->getValue();
272 18
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
273 18
            if ($coefficientName[0] === 'A') {
274 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...
275 18
                $mTdX += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
276
            }
277
        }
278
279 18
        $mTdY = $B0->getValue();
280 18
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
281 18
            if ($coefficientName[0] === 'B') {
282 18
                sscanf($coefficientName, 'Bu%dv%d', $uPower, $vPower);
283 18
                $mTdY += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
284
            }
285
        }
286
287 18
        $xt = $xs - $xso + $xto + $mTdX / $scalingFactorForTargetCRSCoordDifferences->asUnity()->getValue();
288 18
        $yt = $ys - $yso + $yto + $mTdY / $scalingFactorForTargetCRSCoordDifferences->asUnity()->getValue();
289
290 18
        return ['xt' => $xt, 'yt' => $yt];
291
    }
292
293
    /**
294
     * Reversible polynomial.
295
     */
296 36
    protected function reversiblePolynomialUnitless(
297
        float $xs,
298
        float $ys,
299
        Angle $ordinate1OfEvaluationPoint,
300
        Angle $ordinate2OfEvaluationPoint,
301
        Scale $scalingFactorForCoordDifferences,
302
        Scale $A0,
303
        Scale $B0,
304
        array $powerCoefficients
305
    ): array {
306 36
        $xo = $ordinate1OfEvaluationPoint->getValue();
307 36
        $yo = $ordinate2OfEvaluationPoint->getValue();
308
309 36
        $U = $scalingFactorForCoordDifferences->asUnity()->getValue() * ($xs - $xo);
310 36
        $V = $scalingFactorForCoordDifferences->asUnity()->getValue() * ($ys - $yo);
311
312 36
        $mTdX = $A0->getValue();
313 36
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
314 36
            if ($coefficientName[0] === 'A') {
315 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...
316 36
                $mTdX += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
317
            }
318
        }
319
320 36
        $mTdY = $B0->getValue();
321 36
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
322 36
            if ($coefficientName[0] === 'B') {
323 36
                sscanf($coefficientName, 'Bu%dv%d', $uPower, $vPower);
324 36
                $mTdY += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
325
            }
326
        }
327
328 36
        $xt = $xs + $mTdX * $scalingFactorForCoordDifferences->asUnity()->getValue();
329 36
        $yt = $ys + $mTdY * $scalingFactorForCoordDifferences->asUnity()->getValue();
330
331 36
        return ['xt' => $xt, 'yt' => $yt];
332
    }
333
334
    /**
335
     * Floating point vagaries mean that it's possible for inputs to be e.g. 1.00000000000001 which makes PHP give a
336
     * silent NaN as output so inputs need to be capped. atan/atan2 are not affected, they seem to cap internally.
337
     */
338
    protected static function acos(float $num): float
339
    {
340
        if ($num > 1.0) {
341
            $num = 1.0;
342
        } elseif ($num < -1) {
343
            $num = -1.0;
344
        }
345
346
        return acos($num);
347
    }
348
349
    /**
350
     * Floating point vagaries mean that it's possible for inputs to be e.g. 1.00000000000001 which makes PHP give a
351
     * silent NaN as output so inputs need to be capped. atan/atan2 are not affected, they seem to cap internally.
352
     */
353 495
    protected static function asin(float $num): float
354
    {
355 495
        if ($num > 1.0) {
356
            $num = 1.0;
357 495
        } elseif ($num < -1.0) {
358
            $num = -1.0;
359
        }
360
361 495
        return asin($num);
362
    }
363
364
    abstract public function getCRS(): CoordinateReferenceSystem;
365
366
    abstract public function getCoordinateEpoch(): ?DateTimeImmutable;
367
368
    abstract public function calculateDistance(self $to): Length;
369
}
370