Passed
Push — master ( 705a37...2b8bad )
by Doug
17:25 queued 45s
created

GeocentricPoint::getCoordinateEpoch()   A

Complexity

Conditions 1
Paths 1

Size

Total Lines 3
Code Lines 1

Duplication

Lines 0
Ratio 0 %

Code Coverage

Tests 2
CRAP Score 1

Importance

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

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

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