|
1
|
|
|
<?php |
|
2
|
|
|
/* |
|
3
|
|
|
* Calculate the height of the WGS84 geoid above the |
|
4
|
|
|
* ellipsoid at any given latitude and longitude |
|
5
|
|
|
* |
|
6
|
|
|
* Copyright (c) Charles Karney (2009-2015) <[email protected]> and licensed |
|
7
|
|
|
* under the MIT/X11 License. For more information, see |
|
8
|
|
|
* https://geographiclib.sourceforge.io/ |
|
9
|
|
|
*/ |
|
10
|
|
|
/* |
|
11
|
|
|
* Translated to PHP of GeographicLib/src/Geoid.cpp |
|
12
|
|
|
* by Ycarus <[email protected]> in 2017 |
|
13
|
|
|
*/ |
|
14
|
|
|
class GeoidHeight { |
|
15
|
|
|
public $c0 = 240; |
|
16
|
|
|
public $c3 = [[9, -18, -88, 0, 96, 90, 0, 0, -60, -20], [-9, 18, 8, 0, -96, 30, 0, 0, 60, -20], [9, -88, -18, 90, 96, 0, -20, -60, 0, 0], [186, -42, -42, -150, -96, -150, 60, 60, 60, 60], [54, 162, -78, 30, -24, -90, -60, 60, -60, 60], [-9, -32, 18, 30, 24, 0, 20, -60, 0, 0], [-9, 8, 18, 30, -96, 0, -20, 60, 0, 0], [54, -78, 162, -90, -24, 30, 60, -60, 60, -60], [-54, 78, 78, 90, 144, 90, -60, -60, -60, -60], [9, -8, -18, -30, -24, 0, 20, 60, 0, 0], [-9, 18, -32, 0, 24, 30, 0, 0, -60, 20], [9, -18, -8, 0, -24, -30, 0, 0, 60, 20]]; |
|
17
|
|
|
public $c0n = 372; |
|
18
|
|
|
public $c3n = [[0, 0, -131, 0, 138, 144, 0, 0, -102, -31], [0, 0, 7, 0, -138, 42, 0, 0, 102, -31], [62, 0, -31, 0, 0, -62, 0, 0, 0, 31], [124, 0, -62, 0, 0, -124, 0, 0, 0, 62], [124, 0, -62, 0, 0, -124, 0, 0, 0, 62], [62, 0, -31, 0, 0, -62, 0, 0, 0, 31], [0, 0, 45, 0, -183, -9, 0, 93, 18, 0], [0, 0, 216, 0, 33, 87, 0, -93, 12, -93], [0, 0, 156, 0, 153, 99, 0, -93, -12, -93], [0, 0, -45, 0, -3, 9, 0, 93, -18, 0], [0, 0, -55, 0, 48, 42, 0, 0, -84, 31], [0, 0, -7, 0, -48, -42, 0, 0, 84, 31]]; |
|
19
|
|
|
public $c0s = 372; |
|
20
|
|
|
public $c3s = [[18, -36, -122, 0, 120, 135, 0, 0, -84, -31], [-18, 36, -2, 0, -120, 51, 0, 0, 84, -31], [36, -165, -27, 93, 147, -9, 0, -93, 18, 0], [210, 45, -111, -93, -57, -192, 0, 93, 12, 93], [162, 141, -75, -93, -129, -180, 0, 93, -12, 93], [-36, -21, 27, 93, 39, 9, 0, -93, -18, 0], [0, 0, 62, 0, 0, 31, 0, 0, 0, -31], [0, 0, 124, 0, 0, 62, 0, 0, 0, -62], [0, 0, 124, 0, 0, 62, 0, 0, 0, -62], [0, 0, 62, 0, 0, 31, 0, 0, 0, -31], [-18, 36, -64, 0, 66, 51, 0, 0, -102, 31], [18, -36, 2, 0, -66, -51, 0, 0, 102, 31]]; |
|
21
|
|
|
|
|
22
|
|
|
function __construct($name='egm2008-1.pgm') { |
|
|
|
|
|
|
23
|
|
|
if ($name == '') $name = dirname(__FILE__).'/install/tmp/egm2008-1.pgm'; |
|
24
|
|
|
$this->offset = null; |
|
|
|
|
|
|
25
|
|
|
$this->scale = null; |
|
|
|
|
|
|
26
|
|
|
|
|
27
|
|
|
$f = @fopen($name,"r"); |
|
28
|
|
|
if ($f === FALSE) { |
|
29
|
|
|
throw new Exception("Can't open ".$name); |
|
30
|
|
|
} |
|
31
|
|
|
$line = fgets($f,4096); |
|
32
|
|
|
if (trim($line) != 'P5') { |
|
33
|
|
|
throw new Exception('No PGM header'); |
|
34
|
|
|
} |
|
35
|
|
|
$headerlen = strlen($line); |
|
36
|
|
|
while (true) { |
|
37
|
|
|
$line = fgets($f,4096); |
|
38
|
|
|
if ((strlen($line) == 0)) { |
|
39
|
|
|
throw new Exception('EOF before end of file header'); |
|
40
|
|
|
} |
|
41
|
|
|
$headerlen += strlen($line); |
|
42
|
|
|
if (strpos($line,'# Offset ') !== FALSE) { |
|
43
|
|
|
try { |
|
44
|
|
|
$this->offset = substr($line, 9); |
|
45
|
|
|
} catch(ValueError $e) { |
|
|
|
|
|
|
46
|
|
|
throw new Exception('Error reading offset '.$e); |
|
47
|
|
|
} |
|
48
|
|
|
} else if (strpos($line,'# Scale ') !== FALSE) { |
|
49
|
|
|
try { |
|
50
|
|
|
$this->scale = substr($line, 8); |
|
51
|
|
|
} catch(ValueError $e) { |
|
|
|
|
|
|
52
|
|
|
throw new Exception('Error reading scale '.$e); |
|
53
|
|
|
} |
|
54
|
|
|
} else if ((strpos($line,'#') === FALSE)) { |
|
55
|
|
|
try { |
|
56
|
|
|
list($this->width, $this->height) = preg_split('/\s+/',$line); |
|
|
|
|
|
|
57
|
|
|
} catch(ValueError $e) { |
|
|
|
|
|
|
58
|
|
|
throw new Exception('Bad PGM width&height line'. $e); |
|
59
|
|
|
} |
|
60
|
|
|
break; |
|
61
|
|
|
} |
|
62
|
|
|
} |
|
63
|
|
|
$line = fgets($f,4096); |
|
64
|
|
|
$headerlen += strlen($line); |
|
65
|
|
|
$levels = (int)$line; |
|
66
|
|
|
$this->width = (int)$this->width; |
|
67
|
|
|
$this->height = (int)$this->height; |
|
68
|
|
|
if (($levels != 65535)) { |
|
69
|
|
|
throw new Exception('PGM file must have 65535 gray levels ('.$levels.')'); |
|
70
|
|
|
} |
|
71
|
|
|
if (($this->offset == null)) { |
|
|
|
|
|
|
72
|
|
|
throw new Exception('PGM file does not contain offset'); |
|
73
|
|
|
} |
|
74
|
|
|
if (($this->scale == null)) { |
|
|
|
|
|
|
75
|
|
|
throw new Exception('PGM file does not contain scale'); |
|
76
|
|
|
} |
|
77
|
|
|
if (($this->width < 2) || ($this->height < 2)) { |
|
78
|
|
|
throw new Exception('Raster size too small'); |
|
79
|
|
|
} |
|
80
|
|
|
|
|
81
|
|
|
$fullsize = filesize($name); |
|
82
|
|
|
if ((($fullsize - $headerlen) != (($this->width * $this->height) * 2))) { |
|
83
|
|
|
throw new Exception('File has the wrong length'); |
|
84
|
|
|
} |
|
85
|
|
|
|
|
86
|
|
|
$this->headerlen = $headerlen; |
|
|
|
|
|
|
87
|
|
|
$this->raw= $f; |
|
|
|
|
|
|
88
|
|
|
$this->rlonres = ($this->width / 360.0); |
|
|
|
|
|
|
89
|
|
|
$this->rlatres = (($this->height - 1) / 180.0); |
|
|
|
|
|
|
90
|
|
|
$this->ix = null; |
|
|
|
|
|
|
91
|
|
|
$this->iy = null; |
|
|
|
|
|
|
92
|
|
|
} |
|
93
|
|
|
|
|
94
|
|
|
function _rawval($ix,$iy) { |
|
|
|
|
|
|
95
|
|
|
if (($iy < 0)) { |
|
96
|
|
|
$iy = -$iy; |
|
97
|
|
|
$ix += ($this->width / 2); |
|
98
|
|
|
} else if (($iy >= $this->height)) { |
|
99
|
|
|
$iy = ((2 * ($this->height - 1)) - $iy); |
|
100
|
|
|
$ix += ($this->width / 2); |
|
101
|
|
|
} |
|
102
|
|
|
if (($ix < 0)) { |
|
103
|
|
|
$ix += $this->width; |
|
104
|
|
|
} else if (($ix >= $this->width)) { |
|
105
|
|
|
$ix -= $this->width; |
|
106
|
|
|
} |
|
107
|
|
|
$k = ((($iy * $this->width) + $ix) * 2) + $this->headerlen; |
|
108
|
|
|
fseek($this->raw,$k); |
|
109
|
|
|
return unpack('n',fread($this->raw,2))[1]; |
|
110
|
|
|
} |
|
111
|
|
|
|
|
112
|
|
|
function get($lat,$lon,$cubic=true) { |
|
|
|
|
|
|
113
|
|
|
if (($lon < 0)) { |
|
114
|
|
|
$lon += 360; |
|
115
|
|
|
} |
|
116
|
|
|
$fy = ((90 - $lat) * $this->rlatres); |
|
117
|
|
|
$fx = ($lon * $this->rlonres); |
|
118
|
|
|
$iy = (int)$fy; |
|
119
|
|
|
$ix = (int)$fx; |
|
120
|
|
|
$fx -= $ix; |
|
121
|
|
|
$fy -= $iy; |
|
122
|
|
|
if (($iy == ($this->height - 1))) { |
|
123
|
|
|
$iy -= 1; |
|
124
|
|
|
} |
|
125
|
|
|
if (($ix != $this->ix) || ($iy != $this->iy)) { |
|
126
|
|
|
$this->ix = $ix; |
|
127
|
|
|
$this->iy = $iy; |
|
128
|
|
|
if (!($cubic)) { |
|
129
|
|
|
$this->v00 = $this->_rawval($ix, $iy); |
|
|
|
|
|
|
130
|
|
|
$this->v01 = $this->_rawval(($ix + 1), $iy); |
|
|
|
|
|
|
131
|
|
|
$this->v10 = $this->_rawval($ix, ($iy + 1)); |
|
|
|
|
|
|
132
|
|
|
$this->v11 = $this->_rawval(($ix + 1), ($iy + 1)); |
|
|
|
|
|
|
133
|
|
|
} else { |
|
134
|
|
|
$v = [$this->_rawval($ix, ($iy - 1)), $this->_rawval(($ix + 1), ($iy - 1)), $this->_rawval(($ix - 1), $iy), $this->_rawval($ix, $iy), $this->_rawval(($ix + 1), $iy), $this->_rawval(($ix + 2), $iy), $this->_rawval(($ix - 1), ($iy + 1)), $this->_rawval($ix, ($iy + 1)), $this->_rawval(($ix + 1), ($iy + 1)), $this->_rawval(($ix + 2), ($iy + 1)), $this->_rawval($ix, ($iy + 2)), $this->_rawval(($ix + 1), ($iy + 2))]; |
|
135
|
|
|
if (($iy == 0)) { |
|
136
|
|
|
$c3x = $this->c3n; |
|
137
|
|
|
$c0x = $this->c0n; |
|
138
|
|
|
} else if (($iy == ($this->height - 2))) { |
|
139
|
|
|
$c3x = $this->c3s; |
|
140
|
|
|
$c0x = $this->c0s; |
|
141
|
|
|
} else { |
|
142
|
|
|
$c3x = $this->c3; |
|
143
|
|
|
$c0x = $this->c0; |
|
144
|
|
|
} |
|
145
|
|
|
for ($i = 0; $i < 10;++$i) { |
|
146
|
|
|
$this->t[$i] = 0; |
|
|
|
|
|
|
147
|
|
|
for ($j = 0; $j < 12; ++$j) { |
|
148
|
|
|
$this->t[$i] += $v[$j]*$c3x[$j][$i]; |
|
149
|
|
|
} |
|
150
|
|
|
$this->t[$i] /= $c0x; |
|
151
|
|
|
} |
|
152
|
|
|
} |
|
153
|
|
|
} |
|
154
|
|
|
if (!($cubic)) { |
|
155
|
|
|
$a = (((1 - $fx) * $this->v00) + ($fx * $this->v01)); |
|
156
|
|
|
$b = (((1 - $fx) * $this->v10) + ($fx * $this->v11)); |
|
157
|
|
|
$h = (((1 - $fy) * $a) + ($fy * $b)); |
|
158
|
|
|
} else { |
|
159
|
|
|
$h = (($this->t[0] + ($fx * ($this->t[1] + ($fx * ($this->t[3] + ($fx * $this->t[6])))))) + ($fy * (($this->t[2] + ($fx * ($this->t[4] + ($fx * $this->t[7])))) + ($fy * (($this->t[5] + ($fx * $this->t[8])) + ($fy * $this->t[9])))))); |
|
160
|
|
|
} |
|
161
|
|
|
return ((float)$this->offset + ((float)$this->scale * (float)$h)); |
|
162
|
|
|
} |
|
163
|
|
|
} |
|
164
|
|
|
/* |
|
165
|
|
|
$GeoidHeight = new GeoidHeight('../install/tmp/egm96-15.pgm'); |
|
166
|
|
|
$result = $GeoidHeight->get(46.3870,5.2941); |
|
167
|
|
|
var_dump($result); |
|
168
|
|
|
*/ |
|
169
|
|
|
?> |
Ysurac /
FlightAirMap
Loading history...
Adding explicit visibility (
private,protected, orpublic) is generally recommend to communicate to other developers how, and from where this method is intended to be used.