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

Point::resolveParamsByOperation()   B

Complexity

Conditions 9
Paths 40

Size

Total Lines 39
Code Lines 28

Duplication

Lines 0
Ratio 0 %

Code Coverage

Tests 22
CRAP Score 9.2946

Importance

Changes 2
Bugs 0 Features 0
Metric Value
cc 9
eloc 28
c 2
b 0
f 0
nc 40
nop 3
dl 0
loc 39
ccs 22
cts 26
cp 0.8462
crap 9.2946
rs 8.0555
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