Passed
Push — extents ( f35511...ba1d5f )
by Doug
61:44
created

Point::generalPolynomialUnitless()   A

Complexity

Conditions 5
Paths 9

Size

Total Lines 41
Code Lines 19

Duplication

Lines 0
Ratio 0 %

Code Coverage

Tests 22
CRAP Score 5.002

Importance

Changes 0
Metric Value
cc 5
eloc 19
nc 9
nop 11
dl 0
loc 41
ccs 22
cts 23
cp 0.9565
crap 5.002
rs 9.3222
c 0
b 0
f 0

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