forked from geographiclib/geographiclib-octave
-
Notifications
You must be signed in to change notification settings - Fork 0
/
geoid_height.m
176 lines (172 loc) · 7.1 KB
/
geoid_height.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
function N = geoid_height(lat, lon, geoidname, geoiddir)
%GEOID_HEIGHT Compute the height of the geoid above the ellipsoid
%
% N = GEOID_HEIGHT(lat, lon)
% N = GEOID_HEIGHT(lat, lon, geoidname)
% N = GEOID_HEIGHT(lat, lon, geoidname, geoiddir)
% GEOID_HEIGHT([])
% N = GEOID_HEIGHT(lat, lon, geoid)
%
% computes the height, N, of the geoid above the WGS84 ellipsoid. lat
% and lon are the latitude and longitude in degrees; these can be scalars
% or arrays of the same size. N is in meters.
%
% The height of the geoid above the ellipsoid, N, is sometimes called the
% geoid undulation. It can be used to convert a height above the
% ellipsoid, h, to the corresponding height above the geoid (the
% orthometric height, roughly the height above mean sea level), H, using
% the relations
%
% h = N + H; H = -N + h.
%
% The possible geoids are
%
% egm84-30 egm84-15
% egm96-15 egm96-5
% egm2008-5 egm2008-2_5 egm2008-1
%
% The first part of the name is the geoid model. The second part gives
% the resolution of the gridded data (in arc-seconds). Thus egm2008-2_5
% is the egm2008 geoid model at a resolution of 2.5".
%
% By default the egm96-5 geoid is used. This can be overridden by
% specifying geoidname. The geoiddir argument overrides the default
% directory for the model. See the documentation on geoid_load for how
% these arguments are interpreted.
%
% When geoid_height is invoked with a particular geoidname, the geoid
% data is loaded from disk and cached. A subsequent invocation of
% geoid_height with the same geoidname uses the cached data. Use
% GEOID_HEIGHT([]) to clear this cached data.
%
% Alternatively, you can load a geoid with geoid_load and use the
% returned structure as the third argument. This use does not change the
% cached data.
%
% In order to use this routine with Octave, Octave needs to have been
% compiled with a version of GraphicsMagick which supports 16-bit images.
% Also, the first time you uses this routine, you may receive a warning
% message "your version of GraphicsMagick limits images to 16 bits per
% pixel"; this can safely be ignored.
%
% Information on downloading and installing the data for the supported
% geoid models is available at
%
% https://geographiclib.sourceforge.io/C++/doc/geoid.html#geoidinst
%
% GEOID_HEIGHT uses cubic interpolation on gridded data that has been
% quantized at a resolution of 3mm.
%
% See also GEOID_LOAD.
% Copyright (c) Charles Karney (2015-2022) <[email protected]>.
persistent saved_geoid
if nargin == 1 && isempty(lat)
saved_geoid = [];
return
end
narginchk(2, 4)
if nargin == 3 && isstruct(geoidname)
N = geoid_height_int(lat, lon, geoidname, true);
else
if nargin < 3
geoidname = '';
end
if nargin < 4
geoiddir = '';
end
geoidfile = geoid_file(geoidname, geoiddir);
if ~(isstruct(saved_geoid) && strcmp(saved_geoid.file, geoidfile))
saved_geoid = geoid_load_file(geoidfile);
end
N = geoid_height_int(lat, lon, saved_geoid, true);
end
end
function N = geoid_height_int(lat, lon, geoid, cubic)
persistent c0 c3 c0n c3n c0s c3s
if isempty(c3s)
c0 = 240;
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];
c0n = 372;
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];
c0s = 372;
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];
end
try
s = size(lat + lon);
catch
error('lat, lon have incompatible sizes')
end
num = prod(s); Z = -zeros(num,1);
lat = lat(:) + Z; lon = lon(:) + Z;
h = geoid.h; w = geoid.w;
% lat is in [0, h]
flat = min(max((90 - lat) * (h - 1) / 180, 0), (h - 1));
% lon is in [0, w)
flon = mod(lon * w / 360, w);
flon(isnan(flon)) = 0;
ilat = min(floor(flat), h - 2);
ilon = floor(flon);
flat = flat - ilat; flon = flon - ilon;
if ~cubic
ind = imgind(ilon + [0,0,1,1], ilat + [0,1,0,1], w, h);
hf = double(geoid.im(ind));
N = (1 - flon) .* ((1 - flat) .* hf(:,1) + flat .* hf(:,2)) + ...
flon .* ((1 - flat) .* hf(:,3) + flat .* hf(:,4));
else
ind = imgind(repmat(ilon, 1, 12) + ...
repmat([ 0, 1,-1, 0, 1, 2,-1, 0, 1, 2, 0, 1], num, 1), ...
repmat(ilat, 1, 12) + ...
repmat([-1,-1, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2], num, 1), ...
w, h);
hf = double(geoid.im(ind));
hfx = hf * c3 / c0;
hfx(ilat == 0,:) = hf(ilat == 0,:) * c3n / c0n;
hfx(ilat == h-2,:) = hf(ilat == h-2,:) * c3s / c0s;
N = sum(hfx .* [Z+1, flon, flat, flon.^2, flon.*flat, flat.^2, ...
flon.^3, flon.^2.*flat, flon.*flat.^2, flat.^3], ...
2);
end
N = geoid.offset + geoid.scale * N;
N(~(abs(lat) <= 90 & isfinite(lon))) = nan;
N = reshape(N, s);
end
function ind = imgind(ix, iy, w, h)
% return 1-based 1d index to w*h array for 0-based 2d indices (ix,iy)
c = iy < 0; iy(c) = - iy(c); ix(c) = ix(c) + w/2;
c = iy >= h; iy(c) = 2 * (h-1) - iy(c); ix(c) = ix(c) + w/2;
ix = mod(ix, w);
ind = 1 + iy + ix * h;
end