Passed
Push — master ( fefee4...ec3bb4 )
by Doug
39:04
created

GeocentricPoint::geographicGeocentric()   A

Complexity

Conditions 2
Paths 1

Size

Total Lines 7
Code Lines 3

Duplication

Lines 0
Ratio 0 %

Code Coverage

Tests 4
CRAP Score 2

Importance

Changes 0
Metric Value
eloc 3
c 0
b 0
f 0
dl 0
loc 7
ccs 4
cts 4
cp 1
rs 10
cc 2
nc 1
nop 1
crap 2
1
<?php
2
/**
3
 * PHPCoord.
4
 *
5
 * @author Doug Wright
6
 */
7
declare(strict_types=1);
8
9
namespace PHPCoord;
10
11
use DateTime;
12
use DateTimeImmutable;
13
use DateTimeInterface;
14
use PHPCoord\CoordinateOperation\AutoConversion;
15
use PHPCoord\CoordinateOperation\ConvertiblePoint;
16
use PHPCoord\CoordinateOperation\GeocentricValue;
17
use PHPCoord\CoordinateOperation\GeographicValue;
18
use PHPCoord\CoordinateReferenceSystem\Geocentric;
19
use PHPCoord\CoordinateReferenceSystem\Geographic;
20
use PHPCoord\CoordinateReferenceSystem\Geographic2D;
21
use PHPCoord\CoordinateReferenceSystem\Geographic3D;
22
use PHPCoord\CoordinateSystem\Axis;
23
use PHPCoord\Exception\InvalidCoordinateException;
24
use PHPCoord\Exception\InvalidCoordinateReferenceSystemException;
25
use PHPCoord\UnitOfMeasure\Angle\Angle;
26
use PHPCoord\UnitOfMeasure\Angle\Radian;
27
use PHPCoord\UnitOfMeasure\Length\Length;
28
use PHPCoord\UnitOfMeasure\Length\Metre;
29
use PHPCoord\UnitOfMeasure\Rate;
30
use PHPCoord\UnitOfMeasure\Scale\Scale;
31
use PHPCoord\UnitOfMeasure\Scale\Unity;
32
use PHPCoord\UnitOfMeasure\Time\Time;
33
use PHPCoord\UnitOfMeasure\Time\Year;
34
35
use function abs;
36
use function sprintf;
37
38
/**
39
 * Coordinate representing a point in ECEF geocentric form.
40
 */
41
class GeocentricPoint extends Point implements ConvertiblePoint
42
{
43
    use AutoConversion;
44
45
    /**
46
     * X co-ordinate.
47
     */
48
    protected Length $x;
49
50
    /**
51
     * Y co-ordinate.
52
     */
53
    protected Length $y;
54
55
    /**
56
     * Z co-ordinate.
57
     */
58
    protected Length $z;
59
60
    /**
61
     * Coordinate reference system.
62
     */
63
    protected Geocentric $crs;
64
65
    /**
66
     * Coordinate epoch (date for which the specified coordinates represented this point).
67
     */
68
    protected ?DateTimeImmutable $epoch;
69
70 522
    protected function __construct(Geocentric $crs, Length $x, Length $y, Length $z, ?DateTimeInterface $epoch = null)
71
    {
72 522
        $this->crs = $crs;
73 522
        $this->x = $x::convert($x, $this->crs->getCoordinateSystem()->getAxisByName(Axis::GEOCENTRIC_X)->getUnitOfMeasureId());
74 522
        $this->y = $y::convert($y, $this->crs->getCoordinateSystem()->getAxisByName(Axis::GEOCENTRIC_Y)->getUnitOfMeasureId());
75 522
        $this->z = $z::convert($z, $this->crs->getCoordinateSystem()->getAxisByName(Axis::GEOCENTRIC_Z)->getUnitOfMeasureId());
76
77 522
        if ($epoch instanceof DateTime) {
78 54
            $epoch = DateTimeImmutable::createFromMutable($epoch);
79
        }
80 522
        $this->epoch = $epoch;
81
    }
82
83
    /**
84
     * @param Length $x refer to CRS for preferred unit of measure, but any length unit accepted
85
     * @param Length $y refer to CRS for preferred unit of measure, but any length unit accepted
86
     * @param Length $z refer to CRS for preferred unit of measure, but any length unit accepted
87
     */
88 522
    public static function create(Geocentric $crs, Length $x, Length $y, Length $z, ?DateTimeInterface $epoch = null): self
89
    {
90 522
        return new static($crs, $x, $y, $z, $epoch);
91
    }
92
93 297
    public function getX(): Length
94
    {
95 297
        return $this->x;
96
    }
97
98 297
    public function getY(): Length
99
    {
100 297
        return $this->y;
101
    }
102
103 297
    public function getZ(): Length
104
    {
105 297
        return $this->z;
106
    }
107
108 63
    public function getCRS(): Geocentric
109
    {
110 63
        return $this->crs;
111
    }
112
113 54
    public function getCoordinateEpoch(): ?DateTimeImmutable
114
    {
115 54
        return $this->epoch;
116
    }
117
118
    /**
119
     * Calculate surface distance between two points.
120
     */
121 18
    public function calculateDistance(Point $to): Length
122
    {
123
        try {
124 18
            if ($to instanceof ConvertiblePoint) {
125 18
                $to = $to->convert($this->crs);
126
            }
127
        } finally {
128 18
            if ($to->getCRS()->getSRID() !== $this->crs->getSRID()) {
129 9
                throw new InvalidCoordinateReferenceSystemException('Can only calculate distances between two points in the same CRS');
130
            }
131
132
            /* @var GeocentricPoint $to */
133 9
            return static::vincenty($this->asGeographicValue(), $to->asGeographicValue(), $this->getCRS()->getDatum()->getEllipsoid());
134
        }
135
    }
136
137 27
    public function __toString(): string
138
    {
139 27
        return "({$this->x}, {$this->y}, {$this->z})";
140
    }
141
142
    /**
143
     * Geographic/geocentric conversions
144
     * In applications it is often concatenated with the 3- 7- or 10-parameter transformations 9603, 9606, 9607 or
145
     * 9636 to form a geographic to geographic transformation.
146
     */
147 180
    public function geographicGeocentric(
148
        Geographic2D|Geographic3D $to
149
    ): GeographicPoint {
150 180
        $geocentricValue = new GeocentricValue($this->x, $this->y, $this->z, $to->getDatum());
151 180
        $asGeographic = $geocentricValue->asGeographicValue();
152
153 180
        return GeographicPoint::create($to, $asGeographic->getLatitude(), $asGeographic->getLongitude(), $to instanceof Geographic3D ? $asGeographic->getHeight() : null, $this->epoch);
154
    }
155
156
    /**
157
     * Coordinate Frame rotation (geocentric domain)
158
     * This method is a specific case of the Molodensky-Badekas (CF) method (code 1034) in which the evaluation point
159
     * is at the geocentre with coordinate values of zero. Note the analogy with the Position Vector method (code 1033)
160
     * but beware of the differences!
161
     */
162 63
    public function coordinateFrameRotation(
163
        Geocentric $to,
164
        Length $xAxisTranslation,
165
        Length $yAxisTranslation,
166
        Length $zAxisTranslation,
167
        Angle $xAxisRotation,
168
        Angle $yAxisRotation,
169
        Angle $zAxisRotation,
170
        Scale $scaleDifference
171
    ): self {
172 63
        return $this->coordinateFrameMolodenskyBadekas(
173
            $to,
174
            $xAxisTranslation,
175
            $yAxisTranslation,
176
            $zAxisTranslation,
177
            $xAxisRotation,
178
            $yAxisRotation,
179
            $zAxisRotation,
180
            $scaleDifference,
181 63
            new Metre(0),
182 63
            new Metre(0),
183 63
            new Metre(0)
184
        );
185
    }
186
187
    /**
188
     * Molodensky-Badekas (CF geocentric domain)
189
     * See method codes 1039 and 9636 for this operation in other coordinate domains and method code 1061 for opposite
190
     * rotation convention in geocentric domain.
191
     */
192 72
    public function coordinateFrameMolodenskyBadekas(
193
        Geocentric $to,
194
        Length $xAxisTranslation,
195
        Length $yAxisTranslation,
196
        Length $zAxisTranslation,
197
        Angle $xAxisRotation,
198
        Angle $yAxisRotation,
199
        Angle $zAxisRotation,
200
        Scale $scaleDifference,
201
        Length $ordinate1OfEvaluationPoint,
202
        Length $ordinate2OfEvaluationPoint,
203
        Length $ordinate3OfEvaluationPoint
204
    ): self {
205 72
        $xs = $this->x->asMetres()->getValue();
206 72
        $ys = $this->y->asMetres()->getValue();
207 72
        $zs = $this->z->asMetres()->getValue();
208 72
        $tx = $xAxisTranslation->asMetres()->getValue();
209 72
        $ty = $yAxisTranslation->asMetres()->getValue();
210 72
        $tz = $zAxisTranslation->asMetres()->getValue();
211 72
        $rx = $xAxisRotation->asRadians()->getValue();
212 72
        $ry = $yAxisRotation->asRadians()->getValue();
213 72
        $rz = $zAxisRotation->asRadians()->getValue();
214 72
        $M = 1 + $scaleDifference->asUnity()->getValue();
215 72
        $xp = $ordinate1OfEvaluationPoint->asMetres()->getValue();
216 72
        $yp = $ordinate2OfEvaluationPoint->asMetres()->getValue();
217 72
        $zp = $ordinate3OfEvaluationPoint->asMetres()->getValue();
218
219 72
        $xt = $M * ((($xs - $xp) * 1) + (($ys - $yp) * $rz) + (($zs - $zp) * -$ry)) + $tx + $xp;
220 72
        $yt = $M * ((($xs - $xp) * -$rz) + (($ys - $yp) * 1) + (($zs - $zp) * $rx)) + $ty + $yp;
221 72
        $zt = $M * ((($xs - $xp) * $ry) + (($ys - $yp) * -$rx) + (($zs - $zp) * 1)) + $tz + $zp;
222
223 72
        return static::create($to, new Metre($xt), new Metre($yt), new Metre($zt), $this->epoch);
224
    }
225
226
    /**
227
     * Position Vector transformation (geocentric domain)
228
     * This method is a specific case of the Molodensky-Badekas (PV) method (code 1061) in which the evaluation point
229
     * is the geocentre with coordinate values of zero. Note the analogy with the Coordinate Frame method (code 1032)
230
     * but beware of the differences!
231
     */
232 72
    public function positionVectorTransformation(
233
        Geocentric $to,
234
        Length $xAxisTranslation,
235
        Length $yAxisTranslation,
236
        Length $zAxisTranslation,
237
        Angle $xAxisRotation,
238
        Angle $yAxisRotation,
239
        Angle $zAxisRotation,
240
        Scale $scaleDifference
241
    ): self {
242 72
        return $this->positionVectorMolodenskyBadekas(
243
            $to,
244
            $xAxisTranslation,
245
            $yAxisTranslation,
246
            $zAxisTranslation,
247
            $xAxisRotation,
248
            $yAxisRotation,
249
            $zAxisRotation,
250
            $scaleDifference,
251 72
            new Metre(0),
252 72
            new Metre(0),
253 72
            new Metre(0)
254
        );
255
    }
256
257
    /**
258
     * Molodensky-Badekas (PV geocentric domain)
259
     * See method codes 1062 and 1063 for this operation in other coordinate domains and method code 1034 for opposite
260
     * rotation convention in geocentric domain.
261
     */
262 81
    public function positionVectorMolodenskyBadekas(
263
        Geocentric $to,
264
        Length $xAxisTranslation,
265
        Length $yAxisTranslation,
266
        Length $zAxisTranslation,
267
        Angle $xAxisRotation,
268
        Angle $yAxisRotation,
269
        Angle $zAxisRotation,
270
        Scale $scaleDifference,
271
        Length $ordinate1OfEvaluationPoint,
272
        Length $ordinate2OfEvaluationPoint,
273
        Length $ordinate3OfEvaluationPoint
274
    ): self {
275 81
        $xs = $this->x->asMetres()->getValue();
276 81
        $ys = $this->y->asMetres()->getValue();
277 81
        $zs = $this->z->asMetres()->getValue();
278 81
        $tx = $xAxisTranslation->asMetres()->getValue();
279 81
        $ty = $yAxisTranslation->asMetres()->getValue();
280 81
        $tz = $zAxisTranslation->asMetres()->getValue();
281 81
        $rx = $xAxisRotation->asRadians()->getValue();
282 81
        $ry = $yAxisRotation->asRadians()->getValue();
283 81
        $rz = $zAxisRotation->asRadians()->getValue();
284 81
        $M = 1 + $scaleDifference->asUnity()->getValue();
285 81
        $xp = $ordinate1OfEvaluationPoint->asMetres()->getValue();
286 81
        $yp = $ordinate2OfEvaluationPoint->asMetres()->getValue();
287 81
        $zp = $ordinate3OfEvaluationPoint->asMetres()->getValue();
288
289 81
        $xt = $M * ((($xs - $xp) * 1) + (($ys - $yp) * -$rz) + (($zs - $zp) * $ry)) + $tx + $xp;
290 81
        $yt = $M * ((($xs - $xp) * $rz) + (($ys - $yp) * 1) + (($zs - $zp) * -$rx)) + $ty + $yp;
291 81
        $zt = $M * ((($xs - $xp) * -$ry) + (($ys - $yp) * $rx) + (($zs - $zp) * 1)) + $tz + $zp;
292
293 81
        return static::create($to, new Metre($xt), new Metre($yt), new Metre($zt), $this->epoch);
294
    }
295
296
    /**
297
     * Geocentric translations
298
     * This method allows calculation of geocentric coords in the target system by adding the parameter values to the
299
     * corresponding coordinates of the point in the source system. See methods 1035 and 9603 for similar tfms
300
     * operating between other CRSs types.
301
     */
302 9
    public function geocentricTranslation(
303
        Geocentric $to,
304
        Length $xAxisTranslation,
305
        Length $yAxisTranslation,
306
        Length $zAxisTranslation
307
    ): self {
308 9
        return $this->positionVectorTransformation(
309
            $to,
310
            $xAxisTranslation,
311
            $yAxisTranslation,
312
            $zAxisTranslation,
313 9
            new Radian(0),
314 9
            new Radian(0),
315 9
            new Radian(0),
316 9
            new Unity(0)
317
        );
318
    }
319
320
    /**
321
     * Time-dependent Coordinate Frame rotation (geocen)
322
     * Note the analogy with the Time-dependent Position Vector transformation (code 1053) but beware of the
323
     * differences!  The Position Vector convention is used by IAG. See method codes 1057 and 1058 for similar methods
324
     * operating between other CRS types.
325
     */
326 18
    public function timeDependentCoordinateFrameRotation(
327
        Geocentric $to,
328
        Length $xAxisTranslation,
329
        Length $yAxisTranslation,
330
        Length $zAxisTranslation,
331
        Angle $xAxisRotation,
332
        Angle $yAxisRotation,
333
        Angle $zAxisRotation,
334
        Scale $scaleDifference,
335
        Rate $rateOfChangeOfXAxisTranslation,
336
        Rate $rateOfChangeOfYAxisTranslation,
337
        Rate $rateOfChangeOfZAxisTranslation,
338
        Rate $rateOfChangeOfXAxisRotation,
339
        Rate $rateOfChangeOfYAxisRotation,
340
        Rate $rateOfChangeOfZAxisRotation,
341
        Rate $rateOfChangeOfScaleDifference,
342
        Time $parameterReferenceEpoch
343
    ): self {
344 18
        if ($this->epoch === null) {
345 9
            throw new InvalidCoordinateException('This transformation requires an epoch, none given');
346
        }
347
348
        // Points use PHP DateTimes for ease of use, but transformations use decimal years...
349 9
        $pointEpoch = Year::fromDateTime($this->epoch);
350 9
        $yearsToAdjust = $pointEpoch->subtract($parameterReferenceEpoch)->getValue();
351 9
        $xAxisTranslation = $xAxisTranslation->add($rateOfChangeOfXAxisTranslation->getChangePerYear()->multiply($yearsToAdjust));
0 ignored issues
show
Bug introduced by
The method multiply() does not exist on PHPCoord\UnitOfMeasure\UnitOfMeasure. Since it exists in all sub-types, consider adding an abstract or default implementation to PHPCoord\UnitOfMeasure\UnitOfMeasure. ( Ignorable by Annotation )

If this is a false-positive, you can also ignore this issue in your code via the ignore-call  annotation

351
        $xAxisTranslation = $xAxisTranslation->add($rateOfChangeOfXAxisTranslation->getChangePerYear()->/** @scrutinizer ignore-call */ multiply($yearsToAdjust));
Loading history...
352 9
        $yAxisTranslation = $yAxisTranslation->add($rateOfChangeOfYAxisTranslation->getChangePerYear()->multiply($yearsToAdjust));
353 9
        $zAxisTranslation = $zAxisTranslation->add($rateOfChangeOfZAxisTranslation->getChangePerYear()->multiply($yearsToAdjust));
354 9
        $xAxisRotation = $xAxisRotation->add($rateOfChangeOfXAxisRotation->getChangePerYear()->multiply($yearsToAdjust));
355 9
        $yAxisRotation = $yAxisRotation->add($rateOfChangeOfYAxisRotation->getChangePerYear()->multiply($yearsToAdjust));
356 9
        $zAxisRotation = $zAxisRotation->add($rateOfChangeOfZAxisRotation->getChangePerYear()->multiply($yearsToAdjust));
357 9
        $scaleDifference = $scaleDifference->add($rateOfChangeOfScaleDifference->getChangePerYear()->multiply($yearsToAdjust));
358
359 9
        return $this->coordinateFrameRotation($to, $xAxisTranslation, $yAxisTranslation, $zAxisTranslation, $xAxisRotation, $yAxisRotation, $zAxisRotation, $scaleDifference);
360
    }
361
362
    /**
363
     * Time-dependent Position Vector tfm (geocentric)
364
     * Note the analogy with the Time-dependent Coordinate Frame rotation (code 1056) but beware of the differences!
365
     * The Position Vector convention is used by IAG. See method codes 1054 and 1055 for similar methods operating
366
     * between other CRS types.
367
     */
368 36
    public function timeDependentPositionVectorTransformation(
369
        Geocentric $to,
370
        Length $xAxisTranslation,
371
        Length $yAxisTranslation,
372
        Length $zAxisTranslation,
373
        Angle $xAxisRotation,
374
        Angle $yAxisRotation,
375
        Angle $zAxisRotation,
376
        Scale $scaleDifference,
377
        Rate $rateOfChangeOfXAxisTranslation,
378
        Rate $rateOfChangeOfYAxisTranslation,
379
        Rate $rateOfChangeOfZAxisTranslation,
380
        Rate $rateOfChangeOfXAxisRotation,
381
        Rate $rateOfChangeOfYAxisRotation,
382
        Rate $rateOfChangeOfZAxisRotation,
383
        Rate $rateOfChangeOfScaleDifference,
384
        Time $parameterReferenceEpoch
385
    ): self {
386 36
        if ($this->epoch === null) {
387 9
            throw new InvalidCoordinateException('This transformation requires an epoch, none given');
388
        }
389
390
        // Points use PHP DateTimes for ease of use, but transformations use decimal years...
391 27
        $pointEpoch = Year::fromDateTime($this->epoch);
392 27
        $yearsToAdjust = $pointEpoch->subtract($parameterReferenceEpoch)->getValue();
393 27
        $xAxisTranslation = $xAxisTranslation->add($rateOfChangeOfXAxisTranslation->getChangePerYear()->multiply($yearsToAdjust));
394 27
        $yAxisTranslation = $yAxisTranslation->add($rateOfChangeOfYAxisTranslation->getChangePerYear()->multiply($yearsToAdjust));
395 27
        $zAxisTranslation = $zAxisTranslation->add($rateOfChangeOfZAxisTranslation->getChangePerYear()->multiply($yearsToAdjust));
396 27
        $xAxisRotation = $xAxisRotation->add($rateOfChangeOfXAxisRotation->getChangePerYear()->multiply($yearsToAdjust));
397 27
        $yAxisRotation = $yAxisRotation->add($rateOfChangeOfYAxisRotation->getChangePerYear()->multiply($yearsToAdjust));
398 27
        $zAxisRotation = $zAxisRotation->add($rateOfChangeOfZAxisRotation->getChangePerYear()->multiply($yearsToAdjust));
399 27
        $scaleDifference = $scaleDifference->add($rateOfChangeOfScaleDifference->getChangePerYear()->multiply($yearsToAdjust));
400
401 27
        return $this->positionVectorTransformation($to, $xAxisTranslation, $yAxisTranslation, $zAxisTranslation, $xAxisRotation, $yAxisRotation, $zAxisRotation, $scaleDifference);
402
    }
403
404
    /**
405
     * Time-specific Coordinate Frame rotation (geocen)
406
     * Note the analogy with the Time-specific Position Vector method (code 1065) but beware of the differences!
407
     */
408 36
    public function timeSpecificCoordinateFrameRotation(
409
        Geocentric $to,
410
        Length $xAxisTranslation,
411
        Length $yAxisTranslation,
412
        Length $zAxisTranslation,
413
        Angle $xAxisRotation,
414
        Angle $yAxisRotation,
415
        Angle $zAxisRotation,
416
        Scale $scaleDifference,
417
        Time $transformationReferenceEpoch
418
    ): self {
419 36
        if ($this->epoch === null) {
420 9
            throw new InvalidCoordinateException(sprintf('This transformation is only valid for epoch %s, none given', $transformationReferenceEpoch->getValue()));
421
        }
422
423
        // Points use PHP DateTimes for ease of use, but transformations use decimal years...
424 27
        $pointEpoch = Year::fromDateTime($this->epoch);
425
426 27
        if (abs($pointEpoch->getValue() - $transformationReferenceEpoch->getValue()) > 0.001) {
427 9
            throw new InvalidCoordinateException(sprintf('This transformation is only valid for epoch %s, got %s', $transformationReferenceEpoch, $pointEpoch));
428
        }
429
430 18
        return $this->coordinateFrameRotation($to, $xAxisTranslation, $yAxisTranslation, $zAxisTranslation, $xAxisRotation, $yAxisRotation, $zAxisRotation, $scaleDifference);
431
    }
432
433
    /**
434
     * Time-specific Position Vector transform (geocen)
435
     * Note the analogy with the Time-specifc Coordinate Frame method (code 1066) but beware of the differences!
436
     */
437 27
    public function timeSpecificPositionVectorTransformation(
438
        Geocentric $to,
439
        Length $xAxisTranslation,
440
        Length $yAxisTranslation,
441
        Length $zAxisTranslation,
442
        Angle $xAxisRotation,
443
        Angle $yAxisRotation,
444
        Angle $zAxisRotation,
445
        Scale $scaleDifference,
446
        Time $transformationReferenceEpoch
447
    ): self {
448 27
        if ($this->epoch === null) {
449 9
            throw new InvalidCoordinateException(sprintf('This transformation is only valid for epoch %s, none given', $transformationReferenceEpoch->getValue()));
450
        }
451
452
        // Points use PHP DateTimes for ease of use, but transformations use decimal years...
453 18
        $pointEpoch = Year::fromDateTime($this->epoch);
454
455 18
        if (abs($pointEpoch->getValue() - $transformationReferenceEpoch->getValue()) > 0.001) {
456 9
            throw new InvalidCoordinateException(sprintf('This transformation is only valid for epoch %s, got %s', $transformationReferenceEpoch, $pointEpoch));
457
        }
458
459 9
        return $this->positionVectorTransformation($to, $xAxisTranslation, $yAxisTranslation, $zAxisTranslation, $xAxisRotation, $yAxisRotation, $zAxisRotation, $scaleDifference);
460
    }
461
462 18
    public function asGeographicValue(): GeographicValue
463
    {
464 18
        return (new GeocentricValue($this->x, $this->y, $this->z, $this->getCRS()->getDatum()))->asGeographicValue();
465
    }
466
}
467