Passed
Push — master ( bbeaaa...9ed6c0 )
by Doug
25:29
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 DateTimeImmutable;
12
use PHPCoord\CoordinateOperation\CoordinateOperationMethods;
13
use PHPCoord\CoordinateOperation\CoordinateOperations;
14
use PHPCoord\CoordinateOperation\GeographicValue;
15
use PHPCoord\CoordinateReferenceSystem\Compound;
16
use PHPCoord\CoordinateReferenceSystem\CoordinateReferenceSystem;
17
use PHPCoord\CoordinateReferenceSystem\Geocentric;
18
use PHPCoord\CoordinateReferenceSystem\Geographic2D;
19
use PHPCoord\CoordinateReferenceSystem\Geographic3D;
20
use PHPCoord\CoordinateReferenceSystem\Projected;
0 ignored issues
show
Bug introduced by
The type PHPCoord\CoordinateReferenceSystem\Projected 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...
21
use PHPCoord\CoordinateReferenceSystem\Vertical;
22
use PHPCoord\Datum\Ellipsoid;
23
use PHPCoord\UnitOfMeasure\Angle\Angle;
24
use PHPCoord\UnitOfMeasure\Length\Length;
25
use PHPCoord\UnitOfMeasure\Length\Metre;
26
use PHPCoord\UnitOfMeasure\Scale\Coefficient;
27
use PHPCoord\UnitOfMeasure\Scale\Scale;
28
use PHPCoord\UnitOfMeasure\UnitOfMeasure;
29
use PHPCoord\UnitOfMeasure\UnitOfMeasureFactory;
30
use Stringable;
31
32
use function abs;
33
use function acos;
34
use function asin;
35
use function atan;
36
use function atan2;
37
use function cos;
38
use function sin;
39
use function sqrt;
40
use function sscanf;
41
use function str_starts_with;
42
use function tan;
43
44
use const M_PI;
45
46
abstract class Point implements Stringable
47
{
48
    protected const ITERATION_CONVERGENCE_FORMULA = 1e-10;
49
    protected const ITERATION_CONVERGENCE_GRID = 0.0001;
50
    protected const METHODS_REQUIRING_HORIZONTAL_POINT = [
51
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_AND_SLOPE => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_AND_SLOPE,
52
        CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019 => CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019,
53
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX,
54
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_PL_TXT => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_PL_TXT,
55
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_BEV_AT => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_BEV_AT,
56
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_NRCAN_BYN => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_NRCAN_BYN,
57
    ];
58
    protected const METHODS_THAT_REQUIRE_DIRECTION = [
59
        CoordinateOperationMethods::EPSG_SIMILARITY_TRANSFORMATION => CoordinateOperationMethods::EPSG_SIMILARITY_TRANSFORMATION,
60
        CoordinateOperationMethods::EPSG_AFFINE_PARAMETRIC_TRANSFORMATION => CoordinateOperationMethods::EPSG_AFFINE_PARAMETRIC_TRANSFORMATION,
61
        CoordinateOperationMethods::EPSG_NADCON5_2D => CoordinateOperationMethods::EPSG_NADCON5_2D,
62
        CoordinateOperationMethods::EPSG_NADCON5_3D => CoordinateOperationMethods::EPSG_NADCON5_3D,
63
        CoordinateOperationMethods::EPSG_NTV2 => CoordinateOperationMethods::EPSG_NTV2,
64
        CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019 => CoordinateOperationMethods::EPSG_ZERO_TIDE_HEIGHT_TO_MEAN_TIDE_HEIGHT_EVRF2019,
65
        CoordinateOperationMethods::EPSG_GEOCENTRIC_TRANSLATION_BY_GRID_INTERPOLATION_IGN => CoordinateOperationMethods::EPSG_GEOCENTRIC_TRANSLATION_BY_GRID_INTERPOLATION_IGN,
66
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_GTX,
67
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_PL_TXT => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_PL_TXT,
68
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_BEV_AT => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_BEV_AT,
69
        CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_NRCAN_BYN => CoordinateOperationMethods::EPSG_VERTICAL_OFFSET_BY_GRID_INTERPOLATION_NRCAN_BYN,
70
    ];
71
72
    private static array $gridCache = [];
73
74
    /**
75
     * @internal
76
     */
77 5224
    public function performOperation(string $srid, Compound|Geocentric|Geographic2D|Geographic3D|Projected|Vertical $to, bool $inReverse, array $additionalParams = []): self
78
    {
79 5224
        $operation = CoordinateOperations::getOperationData($srid);
80
81 5224
        if ($operation['method'] === CoordinateOperationMethods::EPSG_ALIAS) {
82
            $point = clone $this;
83
            $point->crs = $to;
0 ignored issues
show
Bug Best Practice introduced by
The property crs does not exist. Although not strictly required by PHP, it is generally a best practice to declare properties explicitly.
Loading history...
84
85
            return $point;
86
        } else {
87 5224
            $method = CoordinateOperationMethods::getFunctionName($operation['method']);
88 5224
            $params = self::resolveParamsByOperation($srid, $operation['method'], $inReverse);
89
90 5224
            if (isset(self::METHODS_REQUIRING_HORIZONTAL_POINT[$operation['method']])) {
91 9
                $params['horizontalPoint'] = $additionalParams['horizontalPoint'];
92
            }
93
94 5224
            return $this->$method($to, ...$params);
95
        }
96
    }
97
98 5224
    protected static function resolveParamsByOperation(string $operationSrid, string $methodSrid, bool $inReverse): array
99
    {
100 5224
        $params = [];
101 5224
        $powerCoefficients = [];
102 5224
        foreach (CoordinateOperations::getParamData($operationSrid) as $paramName => $paramData) {
103 5224
            if (isset($paramData['fileProvider'])) {
104 12
                $params[$paramName] = static::$gridCache[$paramData['fileProvider']] ??= (new $paramData['fileProvider']())->provideGrid();
0 ignored issues
show
Bug introduced by
Since $gridCache is declared private, accessing it with static will lead to errors in possible sub-classes; you can either use self, or increase the visibility of $gridCache to at least protected.
Loading history...
105
            } else {
106 5224
                if ($inReverse && $paramData['reverses']) {
107 528
                    $paramData['value'] *= -1;
108
                }
109 5224
                if ($paramData['uom']) {
110 5224
                    $param = UnitOfMeasureFactory::makeUnit($paramData['value'], $paramData['uom']);
111
                } else {
112 11
                    $param = $paramData['value'];
113
                }
114 5224
                if (str_starts_with($paramName, 'Au') || str_starts_with($paramName, 'Bu')) {
115 27
                    $powerCoefficients[$paramName] = $param;
116
                } else {
117 5224
                    $params[$paramName] = $param;
118
                }
119
            }
120
        }
121 5224
        if ($powerCoefficients) {
122 27
            $params['powerCoefficients'] = $powerCoefficients;
123
        }
124 5224
        if (isset(self::METHODS_THAT_REQUIRE_DIRECTION[$methodSrid])) {
125 6
            $params['inReverse'] = $inReverse;
126
        }
127
128 5224
        return $params;
129
    }
130
131 802
    protected static function sign(float $number): int
132
    {
133 802
        if ($number < 0) {
134
            return -1;
135
        }
136
137 802
        return 1;
138
    }
139
140
    /**
141
     * Calculate surface distance between two points.
142
     */
143 162
    protected static function vincenty(GeographicValue $from, GeographicValue $to, Ellipsoid $ellipsoid): Length
144
    {
145 162
        $a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue();
146 162
        $b = $ellipsoid->getSemiMinorAxis()->asMetres()->getValue();
147 162
        $f = $ellipsoid->getFlattening();
148 162
        $U1 = atan((1 - $f) * tan($from->getLatitude()->asRadians()->getValue()));
149 162
        $U2 = atan((1 - $f) * tan($to->getLatitude()->asRadians()->getValue()));
150 162
        $L = $to->getLongitude()->subtract($from->getLongitude())->asRadians()->getValue();
151
152 162
        $lambda = $L;
153
        do {
154 162
            $lambdaN = $lambda;
155
156 162
            $sinSigma = sqrt((cos($U2) * sin($lambda)) ** 2 + (cos($U1) * sin($U2) - sin($U1) * cos($U2) * cos($lambda)) ** 2);
157 162
            $cosSigma = sin($U1) * sin($U2) + cos($U1) * cos($U2) * cos($lambda);
158 162
            $sigma = atan2($sinSigma, $cosSigma);
159
160 162
            $sinAlpha = $sinSigma ? (cos($U1) * cos($U2) * sin($lambda) / $sinSigma) : 0;
161 162
            $cosSqAlpha = (1 - $sinAlpha ** 2);
162 162
            $cos2SigmaM = $cosSqAlpha ? $cosSigma - (2 * sin($U1) * sin($U2) / $cosSqAlpha) : 0;
163 162
            $C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
164 162
            $lambda = $L + (1 - $C) * $f * $sinAlpha * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (-1 + 2 * $cos2SigmaM ** 2)));
165 162
        } while (abs($lambda - $lambdaN) >= static::ITERATION_CONVERGENCE_FORMULA && abs($lambda) < M_PI);
166
167
        // Antipodal case
168 162
        if (abs($lambda) >= M_PI) {
169 27
            if ($L >= 0) {
170 27
                $LPrime = M_PI - $L;
171
            } else {
172
                $LPrime = -M_PI - $L;
173
            }
174
175 27
            $lambdaPrime = 0;
176 27
            $sigma = M_PI - abs($U1 + $U2);
177 27
            $sinSigma = sin($sigma);
178 27
            $cosSqAlpha = 0.5;
179 27
            $sinAlpha = 0;
180
181
            do {
182 27
                $sinAlphaN = $sinAlpha;
183
184 27
                $C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
185 27
                $cos2SigmaM = $cosSqAlpha ? cos($sigma) - 2 * sin($U1) * sin($U2) / $cosSqAlpha : 0;
186 27
                $D = (1 - $C) * $f * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * cos($sigma) * (-1 + 2 * $cos2SigmaM ** 2)));
187 27
                $sinAlpha = ($LPrime - $lambdaPrime) / $D;
188 27
                $cosSqAlpha = (1 - $sinAlpha ** 2);
189 27
                $sinLambdaPrime = ($sinAlpha * $sinSigma) / (cos($U1) * cos($U2));
190 27
                $lambdaPrime = self::asin($sinLambdaPrime);
191 27
                $sinSqSigma = (cos($U2) * $sinLambdaPrime) ** 2 + (cos($U1) * sin($U2) + sin($U1) * cos($U2) * cos($lambdaPrime)) ** 2;
192 27
                $sinSigma = sqrt($sinSqSigma);
193 27
            } while (abs($sinAlpha - $sinAlphaN) >= static::ITERATION_CONVERGENCE_FORMULA);
194
        }
195
196 162
        $E = sqrt(1 + (($a ** 2 - $b ** 2) / $b ** 2) * $cosSqAlpha);
197 162
        $F = ($E - 1) / ($E + 1);
198
199 162
        $A = (1 + $F ** 2 / 4) / (1 - $F);
200 162
        $B = $F * (1 - 3 / 8 * $F ** 2);
201
202 162
        $deltaSigma = $B * $sinSigma * ($cos2SigmaM + $B / 4 * ($cosSigma * (-1 + 2 * $cos2SigmaM ** 2) - $B / 6 * $cos2SigmaM * (-3 + 4 * $sinSigma ** 2) * (-3 + 4 * $cos2SigmaM ** 2)));
203
204 162
        return new Metre($b * $A * ($sigma - $deltaSigma));
205
    }
206
207
    /**
208
     * General polynomial.
209
     * @param Coefficient[] $powerCoefficients
210
     */
211 18
    protected function generalPolynomialUnitless(
212
        float $xs,
213
        float $ys,
214
        UnitOfMeasure $ordinate1OfEvaluationPointInSourceCRS,
215
        UnitOfMeasure $ordinate2OfEvaluationPointInSourceCRS,
216
        UnitOfMeasure $ordinate1OfEvaluationPointInTargetCRS,
217
        UnitOfMeasure $ordinate2OfEvaluationPointInTargetCRS,
218
        Scale $scalingFactorForSourceCRSCoordDifferences,
219
        Scale $scalingFactorForTargetCRSCoordDifferences,
220
        Scale $A0,
221
        Scale $B0,
222
        array $powerCoefficients
223
    ): array {
224 18
        $xso = $ordinate1OfEvaluationPointInSourceCRS->getValue();
225 18
        $yso = $ordinate2OfEvaluationPointInSourceCRS->getValue();
226 18
        $xto = $ordinate1OfEvaluationPointInTargetCRS->getValue();
227 18
        $yto = $ordinate2OfEvaluationPointInTargetCRS->getValue();
228
229 18
        $U = $scalingFactorForSourceCRSCoordDifferences->asUnity()->getValue() * ($xs - $xso);
230 18
        $V = $scalingFactorForSourceCRSCoordDifferences->asUnity()->getValue() * ($ys - $yso);
231
232 18
        $mTdX = $A0->getValue();
233 18
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
234 18
            if ($coefficientName[0] === 'A') {
235 18
                sscanf($coefficientName, 'Au%dv%d', $uPower, $vPower);
236 18
                $mTdX += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
237
            }
238
        }
239
240 18
        $mTdY = $B0->getValue();
241 18
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
242 18
            if ($coefficientName[0] === 'B') {
243 18
                sscanf($coefficientName, 'Bu%dv%d', $uPower, $vPower);
244 18
                $mTdY += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
245
            }
246
        }
247
248 18
        $xt = $xs - $xso + $xto + $mTdX / $scalingFactorForTargetCRSCoordDifferences->asUnity()->getValue();
249 18
        $yt = $ys - $yso + $yto + $mTdY / $scalingFactorForTargetCRSCoordDifferences->asUnity()->getValue();
250
251 18
        return ['xt' => $xt, 'yt' => $yt];
252
    }
253
254
    /**
255
     * Reversible polynomial.
256
     */
257 36
    protected function reversiblePolynomialUnitless(
258
        float $xs,
259
        float $ys,
260
        Angle $ordinate1OfEvaluationPoint,
261
        Angle $ordinate2OfEvaluationPoint,
262
        Scale $scalingFactorForCoordDifferences,
263
        Scale $A0,
264
        Scale $B0,
265
        array $powerCoefficients
266
    ): array {
267 36
        $xo = $ordinate1OfEvaluationPoint->getValue();
268 36
        $yo = $ordinate2OfEvaluationPoint->getValue();
269
270 36
        $U = $scalingFactorForCoordDifferences->asUnity()->getValue() * ($xs - $xo);
271 36
        $V = $scalingFactorForCoordDifferences->asUnity()->getValue() * ($ys - $yo);
272
273 36
        $mTdX = $A0->getValue();
274 36
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
275 36
            if ($coefficientName[0] === 'A') {
276 36
                sscanf($coefficientName, 'Au%dv%d', $uPower, $vPower);
277 36
                $mTdX += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
278
            }
279
        }
280
281 36
        $mTdY = $B0->getValue();
282 36
        foreach ($powerCoefficients as $coefficientName => $coefficientValue) {
283 36
            if ($coefficientName[0] === 'B') {
284 36
                sscanf($coefficientName, 'Bu%dv%d', $uPower, $vPower);
285 36
                $mTdY += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower;
286
            }
287
        }
288
289 36
        $xt = $xs + $mTdX * $scalingFactorForCoordDifferences->asUnity()->getValue();
290 36
        $yt = $ys + $mTdY * $scalingFactorForCoordDifferences->asUnity()->getValue();
291
292 36
        return ['xt' => $xt, 'yt' => $yt];
293
    }
294
295
    /**
296
     * Floating point vagaries mean that it's possible for inputs to be e.g. 1.00000000000001 which makes PHP give a
297
     * silent NaN as output so inputs need to be capped. atan/atan2 are not affected, they seem to cap internally.
298
     */
299
    protected static function acos(float $num): float
300
    {
301
        if ($num > 1.0) {
302
            $num = 1.0;
303
        } elseif ($num < -1) {
304
            $num = -1.0;
305
        }
306
307
        return acos($num);
308
    }
309
310
    /**
311
     * Floating point vagaries mean that it's possible for inputs to be e.g. 1.00000000000001 which makes PHP give a
312
     * silent NaN as output so inputs need to be capped. atan/atan2 are not affected, they seem to cap internally.
313
     */
314 2652
    protected static function asin(float $num): float
315
    {
316 2652
        if ($num > 1.0) {
317
            $num = 1.0;
318 2652
        } elseif ($num < -1.0) {
319
            $num = -1.0;
320
        }
321
322 2652
        return asin($num);
323
    }
324
325
    abstract public function getCRS(): CoordinateReferenceSystem;
326
327
    abstract public function getCoordinateEpoch(): ?DateTimeImmutable;
328
329
    abstract public function calculateDistance(self $to): Length;
330
}
331