1
|
|
|
<?php |
2
|
|
|
|
3
|
|
|
/** |
4
|
|
|
* PHPCoord. |
5
|
|
|
* |
6
|
|
|
* @author Doug Wright |
7
|
|
|
*/ |
8
|
|
|
declare(strict_types=1); |
9
|
|
|
|
10
|
|
|
namespace PHPCoord\Point; |
11
|
|
|
|
12
|
|
|
use DateTime; |
13
|
|
|
use DateTimeImmutable; |
14
|
|
|
use DateTimeInterface; |
15
|
|
|
use PHPCoord\CoordinateOperation\AutoConversion; |
16
|
|
|
use PHPCoord\CoordinateOperation\ConvertiblePoint; |
17
|
|
|
use PHPCoord\CoordinateOperation\GeographicGeoidHeightGrid; |
18
|
|
|
use PHPCoord\CoordinateOperation\OSTNOSGM15Grid; |
19
|
|
|
use PHPCoord\CoordinateReferenceSystem\Compound; |
20
|
|
|
use PHPCoord\CoordinateReferenceSystem\CoordinateReferenceSystem; |
21
|
|
|
use PHPCoord\CoordinateReferenceSystem\Geocentric; |
22
|
|
|
use PHPCoord\CoordinateReferenceSystem\Geographic2D; |
23
|
|
|
use PHPCoord\CoordinateReferenceSystem\Geographic3D; |
24
|
|
|
use PHPCoord\CoordinateReferenceSystem\Projected; |
25
|
|
|
use PHPCoord\CoordinateReferenceSystem\Vertical; |
26
|
|
|
use PHPCoord\CoordinateSystem\Cartesian; |
27
|
|
|
use PHPCoord\Datum\Datum; |
28
|
|
|
use PHPCoord\Exception\InvalidCoordinateReferenceSystemException; |
29
|
|
|
use PHPCoord\Exception\UnknownConversionException; |
30
|
|
|
use PHPCoord\UnitOfMeasure\Angle\Angle; |
31
|
|
|
use PHPCoord\UnitOfMeasure\Angle\Degree; |
32
|
|
|
use PHPCoord\UnitOfMeasure\Length\Length; |
33
|
|
|
use PHPCoord\UnitOfMeasure\Length\Metre; |
34
|
|
|
use PHPCoord\UnitOfMeasure\Scale\Unity; |
35
|
|
|
use Throwable; |
36
|
|
|
|
37
|
|
|
use function assert; |
38
|
|
|
|
39
|
|
|
/** |
40
|
|
|
* Coordinate representing a point expressed in 2 different CRSs (2D horizontal + 1D Vertical). |
41
|
|
|
*/ |
42
|
|
|
class CompoundPoint extends Point implements ConvertiblePoint |
43
|
|
|
{ |
44
|
|
|
use AutoConversion { |
45
|
|
|
convert as protected autoConvert; |
46
|
|
|
} |
47
|
|
|
|
48
|
|
|
/** |
49
|
|
|
* Horizontal point. |
50
|
|
|
*/ |
51
|
|
|
protected GeographicPoint|ProjectedPoint $horizontalPoint; |
52
|
|
|
|
53
|
|
|
/** |
54
|
|
|
* Vertical point. |
55
|
|
|
*/ |
56
|
|
|
protected VerticalPoint $verticalPoint; |
57
|
|
|
|
58
|
|
|
/** |
59
|
|
|
* Coordinate reference system. |
60
|
|
|
*/ |
61
|
|
|
protected Compound $crs; |
62
|
|
|
|
63
|
|
|
/** |
64
|
|
|
* Coordinate epoch (date for which the specified coordinates represented this point). |
65
|
|
|
*/ |
66
|
|
|
protected ?DateTimeImmutable $epoch; |
67
|
81 |
|
|
68
|
|
|
protected function __construct(Compound $crs, GeographicPoint|ProjectedPoint $horizontalPoint, VerticalPoint $verticalPoint, ?DateTimeInterface $epoch = null) |
69
|
81 |
|
{ |
70
|
81 |
|
$this->horizontalPoint = $horizontalPoint; |
71
|
81 |
|
$this->verticalPoint = $verticalPoint; |
72
|
|
|
$this->crs = $crs; |
73
|
81 |
|
|
74
|
9 |
|
if ($epoch instanceof DateTime) { |
75
|
|
|
$epoch = DateTimeImmutable::createFromMutable($epoch); |
76
|
81 |
|
} |
77
|
|
|
$this->epoch = $epoch; |
78
|
|
|
} |
79
|
81 |
|
|
80
|
|
|
public static function create(Compound $crs, GeographicPoint|ProjectedPoint $horizontalPoint, VerticalPoint $verticalPoint, ?DateTimeInterface $epoch = null): self |
81
|
81 |
|
{ |
82
|
|
|
return new self($crs, $horizontalPoint, $verticalPoint, $epoch); |
83
|
|
|
} |
84
|
79 |
|
|
85
|
|
|
public function getHorizontalPoint(): GeographicPoint|ProjectedPoint |
86
|
79 |
|
{ |
87
|
|
|
return $this->horizontalPoint; |
88
|
|
|
} |
89
|
42 |
|
|
90
|
|
|
public function getVerticalPoint(): VerticalPoint |
91
|
42 |
|
{ |
92
|
|
|
return $this->verticalPoint; |
93
|
|
|
} |
94
|
73 |
|
|
95
|
|
|
public function getCRS(): Compound |
96
|
73 |
|
{ |
97
|
|
|
return $this->crs; |
98
|
|
|
} |
99
|
30 |
|
|
100
|
|
|
public function getCoordinateEpoch(): ?DateTimeImmutable |
101
|
30 |
|
{ |
102
|
|
|
return $this->epoch; |
103
|
|
|
} |
104
|
|
|
|
105
|
|
|
/** |
106
|
|
|
* Calculate distance between two points. |
107
|
18 |
|
*/ |
108
|
|
|
public function calculateDistance(Point $to): Length |
109
|
|
|
{ |
110
|
18 |
|
try { |
111
|
18 |
|
if ($to instanceof ConvertiblePoint) { |
112
|
|
|
$to = $to->convert($this->horizontalPoint->getCRS()); |
113
|
|
|
} |
114
|
18 |
|
} finally { |
115
|
9 |
|
if ($to->getCRS()->getSRID() !== $this->horizontalPoint->getCRS()->getSRID()) { |
116
|
|
|
throw new InvalidCoordinateReferenceSystemException('Can only calculate distances between two points in the same CRS'); |
117
|
|
|
} |
118
|
|
|
|
119
|
9 |
|
/* @var CompoundPoint $to */ |
120
|
|
|
return $this->horizontalPoint->calculateDistance($to); |
121
|
|
|
} |
122
|
|
|
} |
123
|
46 |
|
|
124
|
|
|
public function convert(Compound|Geocentric|Geographic2D|Geographic3D|Projected|Vertical $to, bool $ignoreBoundaryRestrictions = false): Point |
125
|
|
|
{ |
126
|
46 |
|
try { |
127
|
27 |
|
return $this->autoConvert($to, $ignoreBoundaryRestrictions); |
128
|
|
|
} catch (UnknownConversionException $e) { |
129
|
27 |
|
// if 2D target, try again with just the horizontal component |
130
|
18 |
|
if ($to instanceof Geographic2D || $to instanceof Projected) { |
131
|
|
|
return $this->horizontalPoint->convert($to, $ignoreBoundaryRestrictions); |
132
|
|
|
} |
133
|
|
|
|
134
|
9 |
|
// try separate horizontal + vertical conversions and stitch results together |
135
|
|
|
if ($to instanceof Compound) { |
136
|
9 |
|
/** @var GeographicPoint|ProjectedPoint $newHorizontalPoint */ |
137
|
|
|
$newHorizontalPoint = $this->horizontalPoint->convert($to->getHorizontal()); |
138
|
9 |
|
|
139
|
9 |
|
if ($this->getCRS()->getVertical()->getSRID() !== $to->getVertical()->getSRID()) { |
140
|
|
|
$path = $this->findOperationPath($this->getCRS()->getVertical(), $to->getVertical(), $ignoreBoundaryRestrictions); |
141
|
9 |
|
|
142
|
9 |
|
if ($path) { |
|
|
|
|
143
|
9 |
|
$newVerticalPoint = $this->verticalPoint; |
144
|
9 |
|
foreach ($path as $step) { |
145
|
|
|
$target = CoordinateReferenceSystem::fromSRID($step['in_reverse'] ? $step['source_crs'] : $step['target_crs']); |
146
|
9 |
|
/** @var VerticalPoint $newVerticalPoint */ |
147
|
|
|
$newVerticalPoint = $newVerticalPoint->performOperation($step['operation'], $target, $step['in_reverse'], ['horizontalPoint' => $newHorizontalPoint]); |
148
|
|
|
} |
149
|
9 |
|
|
150
|
|
|
return static::create($to, $newHorizontalPoint, $newVerticalPoint, $this->epoch); |
151
|
|
|
} |
152
|
|
|
} |
153
|
|
|
} |
154
|
|
|
// try converting to any/all of the other Compound CRSs that include the same vertical CRS, where the |
155
|
|
|
// horizontal CRS has a 3D equivalent. From there, try converting using the usual mechanisms |
156
|
|
|
$candidateIntermediates = Compound::findFromVertical($this->verticalPoint->getCRS()); |
157
|
|
|
unset($candidateIntermediates[$this->getCRS()->getSRID()]); |
158
|
|
|
|
159
|
|
|
foreach ($candidateIntermediates as $candidateIntermediate) { |
160
|
|
|
try { |
161
|
|
|
if ($candidateIntermediate->getHorizontal() instanceof Geographic2D && $candidateIntermediate->getHorizontal()->getBaseCRS() instanceof Geographic3D) { |
162
|
|
|
/** @var GeographicPoint|ProjectedPoint $candidateHorizontalPoint */ |
163
|
|
|
$candidateHorizontalPoint = $this->horizontalPoint->convert($candidateIntermediate->getHorizontal()); |
164
|
|
|
$candidateIntermediatePoint = self::create( |
165
|
|
|
$candidateIntermediate, |
166
|
|
|
$candidateHorizontalPoint, |
167
|
|
|
$this->verticalPoint, |
168
|
|
|
$this->epoch, |
169
|
|
|
); |
170
|
|
|
|
171
|
|
|
/** @var GeographicPoint $candidateBaseCRSPoint */ |
172
|
|
|
$candidateBaseCRSPoint = $candidateIntermediatePoint->convert($candidateIntermediate->getHorizontal()->getBaseCRS()); |
173
|
|
|
|
174
|
|
|
return $candidateBaseCRSPoint->convert($to); |
175
|
|
|
} |
176
|
|
|
} catch (Throwable) { |
|
|
|
|
177
|
|
|
} |
178
|
|
|
} |
179
|
|
|
throw $e; |
180
|
|
|
} |
181
|
|
|
} |
182
|
27 |
|
|
183
|
|
|
public function __toString(): string |
184
|
27 |
|
{ |
185
|
|
|
return "({$this->horizontalPoint}, {$this->verticalPoint})"; |
186
|
|
|
} |
187
|
|
|
|
188
|
|
|
/** |
189
|
|
|
* Geographic2D with Height Offsets. |
190
|
|
|
* This transformation allows calculation of coordinates in the target system by adding the parameter value to the |
191
|
|
|
* coordinate values of the point in the source system. |
192
|
18 |
|
*/ |
193
|
|
|
public function geographic2DWithHeightOffsets( |
194
|
|
|
Geographic3D $to, |
195
|
|
|
Angle $latitudeOffset, |
196
|
|
|
Angle $longitudeOffset, |
197
|
|
|
Length $geoidUndulation |
198
|
18 |
|
): GeographicPoint { |
199
|
18 |
|
assert($this->horizontalPoint instanceof GeographicPoint); |
200
|
18 |
|
$horizontalPoint = $this->horizontalPoint; |
201
|
18 |
|
$toLatitude = $horizontalPoint->getLatitude()->add($latitudeOffset); |
202
|
18 |
|
$toLongitude = $horizontalPoint->getLongitude()->add($longitudeOffset); |
203
|
|
|
$toHeight = $this->verticalPoint->getHeight()->add($geoidUndulation); |
204
|
18 |
|
|
205
|
|
|
return GeographicPoint::create($to, $toLatitude, $toLongitude, $toHeight, $this->epoch); |
206
|
|
|
} |
207
|
|
|
|
208
|
|
|
/** |
209
|
|
|
* Geog3D to Geog2D+GravityRelatedHeight (OSGM-GB). |
210
|
|
|
* Uses ETRS89 / National Grid as an intermediate coordinate system for bi-linear interpolation of gridded grid |
211
|
|
|
* coordinate differences. |
212
|
|
|
*/ |
213
|
|
|
public function geographic3DTo2DPlusGravityHeightOSGM15( |
214
|
|
|
Geographic3D $to, |
215
|
|
|
OSTNOSGM15Grid $geoidHeightCorrectionModelFile |
216
|
|
|
): GeographicPoint { |
217
|
|
|
assert($this->horizontalPoint instanceof GeographicPoint); |
218
|
|
|
$osgb36NationalGrid = Projected::fromSRID(Projected::EPSG_OSGB36_BRITISH_NATIONAL_GRID); |
219
|
|
|
$etrs89NationalGrid = new Projected( |
220
|
|
|
'ETRS89 / National Grid', |
221
|
|
|
Cartesian::fromSRID(Cartesian::EPSG_2D_AXES_EASTING_NORTHING_E_N_ORIENTATIONS_EAST_NORTH_UOM_M), |
222
|
|
|
Datum::fromSRID(Datum::EPSG_EUROPEAN_TERRESTRIAL_REFERENCE_SYSTEM_1989_ENSEMBLE), |
223
|
|
|
$osgb36NationalGrid->getBoundingArea() |
224
|
|
|
); |
225
|
|
|
|
226
|
|
|
$projected = $this->horizontalPoint->transverseMercator($etrs89NationalGrid, new Degree(49), new Degree(-2), new Unity(0.9996012717), new Metre(400000), new Metre(-100000)); |
227
|
|
|
|
228
|
|
|
return GeographicPoint::create( |
229
|
|
|
$to, |
230
|
|
|
$this->horizontalPoint->getLatitude(), |
231
|
|
|
$this->horizontalPoint->getLongitude(), |
232
|
|
|
$this->verticalPoint->getHeight()->add($geoidHeightCorrectionModelFile->getHeightAdjustment($projected)), |
233
|
|
|
$this->getCoordinateEpoch() |
234
|
|
|
); |
235
|
|
|
} |
236
|
|
|
|
237
|
|
|
/** |
238
|
|
|
* Geog3D to Geog2D+GravityRelatedHeight. |
239
|
3 |
|
*/ |
240
|
|
|
public function geographic3DTo2DPlusGravityHeightFromGrid( |
241
|
|
|
Geographic3D $to, |
242
|
|
|
GeographicGeoidHeightGrid $geoidHeightCorrectionModelFile |
243
|
3 |
|
): GeographicPoint { |
244
|
|
|
assert($this->horizontalPoint instanceof GeographicPoint); |
245
|
3 |
|
|
246
|
3 |
|
return GeographicPoint::create( |
247
|
3 |
|
$to, |
248
|
3 |
|
$this->horizontalPoint->getLatitude(), |
249
|
3 |
|
$this->horizontalPoint->getLongitude(), |
250
|
3 |
|
$this->verticalPoint->getHeight()->add($geoidHeightCorrectionModelFile->getHeightAdjustment($this->horizontalPoint)), |
251
|
3 |
|
$this->getCoordinateEpoch() |
252
|
|
|
); |
253
|
|
|
} |
254
|
|
|
} |
255
|
|
|
|
This check marks implicit conversions of arrays to boolean values in a comparison. While in PHP an empty array is considered to be equal (but not identical) to false, this is not always apparent.
Consider making the comparison explicit by using
empty(..)
or! empty(...)
instead.