1 /**********************************************************
2 * Version $Id: FillSinks_WL_XXL.cpp 1921 2014-01-09 10:24:11Z oconrad $
3 *********************************************************/
4
5 ///////////////////////////////////////////////////////////
6 // //
7 // SAGA //
8 // //
9 // System for Automated Geoscientific Analyses //
10 // //
11 // Tool Library //
12 // ta_preprocessor //
13 // //
14 //-------------------------------------------------------//
15 // //
16 // FillSinks_WL.cpp //
17 // //
18 // Copyright (C) 2007 by //
19 // Volker Wichmann //
20 // //
21 //-------------------------------------------------------//
22 // //
23 // This file is part of 'SAGA - System for Automated //
24 // Geoscientific Analyses'. SAGA is free software; you //
25 // can redistribute it and/or modify it under the terms //
26 // of the GNU General Public License as published by the //
27 // Free Software Foundation, either version 2 of the //
28 // License, or (at your option) any later version. //
29 // //
30 // SAGA is distributed in the hope that it will be //
31 // useful, but WITHOUT ANY WARRANTY; without even the //
32 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
33 // PARTICULAR PURPOSE. See the GNU General Public //
34 // License for more details. //
35 // //
36 // You should have received a copy of the GNU General //
37 // Public License along with this program; if not, see //
38 // <http://www.gnu.org/licenses/>. //
39 // //
40 //-------------------------------------------------------//
41 // //
42 // e-mail: reklovw@web.de //
43 // //
44 ///////////////////////////////////////////////////////////
45
46
47 ///////////////////////////////////////////////////////////
48 // //
49 // //
50 // //
51 ///////////////////////////////////////////////////////////
52
53 //---------------------------------------------------------
54 #include "FillSinks_WL.h"
55
56
57 ///////////////////////////////////////////////////////////
58 // //
59 // Construction/Destruction //
60 // //
61 ///////////////////////////////////////////////////////////
62
CFillSinks_WL_XXL(void)63 CFillSinks_WL_XXL::CFillSinks_WL_XXL(void)
64 {
65
66 Set_Name (_TL("Fill Sinks XXL (Wang & Liu)"));
67 Set_Author (_TL("Copyrights (c) 2007 by Volker Wichmann"));
68 Set_Description (_TW(
69 "This tool uses an algorithm proposed by Wang & Liu to identify and fill surface depressions in "
70 "digital elevation models.\n"
71 "The method was enhanced to allow the creation of hydrologic sound elevation models, i.e. not only to "
72 "fill the depression(s) but also to preserve a downward slope along the flow path. If desired, this is accomplished "
73 "by preserving a minimum slope gradient (and thus elevation difference) between cells.\n"
74 "This version of the tool is designed to work on large data sets (e.g. LIDAR data), with smaller "
75 "datasets you might like to check out the fully featured standard version of the tool.\n\n\n"
76 "References:\n"
77 "Wang, L. & H. Liu (2006): An efficient method for identifying and filling surface depressions in "
78 "digital elevation models for hydrologic analysis and modelling. International Journal of Geographical "
79 "Information Science, Vol. 20, No. 2: 193-213.\n"
80 ));
81
82
83 Parameters.Add_Grid(
84 NULL, "ELEV", _TL("DEM"),
85 _TL("Digital elevation model"),
86 PARAMETER_INPUT
87 );
88
89 Parameters.Add_Grid(
90 NULL, "FILLED", _TL("Filled DEM"),
91 _TL("Depression-free digital elevation model"),
92 PARAMETER_OUTPUT
93 );
94
95 Parameters.Add_Value(
96 NULL, "MINSLOPE", _TL("Minimum Slope [Degree]"),
97 _TL("Minimum slope gradient to preserve from cell to cell; with a value of zero sinks are filled up to the spill elevation (which results in flat areas). Unit [Degree]"),
98 PARAMETER_TYPE_Double, 0.1, 0.0, true
99 );
100
101 }
102
103 //---------------------------------------------------------
~CFillSinks_WL_XXL(void)104 CFillSinks_WL_XXL::~CFillSinks_WL_XXL(void)
105 {}
106
107
108 ///////////////////////////////////////////////////////////
109 // //
110 // //
111 // //
112 ///////////////////////////////////////////////////////////
113
114
115
On_Execute(void)116 bool CFillSinks_WL_XXL::On_Execute(void)
117 {
118 CSG_Grid *pElev, *pFilled;
119 PriorityQ theQueue;
120 CFillSinks_WL_Node tempNode;
121
122 int x, y, ix, iy, i;
123 double z, iz, progress;
124 double minslope, mindiff[8];
125 bool preserve;
126
127
128 pElev = Parameters("ELEV")->asGrid();
129 pFilled = Parameters("FILLED")->asGrid();
130 minslope = Parameters("MINSLOPE")->asDouble();
131
132 pFilled->Fmt_Name("%s [%s]", pElev->Get_Name(), _TL("no sinks"));
133
134
135 if( minslope > 0.0 )
136 {
137 minslope = tan(minslope * M_DEG_TO_RAD);
138 for(i=0; i<8; i++)
139 mindiff[i] = minslope * Get_Length(i);
140 preserve = true;
141 }
142 else
143 preserve = false;
144
145
146 pFilled->Assign_NoData();
147
148
149 // fill priority queue with boundary cells, i.e. seed cells
150 for(y=0; y<Get_NY(); y++)
151 {
152 for(x=0; x<Get_NX(); x++)
153 {
154 if( !pElev->is_NoData(x, y) )
155 {
156 for(i=0; i<8; i++)
157 {
158 ix = Get_xTo(i, x);
159 iy = Get_yTo(i, y);
160 if( !is_InGrid(ix, iy) || pElev->is_NoData(ix, iy) )
161 {
162 z = pElev->asDouble(x, y);
163
164 tempNode.x = x;
165 tempNode.y = y;
166 tempNode.spill = z;
167 theQueue.push( tempNode );
168
169 pFilled->Set_Value(x, y, z);
170 break;
171 }
172 }
173 }
174 }
175 }
176
177
178 // process queue
179 progress = 0.0;
180
181 while( !theQueue.empty() )
182 {
183 PriorityQ::value_type tempNode = theQueue.top();
184
185 x = tempNode.x;
186 y = tempNode.y;
187 theQueue.pop();
188
189 z = pFilled->asDouble(x, y);
190
191 for(i=0; i<8; i++)
192 {
193 ix = Get_xTo(i, x);
194 iy = Get_yTo(i, y);
195 if( is_InGrid(ix, iy) && !pElev->is_NoData(ix, iy) && pFilled->is_NoData(ix, iy) )
196 {
197 iz = pElev->asDouble(ix, iy);
198
199 if( preserve )
200 {
201 if( iz < (z + mindiff[i]) )
202 iz = z + mindiff[i];
203 }
204 else if( iz < z )
205 iz = z;
206
207 tempNode.x = ix;
208 tempNode.y = iy;
209 tempNode.spill = iz;
210 theQueue.push( tempNode );
211
212 pFilled->Set_Value(ix, iy, iz);
213 }
214 }
215
216 progress += 1.0;
217 if( ((int)progress) % 10000 == 0 )
218 Set_Progress(progress, (double)pElev->Get_NCells());
219 //DataObject_Update(pFilled, pElev->Get_ZMin(), pElev->Get_ZMax(), true);
220
221 }
222
223 return (true);
224 }
225
226