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 PHPCoord\CoordinateReferenceSystem\Geographic; |
14
|
|
|
use PHPCoord\GeographicPoint; |
15
|
|
|
use PHPCoord\UnitOfMeasure\Angle\Angle; |
16
|
|
|
use PHPCoord\UnitOfMeasure\Angle\ArcSecond; |
17
|
|
|
use function round; |
18
|
|
|
use SplFileObject; |
19
|
|
|
use function unpack; |
20
|
|
|
use function usort; |
21
|
|
|
|
22
|
|
|
class NTv2Grid extends SplFileObject |
23
|
|
|
{ |
24
|
|
|
private const RECORD_SIZE = 16; |
25
|
|
|
private const ITERATION_CONVERGENCE = 0.0001; |
26
|
|
|
private const FLAG_WITHIN_LIMITS = 1; |
27
|
|
|
private const FLAG_ON_UPPER_LATITUDE = 2; |
28
|
|
|
private const FLAG_ON_UPPER_LONGITUDE = 3; |
29
|
|
|
private const FLAG_ON_UPPER_LATITUDE_AND_LONGITUDE = 4; |
30
|
|
|
|
31
|
|
|
private string $integerFormatChar = 'V'; |
32
|
|
|
private string $doubleFormatChar = 'e'; |
33
|
|
|
private string $floatFormatChar = 'g'; |
34
|
|
|
|
35
|
|
|
private array $subFileMetaData = []; |
36
|
|
|
|
37
|
5 |
|
public function __construct($filename) |
38
|
|
|
{ |
39
|
5 |
|
parent::__construct($filename); |
40
|
|
|
|
41
|
5 |
|
$this->readHeader(); |
42
|
5 |
|
} |
43
|
|
|
|
44
|
3 |
|
public function applyForwardAdjustment(GeographicPoint $point, Geographic $to): GeographicPoint |
45
|
|
|
{ |
46
|
3 |
|
$adjustment = $this->getAdjustment($point->getLatitude(), $point->getLongitude()); |
47
|
|
|
|
48
|
3 |
|
$latitude = $point->getLatitude()->add($adjustment[0]); |
49
|
3 |
|
$longitude = $point->getLongitude()->add($adjustment[1]); |
50
|
|
|
|
51
|
3 |
|
return GeographicPoint::create($latitude, $longitude, $point->getHeight(), $to, $point->getCoordinateEpoch()); |
52
|
|
|
} |
53
|
|
|
|
54
|
2 |
|
public function applyReverseAdjustment(GeographicPoint $point, Geographic $to): GeographicPoint |
55
|
|
|
{ |
56
|
2 |
|
$adjustment = [new ArcSecond(0), new ArcSecond(0)]; |
57
|
2 |
|
$latitude = $point->getLatitude(); |
58
|
2 |
|
$longitude = $point->getLongitude(); |
59
|
|
|
|
60
|
|
|
do { |
61
|
2 |
|
$prevAdjustment = $adjustment; |
62
|
2 |
|
$adjustment = $this->getAdjustment($latitude, $longitude); |
63
|
2 |
|
$latitude = $point->getLatitude()->subtract($adjustment[0]); |
64
|
2 |
|
$longitude = $point->getLongitude()->subtract($adjustment[1]); |
65
|
2 |
|
} while (abs($adjustment[0]->subtract($prevAdjustment[0])->getValue()) > self::ITERATION_CONVERGENCE && abs($adjustment[1]->subtract($prevAdjustment[1])->getValue()) > self::ITERATION_CONVERGENCE); |
66
|
|
|
|
67
|
2 |
|
return GeographicPoint::create($latitude, $longitude, $point->getHeight(), $to, $point->getCoordinateEpoch()); |
68
|
|
|
} |
69
|
|
|
|
70
|
|
|
/** |
71
|
|
|
* @return ArcSecond[] |
72
|
|
|
*/ |
73
|
5 |
|
private function getAdjustment(Angle $latitude, Angle $longitude): array |
74
|
|
|
{ |
75
|
|
|
// NTv2 is longitude positive *west* |
76
|
5 |
|
$longitude = $longitude->multiply(-1); |
77
|
|
|
|
78
|
5 |
|
$latitudeAsSeconds = Angle::convert($latitude, Angle::EPSG_ARC_SECOND)->getValue(); |
79
|
5 |
|
$longitudeAsSeconds = Angle::convert($longitude, Angle::EPSG_ARC_SECOND)->getValue(); |
80
|
5 |
|
$gridToUse = $this->determineBestGrid($latitudeAsSeconds, $longitudeAsSeconds); |
81
|
|
|
|
82
|
5 |
|
$offsets = $gridToUse->interpolateBilinear($longitudeAsSeconds, $latitudeAsSeconds); |
83
|
|
|
|
84
|
5 |
|
return [new ArcSecond($offsets[0]), new ArcSecond(-$offsets[1])]; // NTv2 is longitude positive *west* |
85
|
|
|
} |
86
|
|
|
|
87
|
5 |
|
private function readHeader(): void |
88
|
|
|
{ |
89
|
5 |
|
$this->fseek(0); |
90
|
5 |
|
$rawData = $this->fread(11 * self::RECORD_SIZE); |
91
|
5 |
|
if (unpack('VNUM_OREC', $rawData, 8)['NUM_OREC'] !== 11) { |
92
|
|
|
$this->integerFormatChar = 'N'; |
93
|
|
|
$this->doubleFormatChar = 'E'; |
94
|
|
|
$this->floatFormatChar = 'G'; |
95
|
|
|
} |
96
|
|
|
|
97
|
5 |
|
$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); |
98
|
|
|
|
99
|
5 |
|
assert($data['GS_TYPE'] === 'SECONDS'); |
100
|
|
|
|
101
|
5 |
|
$subFileStart = 11 * self::RECORD_SIZE; |
102
|
5 |
|
for ($i = 0; $i < $data['NUM_FILE']; ++$i) { |
103
|
5 |
|
$this->fseek($subFileStart); |
104
|
5 |
|
$subFileRawData = $this->fread(11 * self::RECORD_SIZE); |
105
|
5 |
|
$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); |
106
|
5 |
|
$subFileData['offsetStart'] = $subFileStart; |
107
|
|
|
|
108
|
|
|
//apply rounding to eliminate fp issues when being deserialized |
109
|
5 |
|
$subFileData['S_LAT'] = round($subFileData['S_LAT'], 5); |
110
|
5 |
|
$subFileData['N_LAT'] = round($subFileData['N_LAT'], 5); |
111
|
5 |
|
$subFileData['E_LONG'] = round($subFileData['E_LONG'], 5); |
112
|
5 |
|
$subFileData['W_LONG'] = round($subFileData['W_LONG'], 5); |
113
|
5 |
|
$this->subFileMetaData[$subFileData['SUB_NAME']] = $subFileData; |
114
|
|
|
|
115
|
5 |
|
$subFileStart += 11 * self::RECORD_SIZE + $subFileData['GS_COUNT'] * self::RECORD_SIZE; |
116
|
|
|
} |
117
|
5 |
|
} |
118
|
|
|
|
119
|
5 |
|
private function determineBestGrid(float $latitude, float $longitude): NTv2SubGrid |
120
|
|
|
{ |
121
|
5 |
|
$possibleGrids = []; |
122
|
5 |
|
foreach ($this->subFileMetaData as $subFileMetaDatum) { |
123
|
5 |
|
if ($latitude === $subFileMetaDatum['N_LAT'] && $longitude === $subFileMetaDatum['W_LONG']) { |
124
|
|
|
$possibleGrids[] = [self::FLAG_ON_UPPER_LATITUDE_AND_LONGITUDE, $subFileMetaDatum]; |
125
|
5 |
|
} elseif ($longitude === $subFileMetaDatum['W_LONG']) { |
126
|
|
|
$possibleGrids[] = [self::FLAG_ON_UPPER_LONGITUDE, $subFileMetaDatum]; |
127
|
5 |
|
} elseif ($latitude === $subFileMetaDatum['N_LAT']) { |
128
|
|
|
$possibleGrids[] = [self::FLAG_ON_UPPER_LATITUDE, $subFileMetaDatum]; |
129
|
5 |
|
} elseif ($latitude >= $subFileMetaDatum['S_LAT'] && $latitude <= $subFileMetaDatum['N_LAT'] && $longitude >= $subFileMetaDatum['E_LONG'] && $longitude <= $subFileMetaDatum['W_LONG']) { |
130
|
5 |
|
$possibleGrids[] = [self::FLAG_WITHIN_LIMITS, $subFileMetaDatum]; |
131
|
|
|
} |
132
|
|
|
} |
133
|
|
|
|
134
|
5 |
|
usort($possibleGrids, static function ($a, $b) { |
135
|
3 |
|
return $a[0] <=> $b[0] ?: $a[1]['LAT_INC'] <=> $b[1]['LAT_INC'] ?: $a[2]['LONG_INC'] <=> $b[2]['LONG_INC']; |
136
|
5 |
|
}); |
137
|
|
|
|
138
|
5 |
|
$gridToUse = $possibleGrids[0][1]; |
139
|
|
|
|
140
|
5 |
|
return new NTv2SubGrid( |
141
|
5 |
|
$this->getPathname(), |
142
|
5 |
|
$gridToUse['offsetStart'], |
143
|
5 |
|
$gridToUse['S_LAT'], |
144
|
5 |
|
$gridToUse['N_LAT'], |
145
|
5 |
|
$gridToUse['E_LONG'], |
146
|
5 |
|
$gridToUse['W_LONG'], |
147
|
5 |
|
$gridToUse['LAT_INC'], |
148
|
5 |
|
$gridToUse['LONG_INC'], |
149
|
5 |
|
$this->floatFormatChar |
150
|
|
|
); |
151
|
|
|
} |
152
|
|
|
} |
153
|
|
|
|