Passed
Push — 4.0.x ( 414ad0...5c8c8d )
by Doug
05:04 queued 12s
created

GeocentricPoint::calculateDistance()   A

Complexity

Conditions 3
Paths 6

Size

Total Lines 13
Code Lines 7

Duplication

Lines 0
Ratio 0 %

Code Coverage

Tests 6
CRAP Score 3.0261

Importance

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

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