|
1
|
|
|
<?php |
|
2
|
|
|
/** |
|
3
|
|
|
* PHPCoord. |
|
4
|
|
|
* |
|
5
|
|
|
* @author Doug Wright |
|
6
|
|
|
*/ |
|
7
|
|
|
declare(strict_types=1); |
|
8
|
|
|
|
|
9
|
|
|
namespace PHPCoord\CoordinateOperation; |
|
10
|
|
|
|
|
11
|
|
|
use function abs; |
|
12
|
|
|
use function assert; |
|
13
|
|
|
use InvalidArgumentException; |
|
14
|
|
|
use PHPCoord\CoordinateReferenceSystem\Geographic; |
|
15
|
|
|
use PHPCoord\GeographicPoint; |
|
16
|
|
|
use PHPCoord\UnitOfMeasure\Angle\Angle; |
|
17
|
|
|
use PHPCoord\UnitOfMeasure\Angle\ArcSecond; |
|
18
|
|
|
use function round; |
|
19
|
|
|
use SplFileObject; |
|
20
|
|
|
use function unpack; |
|
21
|
|
|
use function usort; |
|
22
|
|
|
|
|
23
|
|
|
class NTv2Grid extends SplFileObject |
|
24
|
|
|
{ |
|
25
|
|
|
private const RECORD_SIZE = 16; |
|
26
|
|
|
private const ITERATION_CONVERGENCE = 0.0001; |
|
27
|
|
|
private const FLAG_WITHIN_LIMITS = 1; |
|
28
|
|
|
private const FLAG_ON_UPPER_LATITUDE = 2; |
|
29
|
|
|
private const FLAG_ON_UPPER_LONGITUDE = 3; |
|
30
|
|
|
private const FLAG_ON_UPPER_LATITUDE_AND_LONGITUDE = 4; |
|
31
|
|
|
|
|
32
|
|
|
private string $integerFormatChar = 'V'; |
|
33
|
|
|
private string $doubleFormatChar = 'e'; |
|
34
|
|
|
private string $floatFormatChar = 'g'; |
|
35
|
|
|
|
|
36
|
|
|
private array $subFileMetaData = []; |
|
37
|
|
|
|
|
38
|
|
|
private int $numberOfOverviewHeaderRecords; |
|
|
|
|
|
|
39
|
|
|
private int $numberOfSubFileHeaderRecords; |
|
|
|
|
|
|
40
|
|
|
private int $numberOfGridShiftSubFiles; |
|
|
|
|
|
|
41
|
|
|
|
|
42
|
|
|
private float $lowerLatitudeLimit; |
|
|
|
|
|
|
43
|
|
|
private float $upperLatitudeLimit; |
|
|
|
|
|
|
44
|
|
|
private float $lowerLongitudeLimit; |
|
|
|
|
|
|
45
|
|
|
private float $upperLongitudeLimit; |
|
|
|
|
|
|
46
|
|
|
private float $latitudeGridInterval; |
|
|
|
|
|
|
47
|
|
|
private float $longitudeGridInterval; |
|
|
|
|
|
|
48
|
|
|
private string $gridShiftDataType; |
|
|
|
|
|
|
49
|
|
|
private string $systemFrom; |
|
|
|
|
|
|
50
|
|
|
private string $systemTo; |
|
|
|
|
|
|
51
|
|
|
private float $semiMinorAxisFrom; |
|
|
|
|
|
|
52
|
|
|
private float $semiMinorAxisTo; |
|
|
|
|
|
|
53
|
|
|
|
|
54
|
|
|
public function __construct($filename) |
|
55
|
|
|
{ |
|
56
|
|
|
parent::__construct($filename); |
|
57
|
|
|
|
|
58
|
|
|
$this->readHeader(); |
|
59
|
|
|
} |
|
60
|
|
|
|
|
61
|
|
|
public function applyForwardAdjustment(GeographicPoint $point, Geographic $to): GeographicPoint |
|
62
|
|
|
{ |
|
63
|
|
|
$adjustment = $this->getAdjustment($point->getLatitude(), $point->getLongitude()); |
|
64
|
|
|
|
|
65
|
|
|
$latitude = $point->getLatitude()->add($adjustment[0]); |
|
66
|
|
|
$longitude = $point->getLongitude()->add($adjustment[1]); |
|
67
|
|
|
|
|
68
|
|
|
return GeographicPoint::create($latitude, $longitude, $point->getHeight(), $to, $point->getCoordinateEpoch()); |
|
69
|
|
|
} |
|
70
|
|
|
|
|
71
|
|
|
public function applyReverseAdjustment(GeographicPoint $point, Geographic $to): GeographicPoint |
|
72
|
|
|
{ |
|
73
|
|
|
$adjustment = [new ArcSecond(0), new ArcSecond(0)]; |
|
74
|
|
|
$latitude = $point->getLatitude(); |
|
75
|
|
|
$longitude = $point->getLongitude(); |
|
76
|
|
|
|
|
77
|
|
|
do { |
|
78
|
|
|
$prevAdjustment = $adjustment; |
|
79
|
|
|
$adjustment = $this->getAdjustment($latitude, $longitude); |
|
80
|
|
|
$latitude = $point->getLatitude()->subtract($adjustment[0]); |
|
81
|
|
|
$longitude = $point->getLongitude()->subtract($adjustment[1]); |
|
82
|
|
|
} while (abs($adjustment[0]->subtract($prevAdjustment[0])->getValue()) > self::ITERATION_CONVERGENCE && abs($adjustment[1]->subtract($prevAdjustment[1])->getValue()) > self::ITERATION_CONVERGENCE); |
|
83
|
|
|
|
|
84
|
|
|
return GeographicPoint::create($latitude, $longitude, $point->getHeight(), $to, $point->getCoordinateEpoch()); |
|
85
|
|
|
} |
|
86
|
|
|
|
|
87
|
|
|
/** |
|
88
|
|
|
* @return ArcSecond[] |
|
89
|
|
|
*/ |
|
90
|
|
|
private function getAdjustment(Angle $latitude, Angle $longitude): array |
|
91
|
|
|
{ |
|
92
|
|
|
// NTv2 is longitude positive *west* |
|
93
|
|
|
$longitude = $longitude->multiply(-1); |
|
94
|
|
|
|
|
95
|
|
|
$latitudeAsSeconds = Angle::convert($latitude, Angle::EPSG_ARC_SECOND)->getValue(); |
|
96
|
|
|
$longitudeAsSeconds = Angle::convert($longitude, Angle::EPSG_ARC_SECOND)->getValue(); |
|
97
|
|
|
[$flag, $gridToUse] = $this->determineBestGrid($latitudeAsSeconds, $longitudeAsSeconds); |
|
98
|
|
|
|
|
99
|
|
|
$rowIndex = (int) (($latitudeAsSeconds - $gridToUse['S_LAT']) / $gridToUse['LAT_INC']); |
|
100
|
|
|
$columnIndex = (int) (($longitudeAsSeconds - $gridToUse['E_LONG']) / $gridToUse['LONG_INC']); |
|
101
|
|
|
$gridPointsPerRow = (int) (($gridToUse['W_LONG'] - $gridToUse['E_LONG']) / $gridToUse['LONG_INC']) + 1; |
|
102
|
|
|
$gridPointsPerColumn = (int) (($gridToUse['N_LAT'] - $gridToUse['S_LAT']) / $gridToUse['LAT_INC']) + 1; |
|
103
|
|
|
$numberOfRecords = $gridPointsPerRow * $gridPointsPerColumn; |
|
104
|
|
|
assert($numberOfRecords === $gridToUse['GS_COUNT']); |
|
105
|
|
|
|
|
106
|
|
|
if ($flag === self::FLAG_WITHIN_LIMITS) { |
|
107
|
|
|
$recordIndexLR = $rowIndex * $gridPointsPerRow + $columnIndex; |
|
108
|
|
|
$recordIndexLL = $recordIndexLR + 1; |
|
109
|
|
|
$recordIndexUR = $recordIndexLR + $gridPointsPerRow; |
|
110
|
|
|
$recordIndexUL = $recordIndexUR + 1; |
|
111
|
|
|
} elseif ($flag === self::FLAG_ON_UPPER_LATITUDE) { |
|
112
|
|
|
$recordIndexLR = $rowIndex * $gridPointsPerRow + $columnIndex; |
|
113
|
|
|
$recordIndexLL = $recordIndexLR + 1; |
|
114
|
|
|
$recordIndexUR = $recordIndexLR; |
|
115
|
|
|
$recordIndexUL = $recordIndexUR; |
|
116
|
|
|
} elseif ($flag === self::FLAG_ON_UPPER_LONGITUDE) { |
|
117
|
|
|
$recordIndexLR = $rowIndex * $gridPointsPerRow + $columnIndex; |
|
118
|
|
|
$recordIndexUR = $recordIndexLR + $gridPointsPerRow; |
|
119
|
|
|
$recordIndexLL = $recordIndexLR; |
|
120
|
|
|
$recordIndexUL = $recordIndexUR; |
|
121
|
|
|
} elseif ($flag === self::FLAG_ON_UPPER_LATITUDE_AND_LONGITUDE) { |
|
122
|
|
|
$recordIndexLR = $rowIndex * $gridPointsPerRow + $columnIndex; |
|
123
|
|
|
$recordIndexLL = $recordIndexLR; |
|
124
|
|
|
$recordIndexUR = $recordIndexLR; |
|
125
|
|
|
$recordIndexUL = $recordIndexUR; |
|
126
|
|
|
} |
|
127
|
|
|
|
|
128
|
|
|
$latitudeR = $gridToUse['S_LAT'] + $rowIndex * $gridToUse['LAT_INC']; |
|
129
|
|
|
$longitudeL = $gridToUse['E_LONG'] + $columnIndex * $gridToUse['LONG_INC']; |
|
130
|
|
|
|
|
131
|
|
|
$recordLR = $this->getRecord($gridToUse['offsetStart'] + (11 + $recordIndexLR) * self::RECORD_SIZE); |
|
|
|
|
|
|
132
|
|
|
$recordLL = $this->getRecord($gridToUse['offsetStart'] + (11 + $recordIndexLL) * self::RECORD_SIZE); |
|
|
|
|
|
|
133
|
|
|
$recordUR = $this->getRecord($gridToUse['offsetStart'] + (11 + $recordIndexUR) * self::RECORD_SIZE); |
|
|
|
|
|
|
134
|
|
|
$recordUL = $this->getRecord($gridToUse['offsetStart'] + (11 + $recordIndexUL) * self::RECORD_SIZE); |
|
|
|
|
|
|
135
|
|
|
|
|
136
|
|
|
$x = ($latitudeAsSeconds - $latitudeR) / $gridToUse['LAT_INC']; |
|
137
|
|
|
$y = ($longitudeAsSeconds - $longitudeL) / $gridToUse['LONG_INC']; |
|
138
|
|
|
assert($x >= 0 && $x <= 1); |
|
139
|
|
|
assert($y >= 0 && $y <= 1); |
|
140
|
|
|
|
|
141
|
|
|
$eShiftLatitude = $recordLR['LATITUDE_SHIFT'] + ($recordLL['LATITUDE_SHIFT'] - $recordLR['LATITUDE_SHIFT']) * $y; |
|
142
|
|
|
$fShiftLatitude = $recordUR['LATITUDE_SHIFT'] + ($recordUL['LATITUDE_SHIFT'] - $recordUR['LATITUDE_SHIFT']) * $y; |
|
143
|
|
|
$pShiftLatitude = $eShiftLatitude + ($fShiftLatitude - $eShiftLatitude) * $x; |
|
144
|
|
|
|
|
145
|
|
|
$eShiftLongitude = $recordLR['LONGITUDE_SHIFT'] + ($recordLL['LONGITUDE_SHIFT'] - $recordLR['LONGITUDE_SHIFT']) * $y; |
|
146
|
|
|
$fShiftLongitude = $recordUR['LONGITUDE_SHIFT'] + ($recordUL['LONGITUDE_SHIFT'] - $recordUR['LONGITUDE_SHIFT']) * $y; |
|
147
|
|
|
$pShiftLongitude = $eShiftLongitude + ($fShiftLongitude - $eShiftLongitude) * $x; |
|
148
|
|
|
|
|
149
|
|
|
return [new ArcSecond($pShiftLatitude), new ArcSecond(-$pShiftLongitude)]; // NTv2 is longitude positive *west* |
|
150
|
|
|
} |
|
151
|
|
|
|
|
152
|
|
|
private function getRecord(int $recordIndex): array |
|
153
|
|
|
{ |
|
154
|
|
|
$this->fseek($recordIndex); |
|
155
|
|
|
$rawRecord = $this->fread(self::RECORD_SIZE); |
|
156
|
|
|
$shifts = unpack("{$this->floatFormatChar}LATITUDE_SHIFT/{$this->floatFormatChar}LONGITUDE_SHIFT/{$this->floatFormatChar}LATITUDE_ACCURACY/{$this->floatFormatChar}LONGITUDE_ACCURACY", $rawRecord); |
|
157
|
|
|
|
|
158
|
|
|
return $shifts; |
|
159
|
|
|
} |
|
160
|
|
|
|
|
161
|
|
|
private function readHeader(): void |
|
162
|
|
|
{ |
|
163
|
|
|
$this->fseek(0); |
|
164
|
|
|
$rawData = $this->fread(11 * self::RECORD_SIZE); |
|
165
|
|
|
if (unpack('VNUM_OREC', $rawData, 8)['NUM_OREC'] !== 11) { |
|
166
|
|
|
$this->integerFormatChar = 'N'; |
|
167
|
|
|
$this->doubleFormatChar = 'E'; |
|
168
|
|
|
$this->floatFormatChar = 'G'; |
|
169
|
|
|
} |
|
170
|
|
|
|
|
171
|
|
|
$data = unpack("A8/{$this->integerFormatChar}NUM_OREC/x4/A8/{$this->integerFormatChar}NUM_SREC/x4/A8/{$this->integerFormatChar}NUM_FILE/x4/A8/A8GS_TYPE/A8/A8VERSION/A8/A8SYSTEM_F/A8/A8SYSTEM_T/A8/{$this->doubleFormatChar}MAJOR_F/A8/{$this->doubleFormatChar}MINOR_F/A8/{$this->doubleFormatChar}MAJOR_T/A8/{$this->doubleFormatChar}MINOR_T", $rawData); |
|
172
|
|
|
|
|
173
|
|
|
assert($data['GS_TYPE'] === 'SECONDS'); |
|
174
|
|
|
|
|
175
|
|
|
$subFileStart = 11 * self::RECORD_SIZE; |
|
176
|
|
|
for ($i = 0; $i < $data['NUM_FILE']; ++$i) { |
|
177
|
|
|
$this->fseek($subFileStart); |
|
178
|
|
|
$subFileRawData = $this->fread(11 * self::RECORD_SIZE); |
|
179
|
|
|
$subFileData = unpack("A8/A8SUB_NAME/A8/A8PARENT/A8/A8CREATED/A8/A8UPDATED/A8/{$this->doubleFormatChar}S_LAT/A8/{$this->doubleFormatChar}N_LAT/A8/{$this->doubleFormatChar}E_LONG/A8/{$this->doubleFormatChar}W_LONG/A8/{$this->doubleFormatChar}LAT_INC/A8/{$this->doubleFormatChar}LONG_INC/A8/{$this->integerFormatChar}GS_COUNT/x4", $subFileRawData); |
|
180
|
|
|
$subFileData['offsetStart'] = $subFileStart; |
|
181
|
|
|
|
|
182
|
|
|
//apply rounding to eliminate fp issues when being deserialized |
|
183
|
|
|
$subFileData['S_LAT'] = round($subFileData['S_LAT'], 5); |
|
184
|
|
|
$subFileData['N_LAT'] = round($subFileData['N_LAT'], 5); |
|
185
|
|
|
$subFileData['E_LONG'] = round($subFileData['E_LONG'], 5); |
|
186
|
|
|
$subFileData['W_LONG'] = round($subFileData['W_LONG'], 5); |
|
187
|
|
|
$this->subFileMetaData[$subFileData['SUB_NAME']] = $subFileData; |
|
188
|
|
|
|
|
189
|
|
|
$subFileStart += 11 * self::RECORD_SIZE + $subFileData['GS_COUNT'] * self::RECORD_SIZE; |
|
190
|
|
|
} |
|
191
|
|
|
} |
|
192
|
|
|
|
|
193
|
|
|
private function determineBestGrid(float $latitude, float $longitude): array |
|
194
|
|
|
{ |
|
195
|
|
|
$possibleGrids = []; |
|
196
|
|
|
foreach ($this->subFileMetaData as $subFileName => $subFileMetaDatum) { |
|
197
|
|
|
if ($latitude >= $subFileMetaDatum['S_LAT'] && $latitude <= $subFileMetaDatum['N_LAT'] && $longitude >= $subFileMetaDatum['E_LONG'] && $longitude <= $subFileMetaDatum['W_LONG']) { |
|
198
|
|
|
if ($latitude === $subFileMetaDatum['N_LAT'] && $longitude === $subFileMetaDatum['W_LONG']) { |
|
199
|
|
|
$possibleGrids[] = [self::FLAG_ON_UPPER_LATITUDE_AND_LONGITUDE, $subFileMetaDatum]; |
|
200
|
|
|
} elseif ($longitude === $subFileMetaDatum['W_LONG']) { |
|
201
|
|
|
$possibleGrids[] = [self::FLAG_ON_UPPER_LONGITUDE, $subFileMetaDatum]; |
|
202
|
|
|
} elseif ($latitude === $subFileMetaDatum['N_LAT']) { |
|
203
|
|
|
$possibleGrids[] = [self::FLAG_ON_UPPER_LATITUDE, $subFileMetaDatum]; |
|
204
|
|
|
} else { |
|
205
|
|
|
$possibleGrids[] = [self::FLAG_WITHIN_LIMITS, $subFileMetaDatum]; |
|
206
|
|
|
} |
|
207
|
|
|
} |
|
208
|
|
|
} |
|
209
|
|
|
|
|
210
|
|
|
if (!$possibleGrids) { |
|
211
|
|
|
throw new InvalidArgumentException('Specified coordinates are not within this grid file'); |
|
212
|
|
|
} |
|
213
|
|
|
|
|
214
|
|
|
usort($possibleGrids, static function ($a, $b) { |
|
215
|
|
|
return $a[0] <=> $b[0] ?: $a[1]['LAT_INC'] <=> $b[1]['LAT_INC'] ?: $a[2]['LONG_INC'] <=> $b[2]['LONG_INC']; |
|
216
|
|
|
}); |
|
217
|
|
|
|
|
218
|
|
|
return $possibleGrids[0]; |
|
219
|
|
|
} |
|
220
|
|
|
} |
|
221
|
|
|
|