|
1
|
|
|
<?php |
|
2
|
|
|
/** |
|
3
|
|
|
* PHPCoord. |
|
4
|
|
|
* |
|
5
|
|
|
* @author Doug Wright |
|
6
|
|
|
*/ |
|
7
|
|
|
declare(strict_types=1); |
|
8
|
|
|
|
|
9
|
|
|
namespace PHPCoord\CoordinateOperation; |
|
10
|
|
|
|
|
11
|
|
|
use PHPCoord\GeographicPoint; |
|
12
|
|
|
use SplFileObject; |
|
13
|
|
|
use UnexpectedValueException; |
|
14
|
|
|
use function unpack; |
|
15
|
|
|
|
|
16
|
|
|
class NADCON5Grid extends SplFileObject |
|
17
|
|
|
{ |
|
18
|
|
|
use BiquadraticInterpolation; |
|
19
|
|
|
|
|
20
|
|
|
private string $gridDataType; |
|
21
|
|
|
|
|
22
|
7 |
|
public function __construct($filename) |
|
23
|
|
|
{ |
|
24
|
7 |
|
parent::__construct($filename); |
|
25
|
|
|
|
|
26
|
7 |
|
$header = $this->getHeader(); |
|
27
|
7 |
|
$this->startX = $header['xlonsw']; |
|
28
|
7 |
|
$this->startY = $header['xlatsw']; |
|
29
|
7 |
|
$this->columnGridInterval = $header['dlon']; |
|
30
|
7 |
|
$this->rowGridInterval = $header['dlat']; |
|
31
|
7 |
|
$this->numberOfColumns = $header['nlon']; |
|
32
|
7 |
|
$this->numberOfRows = $header['nlat']; |
|
33
|
|
|
|
|
34
|
7 |
|
switch ($header['ikind']) { |
|
35
|
7 |
|
case 1: |
|
36
|
7 |
|
$this->gridDataType = 'G'; |
|
37
|
7 |
|
break; |
|
38
|
|
|
default: |
|
39
|
|
|
throw new UnexpectedValueException("Unknown ikind: {$header['ikind']}"); |
|
40
|
|
|
} |
|
41
|
7 |
|
} |
|
42
|
|
|
|
|
43
|
7 |
|
public function getForwardAdjustment(GeographicPoint $point): float |
|
44
|
|
|
{ |
|
45
|
7 |
|
return $this->getAdjustment($point->getLatitude()->asDegrees()->getValue(), $point->getLongitude()->asDegrees()->getValue()); |
|
46
|
|
|
} |
|
47
|
|
|
|
|
48
|
|
|
/** |
|
49
|
|
|
* Converted from NOAA FORTRAN. |
|
50
|
|
|
*/ |
|
51
|
7 |
|
private function getAdjustment(float $latitude, float $longitude): float |
|
52
|
|
|
{ |
|
53
|
7 |
|
if ($longitude < 0) { |
|
54
|
7 |
|
$longitude += 360; // NADCON5 uses 0 = 360 = Greenwich |
|
55
|
|
|
} |
|
56
|
|
|
|
|
57
|
|
|
// Method: |
|
58
|
|
|
// Fit a 3x3 window over the random point. The closest |
|
59
|
|
|
// 2x2 points will surround the point. But based on which |
|
60
|
|
|
// quadrant of that 2x2 cell in which the point falls, the |
|
61
|
|
|
// 3x3 window could extend NW, NE, SW or SE from the 2x2 cell. |
|
62
|
|
|
|
|
63
|
|
|
// Find which column should be LEAST when fitting |
|
64
|
|
|
// a 3x3 window around $longitude. |
|
65
|
7 |
|
$diflo = ($longitude - $this->startX); |
|
66
|
7 |
|
$ratlo = $diflo / ($this->columnGridInterval / 2); |
|
67
|
7 |
|
$ilo = (int) ($ratlo) + 1; |
|
68
|
7 |
|
if ($ilo % 2 != 0) { |
|
69
|
|
|
$xIndex = ($ilo + 1) / 2 - 1; |
|
70
|
|
|
} else { |
|
71
|
7 |
|
$xIndex = ($ilo) / 2; |
|
72
|
|
|
} |
|
73
|
|
|
|
|
74
|
|
|
// Fix any edge overlaps |
|
75
|
7 |
|
if ($xIndex < 1) { |
|
76
|
|
|
$xIndex = 1; |
|
77
|
|
|
} |
|
78
|
7 |
|
if ($xIndex > $this->numberOfColumns - 2) { |
|
79
|
|
|
$xIndex = $this->numberOfColumns - 2; |
|
80
|
|
|
} |
|
81
|
|
|
|
|
82
|
|
|
// Find which row should be LEAST when fitting |
|
83
|
|
|
// a 3x3 window around $latitude. |
|
84
|
7 |
|
$difla = ($latitude - $this->startY); |
|
85
|
7 |
|
$ratla = $difla / ($this->rowGridInterval / 2); |
|
86
|
7 |
|
$ila = (int) ($ratla) + 1; |
|
87
|
7 |
|
if ($ila % 2 != 0) { |
|
88
|
4 |
|
$yIndex = ($ila + 1) / 2 - 1; |
|
89
|
|
|
} else { |
|
90
|
3 |
|
$yIndex = ($ila) / 2; |
|
91
|
|
|
} |
|
92
|
|
|
|
|
93
|
|
|
// Fix any edge overlaps |
|
94
|
7 |
|
if ($yIndex < 1) { |
|
95
|
|
|
$yIndex = 1; |
|
96
|
|
|
} |
|
97
|
7 |
|
if ($yIndex > $this->numberOfRows - 2) { |
|
98
|
|
|
$yIndex = $this->numberOfRows - 2; |
|
99
|
|
|
} |
|
100
|
|
|
|
|
101
|
|
|
// In the range of 0(westernmost) to |
|
102
|
|
|
// 2(easternmost) col, where is our |
|
103
|
|
|
// random lon value? That is, x must |
|
104
|
|
|
// be real and fall between 0 and 2. |
|
105
|
7 |
|
$x = ($longitude - $this->columnGridInterval * ($xIndex - 1) - $this->startX) / $this->columnGridInterval; |
|
106
|
|
|
|
|
107
|
|
|
// In the range of 0(southernmost) to |
|
108
|
|
|
// 2(northernmost) row, where is our |
|
109
|
|
|
// random lat value? That is, x must |
|
110
|
|
|
// be real and fall between 0 and 2. |
|
111
|
7 |
|
$y = ($latitude - $this->rowGridInterval * ($yIndex - 1) - $this->startY) / $this->rowGridInterval; |
|
112
|
|
|
|
|
113
|
7 |
|
$corners = [ |
|
114
|
7 |
|
'lowerLeft' => $this->getRecord($xIndex, $yIndex), |
|
115
|
7 |
|
'lowerCentre' => $this->getRecord($xIndex + 1, $yIndex), |
|
116
|
7 |
|
'lowerRight' => $this->getRecord($xIndex + 2, $yIndex), |
|
117
|
7 |
|
'middleLeft' => $this->getRecord($xIndex, $yIndex + 1), |
|
118
|
7 |
|
'middleCentre' => $this->getRecord($xIndex + 1, $yIndex + 1), |
|
119
|
7 |
|
'middleRight' => $this->getRecord($xIndex + 2, $yIndex + 1), |
|
120
|
7 |
|
'upperLeft' => $this->getRecord($xIndex, $yIndex + 2), |
|
121
|
7 |
|
'upperCentre' => $this->getRecord($xIndex + 1, $yIndex + 2), |
|
122
|
7 |
|
'upperRight' => $this->getRecord($xIndex + 2, $yIndex + 2), |
|
123
|
|
|
]; |
|
124
|
|
|
|
|
125
|
|
|
//Interpolate value at lower row |
|
126
|
7 |
|
$y0 = $this->interpolateQuadratic($x, $corners['lowerLeft'], $corners['lowerCentre'], $corners['lowerRight']); |
|
127
|
|
|
//Interpolate value at middle row |
|
128
|
7 |
|
$y1 = $this->interpolateQuadratic($x, $corners['middleLeft'], $corners['middleCentre'], $corners['middleRight']); |
|
129
|
|
|
//Interpolate value at upper row |
|
130
|
7 |
|
$y2 = $this->interpolateQuadratic($x, $corners['upperLeft'], $corners['upperCentre'], $corners['upperRight']); |
|
131
|
|
|
|
|
132
|
|
|
//Interpolate between rows |
|
133
|
7 |
|
return $this->interpolateQuadratic($y, $y0, $y1, $y2); |
|
134
|
|
|
} |
|
135
|
|
|
|
|
136
|
7 |
|
public function getRecord(int $longitudeIndex, int $latitudeIndex): float |
|
137
|
|
|
{ |
|
138
|
7 |
|
$startBytes = 52; |
|
139
|
7 |
|
$dataTypeBytes = $this->gridDataType === 'n' ? 2 : 4; |
|
140
|
7 |
|
$recordLength = 4 + $this->numberOfColumns * $dataTypeBytes + 4; |
|
141
|
|
|
|
|
142
|
7 |
|
$this->fseek($startBytes + $recordLength * ($latitudeIndex - 1)); |
|
143
|
7 |
|
$rawRow = $this->fread($recordLength); |
|
144
|
7 |
|
$row = unpack("Gstartbuffer/{$this->gridDataType}{$this->numberOfColumns}lon/Gendbuffer", $rawRow); |
|
145
|
|
|
|
|
146
|
7 |
|
return $row['lon' . ($longitudeIndex)]; |
|
147
|
|
|
} |
|
148
|
|
|
|
|
149
|
7 |
|
private function getHeader(): array |
|
150
|
|
|
{ |
|
151
|
7 |
|
$this->fseek(0); |
|
152
|
7 |
|
$rawData = $this->fread(52); |
|
153
|
7 |
|
$data = unpack('Gstartbuffer/Exlatsw/Exlonsw/Edlat/Edlon/Nnlat/Nnlon/Nikind/Gendbuffer/', $rawData); |
|
154
|
|
|
|
|
155
|
7 |
|
return $data; |
|
156
|
|
|
} |
|
157
|
|
|
} |
|
158
|
|
|
|
dvdoug /
PHPCoord
