|
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
|
|
|
|
dvdoug /
PHPCoord
