-
Notifications
You must be signed in to change notification settings - Fork 110
/
fillsinks.m
executable file
·160 lines (138 loc) · 4.22 KB
/
fillsinks.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
function DEM = fillsinks(DEM,maxdepth)
%FILLSINKS fill/remove pits, sinks or topographic depressions
%
% Syntax
%
% DEMfs = fillsinks(DEM)
% DEMfs = fillsinks(DEM,maxdepth)
% DEMfs = fillsinks(DEM,sinks)
% DEMfs = fillsinks(DEM,sinks,option)
%
% Description
%
% fillsinks removes topographic depressions in a Digital Elevation
% Model (DEM). Use this function to enable a continuous flow towards
% the DEM edges.
%
% Sinks may, however, be closed basins or dolines and as such they are
% important features of DEMs. In order to account for such sinks,
% fillsinks allows you to specify a maximum depth of sinks, that
% will be filled, or to employ a logical grid (sinks) that is true
% where sinks should remain (minima imposition). Note that for latter
% option, there will be one regional minima for each connected
% component in the sinks grid.
%
% Input
%
% DEM digital elevation model (GRIDobj)
% maxdepth positive scalar with maximum depth of sinks that will be
% filled
% SINKS logical matrix same size as dem with true elements
% referring to sinks (GRIDobj)
%
% Output
%
% DEMfs digital elevation model with filled sinks (GRIDobj)
%
% Example
%
% DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
% DEMf = fillsinks(DEM);
% % display sink depth
% DIFFDEM = DEMf-DEM;
% DIFFDEM.Z(DIFFDEM.Z==0) = nan;
% imageschs(DEM,DIFFDEM.Z);
%
%
% See also: IMFILL, IMRECONSTRUCT, IMIMPOSEMIN, PREPROCESSTOOL
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 28. January, 2013
% Check input arguments
narginchk(1, 2)
if nargin == 1
else
if isscalar(maxdepth) && ~isa(maxdepth,'GRIDobj')
validateattributes(maxdepth,{'numeric'},{'>',0});
md = true;
else
SINKS = maxdepth;
validatealignment(DEM,SINKS);
if isa(SINKS,'GRIDobj')
SINKS = SINKS.Z;
end
md = false;
end
end
dem = DEM.Z;
% START here
% identify nans
Inan = isnan(dem);
% set nans to -inf
dem(Inan) = -inf;
if nargin == 1
% fill depressions using imreconstruct with an 8-neighborhood
marker = -dem;
II = false(size(dem));
II(2:end-1,2:end-1) = true;
marker(II & ~Inan) = -inf;
demfs = -imreconstruct(marker,-dem,8);
elseif nargin==2 && md
% create mask
% complement image
dem = imcomplement(dem);
% create marker
marker = dem;
marker(2:end-1,2:end-1) = -inf;
if any(Inan)
I = (imdilate(Inan,ones(3)) & ~Inan);
marker(I) = dem(I);
end
demfs = imreconstruct(marker,dem);
% difference image between filled and original DEM
D = dem-demfs;
Inan = ~Inan;
Inan = Inan(:);
while any(D(Inan) > maxdepth)
% find maximum difference in each connected sink area
I = D>0;
STATS = regionprops(I,D,'MaxIntensity','PixelIdxList');
for r = 1:numel(STATS)
if STATS(r).MaxIntensity < maxdepth
% do nothing
else
[~,ix] = max(D(STATS(r).PixelIdxList));
ix = ix(1);
marker(STATS(r).PixelIdxList(ix)) = dem(STATS(r).PixelIdxList(ix));
end
end
demfs = imreconstruct(marker,dem);
D = dem-demfs;
end
% complement image again
demfs = imcomplement(demfs);
Inan = reshape(~Inan,DEM.size);
else
% minima imposition
I = true(size(dem));
I(2:end-1,2:end-1) = false;
if any(Inan(:))
I = (imdilate(Inan,ones(3)) & ~Inan) | I;
end
% Refine markers by identifying the lowest value in each sink
STATS = regionprops(SINKS,dem,'MinIntensity','PixelIdxList');
SINKS = false(size(SINKS));
for r = 1:numel(STATS)
ix = find(dem(STATS(r).PixelIdxList)==STATS(r).MinIntensity,1,'first');
SINKS(STATS(r).PixelIdxList(ix)) = true;
end
marker = -inf(size(dem),class(dem));
marker(SINKS) = -dem(SINKS);
marker(I) = -dem(I);
demfs = -imreconstruct(marker,-dem);
end
% nans in the dem are set to nan again
if isfloat(demfs)
demfs(Inan) = nan;
end
DEM.Z = demfs;