Passed
Push — 4.0.x ( 414ad0...5c8c8d )
by Doug
05:04 queued 12s
created

Point::vincenty()   B

Complexity

Conditions 7
Paths 6

Size

Total Lines 62
Code Lines 46

Duplication

Lines 0
Ratio 0 %

Code Coverage

Tests 43
CRAP Score 7.0005

Importance

Changes 0
Metric Value
cc 7
eloc 46
c 0
b 0
f 0
nc 6
nop 3
dl 0
loc 62
ccs 43
cts 44
cp 0.9773
crap 7.0005
rs 8.2448

How to fix   Long Method   

Long Method

Small methods make your code easier to understand, in particular if combined with a good name. Besides, if your method is small, finding a good name is usually much easier.

For example, if you find yourself adding comments to a method's body, this is usually a good sign to extract the commented part to a new method, and use the comment as a starting point when coming up with a good name for this new method.

Commonly applied refactorings include:

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