|
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
|
|
|
private float $startLatitude; |
|
19
|
|
|
private float $startLongitude; |
|
20
|
|
|
private float $latitudeGridSize; |
|
21
|
|
|
private float $longitudeGridSize; |
|
22
|
|
|
private int $numberLatitudes; |
|
23
|
|
|
private int $numberLongitudes; |
|
24
|
|
|
private string $gridDataType; |
|
25
|
|
|
|
|
26
|
|
|
public function __construct($filename) |
|
27
|
|
|
{ |
|
28
|
|
|
parent::__construct($filename); |
|
29
|
|
|
|
|
30
|
|
|
$header = $this->getHeader(); |
|
31
|
|
|
$this->startLatitude = $header['xlatsw']; |
|
32
|
|
|
$this->startLongitude = $header['xlonsw']; |
|
33
|
|
|
$this->latitudeGridSize = $header['dlat']; |
|
34
|
|
|
$this->longitudeGridSize = $header['dlon']; |
|
35
|
|
|
$this->numberLatitudes = $header['nlat']; |
|
36
|
|
|
$this->numberLongitudes = $header['nlon']; |
|
37
|
|
|
|
|
38
|
|
|
switch ($header['ikind']) { |
|
39
|
|
|
case 1: |
|
40
|
|
|
$this->gridDataType = 'G'; |
|
41
|
|
|
break; |
|
42
|
|
|
default: |
|
43
|
|
|
throw new UnexpectedValueException("Unknown ikind: {$header['ikind']}"); |
|
44
|
|
|
} |
|
45
|
|
|
} |
|
46
|
|
|
|
|
47
|
|
|
public function getForwardAdjustment(GeographicPoint $point): float |
|
48
|
|
|
{ |
|
49
|
|
|
return $this->getAdjustment($point->getLatitude()->asDegrees()->getValue(), $point->getLongitude()->asDegrees()->getValue()); |
|
50
|
|
|
} |
|
51
|
|
|
|
|
52
|
|
|
/** |
|
53
|
|
|
* Converted from NOAA FORTRAN. |
|
54
|
|
|
*/ |
|
55
|
|
|
public function getAdjustment(float $latitude, float $longitude): float |
|
56
|
|
|
{ |
|
57
|
|
|
if ($longitude < 0) { |
|
58
|
|
|
$longitude += 360; // NADCON5 uses 0 = 360 = Greenwich |
|
59
|
|
|
} |
|
60
|
|
|
|
|
61
|
|
|
// Method: |
|
62
|
|
|
// Fit a 3x3 window over the random point. The closest |
|
63
|
|
|
// 2x2 points will surround the point. But based on which |
|
64
|
|
|
// quadrant of that 2x2 cell in which the point falls, the |
|
65
|
|
|
// 3x3 window could extend NW, NE, SW or SE from the 2x2 cell. |
|
66
|
|
|
|
|
67
|
|
|
// Find which row should be LEAST when fitting |
|
68
|
|
|
// a 3x3 window around $latitude. |
|
69
|
|
|
$difla = ($latitude - $this->startLatitude); |
|
70
|
|
|
$ratla = $difla / ($this->latitudeGridSize / 2); |
|
71
|
|
|
$ila = (int) ($ratla) + 1; |
|
72
|
|
|
if ($ila % 2 != 0) { |
|
73
|
|
|
$jla = ($ila + 1) / 2 - 1; |
|
74
|
|
|
} else { |
|
75
|
|
|
$jla = ($ila) / 2; |
|
76
|
|
|
} |
|
77
|
|
|
|
|
78
|
|
|
// Fix any edge overlaps |
|
79
|
|
|
if ($jla < 1) { |
|
80
|
|
|
$jla = 1; |
|
81
|
|
|
} |
|
82
|
|
|
if ($jla > $this->numberLatitudes - 2) { |
|
83
|
|
|
$jla = $this->numberLatitudes - 2; |
|
84
|
|
|
} |
|
85
|
|
|
|
|
86
|
|
|
// Find which column should be LEAST when fitting |
|
87
|
|
|
// a 3x3 window around $longitude. |
|
88
|
|
|
$diflo = ($longitude - $this->startLongitude); |
|
89
|
|
|
$ratlo = $diflo / ($this->longitudeGridSize / 2); |
|
90
|
|
|
$ilo = (int) ($ratlo) + 1; |
|
91
|
|
|
if ($ilo % 2 != 0) { |
|
92
|
|
|
$jlo = ($ilo + 1) / 2 - 1; |
|
93
|
|
|
} else { |
|
94
|
|
|
$jlo = ($ilo) / 2; |
|
95
|
|
|
} |
|
96
|
|
|
|
|
97
|
|
|
// Fix any edge overlaps |
|
98
|
|
|
if ($jlo < 1) { |
|
99
|
|
|
$jlo = 1; |
|
100
|
|
|
} |
|
101
|
|
|
if ($jlo > $this->numberLongitudes - 2) { |
|
102
|
|
|
$jlo = $this->numberLongitudes - 2; |
|
103
|
|
|
} |
|
104
|
|
|
|
|
105
|
|
|
// In the range of 0(westernmost) to |
|
106
|
|
|
// 2(easternmost) col, where is our |
|
107
|
|
|
// random lon value? That is, x must |
|
108
|
|
|
// be real and fall between 0 and 2. |
|
109
|
|
|
$x = ($longitude - $this->longitudeGridSize * ($jlo - 1) - $this->startLongitude) / $this->longitudeGridSize; |
|
110
|
|
|
|
|
111
|
|
|
// In the range of 0(southernmost) to |
|
112
|
|
|
// 2(northernmost) row, where is our |
|
113
|
|
|
// random lat value? That is, x must |
|
114
|
|
|
// be real and fall between 0 and 2. |
|
115
|
|
|
$y = ($latitude - $this->latitudeGridSize * ($jla - 1) - $this->startLatitude) / $this->latitudeGridSize; |
|
116
|
|
|
|
|
117
|
|
|
// Now do the interpolation. First, build a parabola |
|
118
|
|
|
// east-west the southermost row and interpolate to longitude |
|
119
|
|
|
// "xlo" (at "x" for 0<x<2). Then do it in the middle |
|
120
|
|
|
// row, then finally the northern row. The last step |
|
121
|
|
|
// is to fit a parabola north-south at the three previous |
|
122
|
|
|
// interpolations, but this time to interpolate to |
|
123
|
|
|
// latitude "xla" (at "y" for 0<y<2). Obviously we |
|
124
|
|
|
// could reverse the order, doing 3 north-south parabolas |
|
125
|
|
|
// followed by one east-east and we'd get the same answer. |
|
126
|
|
|
|
|
127
|
|
|
$fx0 = $this->qterp($x, $this->getRecord($jla, $jlo), $this->getRecord($jla, $jlo + 1), $this->getRecord($jla, $jlo + 2)); |
|
128
|
|
|
$fx1 = $this->qterp($x, $this->getRecord($jla + 1, $jlo), $this->getRecord($jla + 1, $jlo + 1), $this->getRecord($jla + 1, $jlo + 2)); |
|
129
|
|
|
$fx2 = $this->qterp($x, $this->getRecord($jla + 2, $jlo), $this->getRecord($jla + 2, $jlo + 1), $this->getRecord($jla + 2, $jlo + 2)); |
|
130
|
|
|
|
|
131
|
|
|
return $this->qterp($y, $fx0, $fx1, $fx2); |
|
132
|
|
|
} |
|
133
|
|
|
|
|
134
|
|
|
public function getRecord(int $latitudeIndex, int $longitudeIndex): float |
|
135
|
|
|
{ |
|
136
|
|
|
$startBytes = 52; |
|
137
|
|
|
$dataTypeBytes = $this->gridDataType === 'n' ? 2 : 4; |
|
138
|
|
|
$recordLength = 4 + $this->numberLongitudes * $dataTypeBytes + 4; |
|
139
|
|
|
|
|
140
|
|
|
$this->fseek($startBytes + $recordLength * ($latitudeIndex - 1)); |
|
141
|
|
|
$rawRow = $this->fread($recordLength); |
|
142
|
|
|
$row = unpack("Gstartbuffer/{$this->gridDataType}{$this->numberLongitudes}lon/Gendbuffer", $rawRow); |
|
143
|
|
|
|
|
144
|
|
|
return $row['lon' . ($longitudeIndex)]; |
|
145
|
|
|
} |
|
146
|
|
|
|
|
147
|
|
|
private function getHeader(): array |
|
148
|
|
|
{ |
|
149
|
|
|
$this->fseek(0); |
|
150
|
|
|
$rawData = $this->fread(52); |
|
151
|
|
|
$data = unpack('Gstartbuffer/Exlatsw/Exlonsw/Edlat/Edlon/Nnlat/Nnlon/Nikind/Gendbuffer/', $rawData); |
|
152
|
|
|
|
|
153
|
|
|
return $data; |
|
154
|
|
|
} |
|
155
|
|
|
|
|
156
|
|
|
/** |
|
157
|
|
|
* Converted from NOAA FORTRAN. |
|
158
|
|
|
* This function fits a parabola (quadratic) through |
|
159
|
|
|
* three points, *equally* spaced along the x-axis |
|
160
|
|
|
* at indices 0, 1 and 2. The spacing along the |
|
161
|
|
|
* x-axis is "dx" |
|
162
|
|
|
* Thus:. |
|
163
|
|
|
* |
|
164
|
|
|
* f0 = y(x(0)) |
|
165
|
|
|
* f1 = y(x(1)) |
|
166
|
|
|
* f2 = y(x(2)) |
|
167
|
|
|
* Where |
|
168
|
|
|
* x(1) = x(0) + dx |
|
169
|
|
|
* x(2) = x(1) + dx |
|
170
|
|
|
* The input value is some value of "x" that falls |
|
171
|
|
|
* between 0 and 2. The output value is |
|
172
|
|
|
* the parabolic function at x. |
|
173
|
|
|
* |
|
174
|
|
|
* This function uses Newton-Gregory forward polynomial. |
|
175
|
|
|
*/ |
|
176
|
|
|
private function qterp(float $x, float $f0, float $f1, float $f2): float |
|
177
|
|
|
{ |
|
178
|
|
|
$df0 = $f1 - $f0; |
|
179
|
|
|
$df1 = $f2 - $f1; |
|
180
|
|
|
$d2f0 = $df1 - $df0; |
|
181
|
|
|
|
|
182
|
|
|
return $f0 + $x * $df0 + 0.5 * $x * ($x - 1.0) * $d2f0; |
|
183
|
|
|
} |
|
184
|
|
|
} |
|
185
|
|
|
|
dvdoug /
PHPCoord
