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; |
|
|
|
|
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
|
236 |
|
public function performOperation(string $srid, CoordinateReferenceSystem $to, bool $inReverse, array $additionalParams = []): self |
63
|
|
|
{ |
64
|
236 |
|
$operation = CoordinateOperations::getOperationData($srid); |
65
|
236 |
|
$method = CoordinateOperationMethods::getFunctionName($operation['method']); |
66
|
236 |
|
$params = self::resolveParamsByOperation($srid, $operation['method'], $inReverse); |
67
|
|
|
|
68
|
236 |
|
if (isset(self::METHODS_REQUIRING_HORIZONTAL_POINT[$operation['method']])) { |
69
|
9 |
|
$params['horizontalPoint'] = $additionalParams['horizontalPoint']; |
70
|
|
|
} |
71
|
|
|
|
72
|
236 |
|
return $this->$method($to, ...$params); |
73
|
|
|
} |
74
|
|
|
|
75
|
236 |
|
protected static function resolveParamsByOperation(string $operationSrid, string $methodSrid, bool $inReverse): array |
76
|
|
|
{ |
77
|
236 |
|
$params = []; |
78
|
236 |
|
$powerCoefficients = []; |
79
|
236 |
|
foreach (CoordinateOperationParams::getParamData($operationSrid) as $paramName => $paramData) { |
80
|
236 |
|
if (isset($paramData['fileProvider'])) { |
81
|
6 |
|
$params[$paramName] = (new $paramData['fileProvider']())->provideGrid(); |
82
|
|
|
} else { |
83
|
236 |
|
if ($inReverse && $paramData['reverses']) { |
84
|
69 |
|
$paramData['value'] *= -1; |
85
|
|
|
} |
86
|
236 |
|
if ($paramData['uom']) { |
87
|
236 |
|
$param = UnitOfMeasureFactory::makeUnit($paramData['value'], $paramData['uom']); |
88
|
|
|
} else { |
89
|
10 |
|
$param = $paramData['value']; |
90
|
|
|
} |
91
|
236 |
|
if (str_starts_with($paramName, 'Au') || str_starts_with($paramName, 'Bu')) { |
92
|
27 |
|
$powerCoefficients[$paramName] = $param; |
93
|
|
|
} else { |
94
|
236 |
|
$params[$paramName] = $param; |
95
|
|
|
} |
96
|
|
|
} |
97
|
|
|
} |
98
|
236 |
|
if ($powerCoefficients) { |
99
|
27 |
|
$params['powerCoefficients'] = $powerCoefficients; |
100
|
|
|
} |
101
|
236 |
|
if (isset(self::METHODS_THAT_REQUIRE_DIRECTION[$methodSrid])) { |
102
|
2 |
|
$params['inReverse'] = $inReverse; |
103
|
|
|
} |
104
|
|
|
|
105
|
236 |
|
return $params; |
106
|
|
|
} |
107
|
|
|
|
108
|
90 |
|
protected static function sign(float $number): int |
109
|
|
|
{ |
110
|
90 |
|
if ($number < 0) { |
111
|
|
|
return -1; |
112
|
|
|
} |
113
|
|
|
|
114
|
90 |
|
return 1; |
115
|
|
|
} |
116
|
|
|
|
117
|
|
|
/** |
118
|
|
|
* Calculate surface distance between two points. |
119
|
|
|
*/ |
120
|
153 |
|
protected static function vincenty(GeographicValue $from, GeographicValue $to, Ellipsoid $ellipsoid): Length |
121
|
|
|
{ |
122
|
153 |
|
$a = $ellipsoid->getSemiMajorAxis()->asMetres()->getValue(); |
123
|
153 |
|
$b = $ellipsoid->getSemiMinorAxis()->asMetres()->getValue(); |
124
|
153 |
|
$f = $ellipsoid->getInverseFlattening(); |
125
|
153 |
|
$U1 = atan((1 - $f) * tan($from->getLatitude()->asRadians()->getValue())); |
126
|
153 |
|
$U2 = atan((1 - $f) * tan($to->getLatitude()->asRadians()->getValue())); |
127
|
153 |
|
$L = $to->getLongitude()->subtract($from->getLongitude())->asRadians()->getValue(); |
128
|
|
|
|
129
|
153 |
|
$lambda = $L; |
130
|
|
|
do { |
131
|
153 |
|
$lambdaN = $lambda; |
132
|
|
|
|
133
|
153 |
|
$sinSigma = sqrt((cos($U2) * sin($lambda)) ** 2 + (cos($U1) * sin($U2) - sin($U1) * cos($U2) * cos($lambda)) ** 2); |
134
|
153 |
|
$cosSigma = sin($U1) * sin($U2) + cos($U1) * cos($U2) * cos($lambda); |
135
|
153 |
|
$sigma = atan2($sinSigma, $cosSigma); |
136
|
|
|
|
137
|
153 |
|
$sinAlpha = cos($U1) * cos($U2) * sin($lambda) / $sinSigma; |
138
|
153 |
|
$cosSqAlpha = (1 - $sinAlpha ** 2); |
139
|
153 |
|
$cos2SigmaM = $cosSqAlpha ? $cosSigma - (2 * sin($U1) * sin($U2) / $cosSqAlpha) : 0; |
140
|
153 |
|
$C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha)); |
141
|
153 |
|
$lambda = $L + (1 - $C) * $f * $sinAlpha * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (-1 + 2 * $cos2SigmaM ** 2))); |
142
|
153 |
|
} while (abs($lambda - $lambdaN) >= static::ITERATION_CONVERGENCE_FORMULA && abs($lambda) < M_PI); |
143
|
|
|
|
144
|
|
|
// Antipodal case |
145
|
153 |
|
if (abs($lambda) >= M_PI) { |
146
|
27 |
|
if ($L >= 0) { |
147
|
27 |
|
$LPrime = M_PI - $L; |
148
|
|
|
} else { |
149
|
|
|
$LPrime = -M_PI - $L; |
150
|
|
|
} |
151
|
|
|
|
152
|
27 |
|
$lambdaPrime = 0; |
153
|
27 |
|
$sigma = M_PI - abs($U1 + $U2); |
154
|
27 |
|
$sinSigma = sin($sigma); |
155
|
27 |
|
$cosSqAlpha = 0.5; |
156
|
27 |
|
$sinAlpha = 0; |
157
|
|
|
|
158
|
|
|
do { |
159
|
27 |
|
$sinAlphaN = $sinAlpha; |
160
|
|
|
|
161
|
27 |
|
$C = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha)); |
162
|
27 |
|
$cos2SigmaM = cos($sigma) - 2 * sin($U1) * sin($U2) / $cosSqAlpha; |
163
|
27 |
|
$D = (1 - $C) * $f * ($sigma + $C * $sinSigma * ($cos2SigmaM + $C * cos($sigma) * (-1 + 2 * $cos2SigmaM ** 2))); |
164
|
27 |
|
$sinAlpha = ($LPrime - $lambdaPrime) / $D; |
165
|
27 |
|
$cosSqAlpha = (1 - $sinAlpha ** 2); |
166
|
27 |
|
$sinLambdaPrime = ($sinAlpha * $sinSigma) / (cos($U1) * cos($U2)); |
167
|
27 |
|
$lambdaPrime = self::asin($sinLambdaPrime); |
168
|
27 |
|
$sinSqSigma = (cos($U2) * $sinLambdaPrime) ** 2 + (cos($U1) * sin($U2) + sin($U1) * cos($U2) * cos($lambdaPrime)) ** 2; |
169
|
27 |
|
$sinSigma = sqrt($sinSqSigma); |
170
|
27 |
|
} while (abs($sinAlpha - $sinAlphaN) >= static::ITERATION_CONVERGENCE_FORMULA); |
171
|
|
|
} |
172
|
|
|
|
173
|
153 |
|
$E = sqrt(1 + (($a ** 2 - $b ** 2) / $b ** 2) * $cosSqAlpha); |
174
|
153 |
|
$F = ($E - 1) / ($E + 1); |
175
|
|
|
|
176
|
153 |
|
$A = (1 + $F ** 2 / 4) / (1 - $F); |
177
|
153 |
|
$B = $F * (1 - 3 / 8 * $F ** 2); |
178
|
|
|
|
179
|
153 |
|
$deltaSigma = $B * $sinSigma * ($cos2SigmaM + $B / 4 * ($cosSigma * (-1 + 2 * $cos2SigmaM ** 2) - $B / 6 * $cos2SigmaM * (-3 + 4 * $sinSigma ** 2) * (-3 + 4 * $cos2SigmaM ** 2))); |
180
|
|
|
|
181
|
153 |
|
return new Metre($b * $A * ($sigma - $deltaSigma)); |
182
|
|
|
} |
183
|
|
|
|
184
|
|
|
/** |
185
|
|
|
* General polynomial. |
186
|
|
|
* @param Coefficient[] $powerCoefficients |
187
|
|
|
*/ |
188
|
18 |
|
protected function generalPolynomialUnitless( |
189
|
|
|
float $xs, |
190
|
|
|
float $ys, |
191
|
|
|
UnitOfMeasure $ordinate1OfEvaluationPointInSourceCRS, |
192
|
|
|
UnitOfMeasure $ordinate2OfEvaluationPointInSourceCRS, |
193
|
|
|
UnitOfMeasure $ordinate1OfEvaluationPointInTargetCRS, |
194
|
|
|
UnitOfMeasure $ordinate2OfEvaluationPointInTargetCRS, |
195
|
|
|
Scale $scalingFactorForSourceCRSCoordDifferences, |
196
|
|
|
Scale $scalingFactorForTargetCRSCoordDifferences, |
197
|
|
|
Scale $A0, |
198
|
|
|
Scale $B0, |
199
|
|
|
array $powerCoefficients |
200
|
|
|
): array { |
201
|
18 |
|
$xso = $ordinate1OfEvaluationPointInSourceCRS->getValue(); |
202
|
18 |
|
$yso = $ordinate2OfEvaluationPointInSourceCRS->getValue(); |
203
|
18 |
|
$xto = $ordinate1OfEvaluationPointInTargetCRS->getValue(); |
204
|
18 |
|
$yto = $ordinate2OfEvaluationPointInTargetCRS->getValue(); |
205
|
|
|
|
206
|
18 |
|
$U = $scalingFactorForSourceCRSCoordDifferences->asUnity()->getValue() * ($xs - $xso); |
207
|
18 |
|
$V = $scalingFactorForSourceCRSCoordDifferences->asUnity()->getValue() * ($ys - $yso); |
208
|
|
|
|
209
|
18 |
|
$mTdX = $A0->getValue(); |
210
|
18 |
|
foreach ($powerCoefficients as $coefficientName => $coefficientValue) { |
211
|
18 |
|
if ($coefficientName[0] === 'A') { |
212
|
18 |
|
sscanf($coefficientName, 'Au%dv%d', $uPower, $vPower); |
|
|
|
|
213
|
18 |
|
$mTdX += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower; |
214
|
|
|
} |
215
|
|
|
} |
216
|
|
|
|
217
|
18 |
|
$mTdY = $B0->getValue(); |
218
|
18 |
|
foreach ($powerCoefficients as $coefficientName => $coefficientValue) { |
219
|
18 |
|
if ($coefficientName[0] === 'B') { |
220
|
18 |
|
sscanf($coefficientName, 'Bu%dv%d', $uPower, $vPower); |
221
|
18 |
|
$mTdY += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower; |
222
|
|
|
} |
223
|
|
|
} |
224
|
|
|
|
225
|
18 |
|
$xt = $xs - $xso + $xto + $mTdX / $scalingFactorForTargetCRSCoordDifferences->asUnity()->getValue(); |
226
|
18 |
|
$yt = $ys - $yso + $yto + $mTdY / $scalingFactorForTargetCRSCoordDifferences->asUnity()->getValue(); |
227
|
|
|
|
228
|
18 |
|
return ['xt' => $xt, 'yt' => $yt]; |
229
|
|
|
} |
230
|
|
|
|
231
|
|
|
/** |
232
|
|
|
* Reversible polynomial. |
233
|
|
|
*/ |
234
|
36 |
|
protected function reversiblePolynomialUnitless( |
235
|
|
|
float $xs, |
236
|
|
|
float $ys, |
237
|
|
|
Angle $ordinate1OfEvaluationPoint, |
238
|
|
|
Angle $ordinate2OfEvaluationPoint, |
239
|
|
|
Scale $scalingFactorForCoordDifferences, |
240
|
|
|
Scale $A0, |
241
|
|
|
Scale $B0, |
242
|
|
|
array $powerCoefficients |
243
|
|
|
): array { |
244
|
36 |
|
$xo = $ordinate1OfEvaluationPoint->getValue(); |
245
|
36 |
|
$yo = $ordinate2OfEvaluationPoint->getValue(); |
246
|
|
|
|
247
|
36 |
|
$U = $scalingFactorForCoordDifferences->asUnity()->getValue() * ($xs - $xo); |
248
|
36 |
|
$V = $scalingFactorForCoordDifferences->asUnity()->getValue() * ($ys - $yo); |
249
|
|
|
|
250
|
36 |
|
$mTdX = $A0->getValue(); |
251
|
36 |
|
foreach ($powerCoefficients as $coefficientName => $coefficientValue) { |
252
|
36 |
|
if ($coefficientName[0] === 'A') { |
253
|
36 |
|
sscanf($coefficientName, 'Au%dv%d', $uPower, $vPower); |
|
|
|
|
254
|
36 |
|
$mTdX += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower; |
255
|
|
|
} |
256
|
|
|
} |
257
|
|
|
|
258
|
36 |
|
$mTdY = $B0->getValue(); |
259
|
36 |
|
foreach ($powerCoefficients as $coefficientName => $coefficientValue) { |
260
|
36 |
|
if ($coefficientName[0] === 'B') { |
261
|
36 |
|
sscanf($coefficientName, 'Bu%dv%d', $uPower, $vPower); |
262
|
36 |
|
$mTdY += $coefficientValue->getValue() * $U ** $uPower * $V ** $vPower; |
263
|
|
|
} |
264
|
|
|
} |
265
|
|
|
|
266
|
36 |
|
$xt = $xs + $mTdX * $scalingFactorForCoordDifferences->asUnity()->getValue(); |
267
|
36 |
|
$yt = $ys + $mTdY * $scalingFactorForCoordDifferences->asUnity()->getValue(); |
268
|
|
|
|
269
|
36 |
|
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
|
|
|
* 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
|
|
|
/** |
288
|
|
|
* Floating point vagaries mean that it's possible for inputs to be e.g. 1.00000000000001 which makes PHP give a |
289
|
|
|
* silent NaN as output so inputs need to be capped. atan/atan2 are not affected, they seem to cap internally. |
290
|
|
|
*/ |
291
|
482 |
|
protected static function asin(float $num): float |
292
|
|
|
{ |
293
|
482 |
|
if ($num > 1.0) { |
294
|
|
|
$num = 1.0; |
295
|
482 |
|
} elseif ($num < -1.0) { |
296
|
|
|
$num = -1.0; |
297
|
|
|
} |
298
|
|
|
|
299
|
482 |
|
return asin($num); |
300
|
|
|
} |
301
|
|
|
|
302
|
|
|
abstract public function getCRS(): CoordinateReferenceSystem; |
303
|
|
|
|
304
|
|
|
abstract public function getCoordinateEpoch(): ?DateTimeImmutable; |
305
|
|
|
|
306
|
|
|
abstract public function calculateDistance(self $to): Length; |
307
|
|
|
} |
308
|
|
|
|
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:For further information see https://scrutinizer-ci.com/docs/tools/php/php-scrutinizer/#list-dependency-paths