-
Notifications
You must be signed in to change notification settings - Fork 110
/
hillshade.m
executable file
·171 lines (142 loc) · 4.71 KB
/
hillshade.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
function OUT2 = hillshade(DEM,varargin)
%HILLSHADE create hillshading from a digital elevation model (GRIDobj)
%
% Syntax
%
% H = hillshade(DEM)
% H = hillshade(DEM,'pn','pv',...)
%
% Description
%
% Hillshading is a very powerful tool for relief depiction.
% hillshade calculates a shaded relief for a digital elevation model
% based on the angle between the surface and the incoming light beams.
% If no output arguments are defined, the hillshade matrix will be
% plotted with a gray colormap. The hillshading algorithm follows the
% logarithmic approach to shaded relief representation of Katzil and
% Doytsher (2003).
%
% Input
%
% DEM Digital elevation model (class: GRIDobj)
%
% Parameter name/value pairs
%
% 'azimuth' azimuth angle, (default=315)
% 'altitude' altitude angle, (default=60)
% 'exaggerate' elevation exaggeration (default=1). Increase to
% pronounce elevation differences in flat terrain
% 'useblockproc' true or {false}: use block processing
% (see function blockproc)
% 'useparallel' true or {false}: use parallel computing toolbox
% 'blocksize' blocksize for blockproc (default: 5000)
% 'method' 'surfnorm' (default) or 'mdow'
%
%
% Output
%
% H shaded relief (ranges between 0 and 1)
%
%
% Example
%
% DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
% hillshade(DEM)
%
% References
%
% Katzil, Y., Doytsher, Y. (2003): A logarithmic and sub-pixel approach
% to shaded relief representation. Computers & Geosciences, 29,
% 1137-1142.
%
% See also: SURFNORM, IMAGESCHS
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 18. August, 2017
% Parse inputs
p = inputParser;
p.StructExpand = true;
p.KeepUnmatched = false;
p.FunctionName = 'hillshade';
addParamValue(p,'azimuth',315,@(x) isscalar(x) && x>= 0 && x<=360);
addParamValue(p,'altitude',60,@(x) isscalar(x) && x>= 0 && x<=90);
addParamValue(p,'exaggerate',1,@(x) isscalar(x) && x>0);
addParamValue(p,'useparallel',true);
addParamValue(p,'blocksize',2000);
addParamValue(p,'useblockproc',true,@(x) isscalar(x));
addParamValue(p,'method','default');
parse(p,varargin{:});
OUT = DEM;
OUT.Z = [];
cs = DEM.cellsize;
azimuth = p.Results.azimuth;
altitude = p.Results.altitude;
exaggerate = p.Results.exaggerate;
method = validatestring(p.Results.method,{'default','surfnorm','mdow'});
% Large matrix support. Break calculations in chunks using blockproc
if numel(DEM.Z)>(10001*10001) && p.Results.useblockproc
blksiz = bestblk(size(DEM.Z),p.Results.blocksize);
padval = 'symmetric';
Z = DEM.Z;
% The anonymous function must be defined as a variable: see bug 1157095
fun = @(x) hsfun(x,cs,azimuth,altitude,exaggerate,method);
HS = blockproc(Z,blksiz,fun,...
'BorderSize',[1 1],...
'padmethod',padval,...
'UseParallel',p.Results.useparallel);
OUT.Z = HS;
else
OUT.Z = hsfun(DEM.Z,cs,azimuth,altitude,exaggerate,method);
end
OUT.name = 'hillshade';
OUT.zunit = '';
if nargout == 0
OUT.Z = uint8(OUT.Z*255);
imagesc(OUT);
colormap(gray)
else
OUT2 = OUT;
end
end
%% Subfunction
function H = hsfun(Z,cs,azimuth,altitude,exaggerate,method)
if isstruct(Z)
Z = Z.data;
end
switch method
case {'default','surfnorm'}
% correct azimuth so that angles go clockwise from top
azid = azimuth-90;
% use radians
altsource = altitude/180*pi;
azisource = azid/180*pi;
% calculate solar vector
[sx,sy,sz] = sph2cart(azisource,altsource,1);
% calculate surface normals
[Nx,Ny,Nz] = surfnorm(Z/cs*exaggerate);
% calculate cos(angle)
% H = [Nx(:) Ny(:) Nz(:)]*[sx;sy;sz];
% % reshape
% H = reshape(H,size(Nx));
H = Nx*sx + Ny*sy + Nz*sz;
% % usual GIS approach
% H = acos(H);
% % force H to range between 0 and 1
% H = H-min(H(:));
% H = H/max(H(:));
case 'mdow'
% correct azimuth so that angles go clockwise from top
azid = [360 315 225 270] - 90;
azisource = azid/180*pi;
altsource = 30/180*pi;
altsource = repmat(altsource,size(azisource));
% calculate solar vector
[sx,sy,sz] = sph2cart(azisource,altsource,1);
% calculate surface normals
[Nx,Ny,Nz] = surfnorm(Z/cs*exaggerate);
H = bsxfun(@times,Nx(:),sx) + bsxfun(@times,Ny(:),sy) + bsxfun(@times,Nz(:),sz);
H = max(H,0);
H = sum(H,2)./3;
H = reshape(H,size(Z));
end
end