Passed
Push — master ( 2a4317...e180e7 )
by Doug
47:34
created

Point   B

Complexity

Total Complexity 43

Size/Duplication

Total Lines 281
Duplicated Lines 0 %

Test Coverage

Coverage 92.06%

Importance

Changes 8
Bugs 0 Features 0
Metric Value
eloc 146
dl 0
loc 281
ccs 116
cts 126
cp 0.9206
rs 8.96
c 8
b 0
f 0
wmc 43

8 Methods

Rating   Name   Duplication   Size   Complexity  
C resolveParamsByOperation() 0 28 13
A generalPolynomialUnitless() 0 41 5
A reversiblePolynomialUnitless() 0 36 5
A sign() 0 7 2
A acos() 0 9 3
B vincenty() 0 62 9
A asin() 0 9 3
A performOperation() 0 18 3

How to fix   Complexity   

Complex Class

Complex classes like Point often do a lot of different things. To break such a class down, we need to identify a cohesive component within that class. A common approach to find such a component is to look for fields/methods that share the same prefixes, or suffixes.

Once you have determined the fields that belong together, you can apply the Extract Class refactoring. If the component makes sense as a sub-class, Extract Subclass is also a candidate, and is often faster.

While breaking up the class, it is a good idea to analyze how other classes use Point, and based on these observations, apply Extract Interface, too.

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