Passed
Push — master ( 9abb61...fb6dbd )
by Doug
08:29
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 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;
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...
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 Stringable;
44
use function tan;
45
use function ucwords;
46
47
abstract class Point implements Stringable
48
{
49
    protected const ITERATION_CONVERGENCE = 1e-12;
50
51
    /**
52
     * @internal
53
     */
54 19
    public function performOperation(string $srid, CoordinateReferenceSystem $to, bool $inReverse): self
55
    {
56 19
        $operations = self::resolveConcatenatedOperations($srid, $inReverse);
57
58 19
        $point = $this;
59 19
        foreach ($operations as $operationSrid => $operation) {
60 19
            $method = CoordinateOperationMethods::getFunctionName($operation['method']);
61 19
            if (isset($operation['source_crs'])) {
62
                $destCRS = CoordinateReferenceSystem::fromSRID($inReverse ? $operation['source_crs'] : $operation['target_crs']);
63
            } else {
64 19
                $destCRS = $to;
65
            }
66
67 19
            $params = self::resolveParamsByOperation($operationSrid, $operation['method'], $inReverse);
68
69 19
            if (PHP_MAJOR_VERSION >= 8) {
70
                $point = $point->$method($destCRS, ...$params);
71
            } else {
72 19
                $point = $point->$method($destCRS, ...array_values($params));
73
            }
74
        }
75
76 19
        $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
77
78 19
        return $point;
79
    }
80
81 19
    protected static function camelCase(string $string): string
82
    {
83 19
        $string = str_replace([' ', '-'], '', ucwords($string, ' -'));
84 19
        if (!preg_match('/[ABC][uv\d]/', $string)) {
85 19
            $string = lcfirst($string);
86
        }
87
88 19
        return $string;
89
    }
90
91 20
    protected static function resolveConcatenatedOperations(string $operationSrid, bool $inReverse): array
92
    {
93 20
        $operations = [];
94 20
        $operation = CoordinateOperations::getOperationData($operationSrid);
95 20
        if (isset($operation['operations'])) {
96 1
            foreach ($operation['operations'] as $subOperation) {
97 1
                $subOperationData = CoordinateOperations::getOperationData($subOperation['operation']);
98 1
                $subOperationData['source_crs'] = $subOperation['source_crs'];
99 1
                $subOperationData['target_crs'] = $subOperation['target_crs'];
100 1
                $operations[$subOperation['operation']] = $subOperationData;
101
            }
102
        } else {
103 20
            $operations[$operationSrid] = $operation;
104
        }
105
106 20
        if ($inReverse) {
107 10
            $operations = array_reverse($operations, true);
108
        }
109
110 20
        return $operations;
111
    }
112
113 19
    protected static function resolveParamsByOperation(string $operationSrid, string $methodSrid, bool $inReverse): array
114
    {
115 19
        $params = [];
116 19
        $powerCoefficients = [];
117 19
        foreach (CoordinateOperationParams::getParamData($operationSrid) as $paramName => $paramData) {
118 19
            $value = $paramData['value'];
119 19
            if ($inReverse && $paramData['reverses']) {
120 4
                $value *= -1;
121
            }
122 19
            if ($paramData['uom']) {
123 19
                $param = UnitOfMeasureFactory::makeUnit($value, $paramData['uom']);
124
            } else {
125
                $param = $paramData['value'];
126
            }
127 19
            $paramName = static::camelCase($paramName);
128 19
            if (preg_match('/^(Au|Bu)/', $paramName)) {
129 3
                $powerCoefficients[$paramName] = $param;
130
            } else {
131 19
                $params[$paramName] = $param;
132
            }
133
        }
134 19
        if ($powerCoefficients) {
135 3
            $params['powerCoefficients'] = $powerCoefficients;
136
        }
137 19
        if (in_array($methodSrid, [
138
            CoordinateOperationMethods::EPSG_SIMILARITY_TRANSFORMATION,
139
            CoordinateOperationMethods::EPSG_AFFINE_PARAMETRIC_TRANSFORMATION,
140
        ], true)) {
141
            $params['inReverse'] = $inReverse;
142
        }
143
144 19
        return $params;
145
    }
146
147 215
    protected function getAxisByName(string $name): ?Axis
148
    {
149 215
        foreach ($this->getCRS()->getCoordinateSystem()->getAxes() as $axis) {
150 215
            if ($axis->getName() === $name) {
151 215
                return $axis;
152
            }
153
        }
154
155 136
        return null;
156
    }
157
158 10
    protected static function sign(float $number): int
159
    {
160 10
        if ($number < 0) {
161
            return -1;
162
        }
163
164 10
        return 1;
165
    }
166
167
    /**
168
     * Calculate surface distance between two points.
169
     */
170 15
    protected static function vincenty(GeographicValue $from, GeographicValue $to, Ellipsoid $ellipsoid): Length
171
    {
172 15
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
173 15
        $b = $ellipsoid->getSemiMinorAxis()->asMetres()->getValue();
174 15
        $f = $ellipsoid->getInverseFlattening();
175 15
        $U1 = atan((1 - $f) * tan($from->getLatitude()->asRadians()->getValue()));
176 15
        $U2 = atan((1 - $f) * tan($to->getLatitude()->asRadians()->getValue()));
177 15
        $L = $to->getLongitude()->subtract($from->getLongitude())->asRadians()->getValue();
178
179 15
        $lambda = $L;
180
        do {
181 15
            $lambdaN = $lambda;
182
183 15
            $sinSigma = sqrt((cos($U2) * sin($lambda)) ** 2 + (cos($U1) * sin($U2) - sin($U1) * cos($U2) * cos($lambda)) ** 2);
184 15
            $cosSigma = sin($U1) * sin($U2) + cos($U1) * cos($U2) * cos($lambda);
185 15
            $sigma = atan2($sinSigma, $cosSigma);
186
187 15
            $sinAlpha = cos($U1) * cos($U2) * sin($lambda) / $sinSigma;
188 15
            $cosSqAlpha = (1 - $sinAlpha ** 2);
189 15
            $cos2SigmaM = $cosSqAlpha ? $cosSigma - (2 * sin($U1) * sin($U2) / $cosSqAlpha) : 0;
190 15
            $C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
191 15
            $lambda = $L + (1 - $C) * $f * $sinAlpha * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (-1 + 2 * $cos2SigmaM ** 2)));
192 15
        } while (abs($lambda - $lambdaN) >= static::ITERATION_CONVERGENCE && abs($lambda) < M_PI);
193
194
        // Antipodal case
195 15
        if (abs($lambda) >= M_PI) {
196 3
            if ($L >= 0) {
197 3
                $LPrime = M_PI - $L;
198
            } else {
199
                $LPrime = -M_PI - $L;
200
            }
201
202 3
            $lambdaPrime = 0;
203 3
            $sigma = M_PI - abs($U1 + $U2);
204 3
            $sinSigma = sin($sigma);
205 3
            $cosSqAlpha = 0.5;
206 3
            $sinAlpha = 0;
207
208
            do {
209 3
                $sinAlphaN = $sinAlpha;
210
211 3
                $C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
212 3
                $cos2SigmaM = cos($sigma) - 2 * sin($U1) * sin($U2) / $cosSqAlpha;
213 3
                $D = (1 - $C) * $f * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * cos($sigma) * (-1 + 2 * $cos2SigmaM ** 2)));
214 3
                $sinAlpha = ($LPrime - $lambdaPrime) / $D;
215 3
                $cosSqAlpha = (1 - $sinAlpha ** 2);
216 3
                $sinLambdaPrime = ($sinAlpha * $sinSigma) / (cos($U1) * cos($U2));
217 3
                $lambdaPrime = self::asin($sinLambdaPrime);
218 3
                $sinSqSigma = (cos($U2) * $sinLambdaPrime) ** 2 + (cos($U1) * sin($U2) + sin($U1) * cos($U2) * cos($lambdaPrime)) ** 2;
219 3
                $sinSigma = sqrt($sinSqSigma);
220 3
            } while (abs($sinAlpha - $sinAlphaN) >= static::ITERATION_CONVERGENCE);
221
        }
222
223 15
        $E = sqrt(1 + (($a ** 2 - $b ** 2) / $b ** 2) * $cosSqAlpha);
224 15
        $F = ($E - 1) / ($E + 1);
225
226 15
        $A = (1 + $F ** 2 / 4) / (1 - $F);
227 15
        $B = $F * (1 - 3 / 8 * $F ** 2);
228
229 15
        $deltaSigma = $B * $sinSigma * ($cos2SigmaM + $B / 4 * ($cosSigma * (-1 + 2 * $cos2SigmaM ** 2) - $B / 6 * $cos2SigmaM * (-3 + 4 * $sinSigma ** 2) * (-3 + 4 * $cos2SigmaM ** 2)));
230
231 15
        return new Metre($b * $A * ($sigma - $deltaSigma));
232
    }
233
234
    /**
235
     * General polynomial.
236
     * @param Coefficient[] $powerCoefficients
237
     */
238 2
    protected function generalPolynomialUnitless(
239
        float $xs,
240
        float $ys,
241
        UnitOfMeasure $ordinate1OfEvaluationPointInSourceCRS,
242
        UnitOfMeasure $ordinate2OfEvaluationPointInSourceCRS,
243
        UnitOfMeasure $ordinate1OfEvaluationPointInTargetCRS,
244
        UnitOfMeasure $ordinate2OfEvaluationPointInTargetCRS,
245
        Scale $scalingFactorForSourceCRSCoordDifferences,
246
        Scale $scalingFactorForTargetCRSCoordDifferences,
247
        Scale $A0,
248
        Scale $B0,
249
        array $powerCoefficients
250
    ): array {
251 2
        $xso = $ordinate1OfEvaluationPointInSourceCRS->getValue();
252 2
        $yso = $ordinate2OfEvaluationPointInSourceCRS->getValue();
253 2
        $xto = $ordinate1OfEvaluationPointInTargetCRS->getValue();
254 2
        $yto = $ordinate2OfEvaluationPointInTargetCRS->getValue();
255
256 2
        $U = $scalingFactorForSourceCRSCoordDifferences->asUnity()->getValue() * ($xs - $xso);
257 2
        $V = $scalingFactorForSourceCRSCoordDifferences->asUnity()->getValue() * ($ys - $yso);
258
259 2
        $mTdX = $A0->getValue();
260 2
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
261 2
            if ($coefficientName[0] === 'A') {
262 2
                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...
263 2
                $mTdX += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
264
            }
265
        }
266
267 2
        $mTdY = $B0->getValue();
268 2
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
269 2
            if ($coefficientName[0] === 'B') {
270 2
                sscanf($coefficientName, 'Bu%dv%d', $uPower, $vPower);
271 2
                $mTdY += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
272
            }
273
        }
274
275 2
        $xt = $xs - $xso + $xto + $mTdX / $scalingFactorForTargetCRSCoordDifferences->asUnity()->getValue();
276 2
        $yt = $ys - $yso + $yto + $mTdY / $scalingFactorForTargetCRSCoordDifferences->asUnity()->getValue();
277
278 2
        return ['xt' => $xt, 'yt' => $yt];
279
    }
280
281
    /**
282
     * Reversible polynomial.
283
     */
284 4
    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
    ): array {
294 4
        $xo = $ordinate1OfEvaluationPoint->getValue();
295 4
        $yo = $ordinate2OfEvaluationPoint->getValue();
296
297 4
        $U = $scalingFactorForCoordDifferences->asUnity()->getValue() * ($xs - $xo);
298 4
        $V = $scalingFactorForCoordDifferences->asUnity()->getValue() * ($ys - $yo);
299
300 4
        $mTdX = $A0->getValue();
301 4
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
302 4
            if ($coefficientName[0] === 'A') {
303 4
                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...
304 4
                $mTdX += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
305
            }
306
        }
307
308 4
        $mTdY = $B0->getValue();
309 4
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
310 4
            if ($coefficientName[0] === 'B') {
311 4
                sscanf($coefficientName, 'Bu%dv%d', $uPower, $vPower);
312 4
                $mTdY += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
313
            }
314
        }
315
316 4
        $xt = $xs + $mTdX * $scalingFactorForCoordDifferences->asUnity()->getValue();
317 4
        $yt = $ys + $mTdY * $scalingFactorForCoordDifferences->asUnity()->getValue();
318
319 4
        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
     */
341 45
    protected static function asin(float $num): float
342
    {
343 45
        if ($num > 1.0) {
344
            $num = 1.0;
345 45
        } elseif ($num < -1.0) {
346
            $num = -1.0;
347
        }
348
349 45
        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