1 /**********************************************************
2 * Version $Id: flow_massflux.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_hydrology //
13 // //
14 //-------------------------------------------------------//
15 // //
16 // flow_massflux.cpp //
17 // //
18 // Copyright (C) 2009 by //
19 // Olaf Conrad //
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: oconrad@saga-gis.org //
43 // //
44 // contact: Olaf Conrad //
45 // Institute of Geography //
46 // University of Hamburg //
47 // Bundesstr. 55 //
48 // 20146 Hamburg //
49 // Germany //
50 // //
51 ///////////////////////////////////////////////////////////
52
53 ///////////////////////////////////////////////////////////
54 // //
55 // //
56 // //
57 ///////////////////////////////////////////////////////////
58
59 //---------------------------------------------------------
60 #include "flow_massflux.h"
61
62
63 ///////////////////////////////////////////////////////////
64 // //
65 // //
66 // //
67 ///////////////////////////////////////////////////////////
68
69 //---------------------------------------------------------
70 const int xDir[4] = { 1, 1, 0, 0 };
71 const int yDir[4] = { 1, 0, 0, 1 };
72
73
74 ///////////////////////////////////////////////////////////
75 // //
76 // //
77 // //
78 ///////////////////////////////////////////////////////////
79
80 //---------------------------------------------------------
CFlow_MassFlux(void)81 CFlow_MassFlux::CFlow_MassFlux(void)
82 {
83 Set_Name (_TL("Flow Accumulation (Mass-Flux Method)"));
84
85 Set_Author ("O. Conrad (c) 2009");
86
87 Set_Description (_TW(
88 "The Mass-Flux Method (MFM) for the DEM based calculation of flow accumulation "
89 "as proposed by Gruber and Peckham (2008).\n"
90 "\n"
91 "!!!UNDER DEVELOPMENT!!! To be done: solving the streamline resolution problem"
92 ));
93
94 Add_Reference("Gruber, S. & Peckham, S.", "2008",
95 "Land-Surface Parameters and Objects in Hydrology",
96 "In: Hengl, T. & Reuter, H.I. [Eds.]: Geomorphometry: Concepts, Software, Applications. "
97 "Developments in Soil Science, Elsevier, Bd.33, S.293-308.",
98 SG_T("https://www.researchgate.net/profile/Scott_Peckham2/publication/256829608_Chapter_7_Land-Surface_Parameters_and_Objects_in_Hydrology/links/0c960523c979588e91000000.pdf"),
99 SG_T("ResearchGate")
100 );
101
102 //-----------------------------------------------------
103 Parameters.Add_Grid("",
104 "DEM" , _TL("Elevation"),
105 _TL(""),
106 PARAMETER_INPUT
107 );
108
109 Parameters.Add_Grid("",
110 "AREA" , _TL("Flow Accumulation"),
111 _TL(""),
112 PARAMETER_OUTPUT
113 );
114
115 Parameters.Add_Choice("",
116 "METHOD" , _TL("Flow Split Method"),
117 _TL(""),
118 CSG_String::Format("%s|%s",
119 _TL("flow width (original)"),
120 _TL("cell area")
121 ), 0
122 );
123
124 //-----------------------------------------------------
125 Parameters.Add_Node("",
126 "QUARTERS" , _TL("Create Output of Quarter Cell Information"),
127 _TL("")
128 );
129
130 Parameters.Add_Bool("QUARTERS" , "B_SLOPE" , _TL("Slope" ), _TL(""));
131 Parameters.Add_Grid_Output ("", "G_SLOPE" , _TL("Slope" ), _TL(""));
132
133 Parameters.Add_Bool("QUARTERS" , "B_ASPECT", _TL("Aspect" ), _TL(""));
134 Parameters.Add_Grid_Output ("", "G_ASPECT", _TL("Aspect" ), _TL(""));
135
136 Parameters.Add_Bool("QUARTERS" , "B_AREA" , _TL("Flow Accumulation"), _TL(""));
137 Parameters.Add_Grid_Output ("", "G_AREA" , _TL("Flow Accumulation"), _TL(""));
138
139 Parameters.Add_Bool("QUARTERS" , "B_FLOW" , _TL("Flow Lines" ), _TL(""));
140 Parameters.Add_Shapes_Output("", "G_FLOW" , _TL("Flow Lines" ), _TL(""));
141 }
142
143
144 ///////////////////////////////////////////////////////////
145 // //
146 ///////////////////////////////////////////////////////////
147
148 //---------------------------------------------------------
On_Execute(void)149 bool CFlow_MassFlux::On_Execute(void)
150 {
151 int x, y, i, ix, iy;
152 CSG_Grid *pArea;
153
154 //-----------------------------------------------------
155 m_pDEM = Parameters("DEM" )->asGrid();
156 pArea = Parameters("AREA" )->asGrid();
157 m_Method = Parameters("METHOD")->asInt();
158
159 //-----------------------------------------------------
160 if( 1 )
161 {
162 CSG_Grid_System System(0.5 * Get_Cellsize(), Get_XMin() - 0.25 * Get_Cellsize(), Get_YMin() - 0.25 * Get_Cellsize(), 2 * Get_NX(), 2 * Get_NY());
163
164 m_Area .Create(System, SG_DATATYPE_Float);
165 m_dir .Create(System, SG_DATATYPE_Byte);
166 m_dif .Create(System, SG_DATATYPE_Float);
167
168 m_dir .Assign(-1.0);
169 m_Area .Assign( 0.0);
170 m_Area .Set_NoData_Value(0.0);
171
172 Parameters("G_SLOPE" )->Set_Value(m_pSlope = !Parameters("B_SLOPE" )->asBool() ? NULL : SG_Create_Grid(System, SG_DATATYPE_Float));
173 Parameters("G_ASPECT")->Set_Value(m_pAspect = !Parameters("B_ASPECT")->asBool() ? NULL : SG_Create_Grid(System, SG_DATATYPE_Float));
174 Parameters("G_FLOW" )->Set_Value(m_pFlow = !Parameters("B_FLOW" )->asBool() ? NULL : SG_Create_Shapes(SHAPE_TYPE_Line, _TL("Flow Lines")));
175
176 //-------------------------------------------------
177 // Calculate flow portions...
178 for(y=0; y<Get_NY() && Set_Progress(y); y++)
179 {
180 for(x=0; x<Get_NX(); x++)
181 {
182 for(i=0; i<4; i++)
183 {
184 Set_Flow(x, y, i);
185 }
186 }
187 }
188
189 //-------------------------------------------------
190 // Check for consistent flow directions...
191
192 // still missing...
193
194
195 //-------------------------------------------------
196 // Calculate flow accumulation...
197 for(y=0, iy=0; y<Get_NY() && Set_Progress(y); y++, iy+=2)
198 {
199 for(x=0, ix=0; x<Get_NX(); x++, ix+=2)
200 {
201 for(i=0; i<4; i++)
202 {
203 Get_Area(ix + xDir[i], iy + yDir[i]);
204 }
205 }
206 }
207
208 //-------------------------------------------------
209 // Scale flow accumulation to original cell size...
210 for(y=0, iy=0; y<Get_NY() && Set_Progress(y); y++, iy+=2)
211 {
212 for(x=0, ix=0; x<Get_NX(); x++, ix+=2)
213 {
214 double Area = 0.0, d;
215
216 for(i=0; i<4; i++)
217 {
218 if( (d = m_Area.asDouble(ix + xDir[i], iy + yDir[i])) > 0.0 )
219 {
220 Area += d;
221 }
222 }
223
224 if( Area > 0.0 )
225 {
226 pArea->Set_Value(x, y, Area * m_Area.Get_Cellarea());
227 }
228 else
229 {
230 pArea->Set_NoData(x, y);
231 }
232 }
233 }
234
235 //-------------------------------------------------
236 if( Parameters("B_AREA")->asBool() )
237 {
238 Parameters("G_AREA")->Set_Value(SG_Create_Grid(m_Area));
239 }
240
241 m_Area.Destroy();
242 m_dif .Destroy();
243 m_dir .Destroy();
244
245 DataObject_Set_Colors(pArea, 11, SG_COLORS_WHITE_BLUE);
246
247 return( true );
248 }
249
250 //-----------------------------------------------------
251 return( false );
252 }
253
254
255 ///////////////////////////////////////////////////////////
256 // //
257 ///////////////////////////////////////////////////////////
258
259 //---------------------------------------------------------
Set_Flow(int x,int y,int Direction)260 bool CFlow_MassFlux::Set_Flow(int x, int y, int Direction)
261 {
262 int dir, ix, iy, jx, jy;
263 double dif, Z, A, B;
264
265 if( m_pDEM->is_InGrid(x, y)
266 && m_pDEM->is_InGrid(ix = Get_xTo(2 * Direction , x), iy = Get_yTo(2 * Direction , y))
267 && m_pDEM->is_InGrid(jx = Get_xTo(2 * Direction + 2, x), jy = Get_yTo(2 * Direction + 2, y)) )
268 {
269 Z = m_pDEM->asDouble( x, y);
270 A = (m_pDEM->asDouble(ix, iy) - Z) / Get_Cellsize();
271 B = (m_pDEM->asDouble(jx, jy) - Z) / Get_Cellsize();
272
273 dif = A != 0.0 ? M_PI_180 + atan2(B, A) : (B > 0.0 ? M_PI_270 : (B < 0.0 ? M_PI_090 : -1.0));
274
275 if( dif >= 0.0 )
276 {
277 x = 2 * x + xDir[Direction];
278 y = 2 * y + yDir[Direction];
279
280 dif = fmod(dif + M_PI_090 * Direction, M_PI_360);
281
282 if( m_pFlow )
283 {
284 double dScale = 0.50 * m_dir.Get_System().Get_Cellsize();
285 TSG_Point Point = m_dir.Get_System().Get_Grid_to_World(x, y);
286 CSG_Shape *pLine = m_pFlow->Add_Shape();
287
288 pLine->Add_Point(
289 Point.x - dScale * sin(dif),
290 Point.y - dScale * cos(dif), 0
291 );
292 pLine->Add_Point(Point, 0);
293
294 dScale = 0.20 * m_dir.Get_System().Get_Cellsize();
295 pLine->Add_Point(
296 Point.x - dScale * sin(dif + 25.0 * M_DEG_TO_RAD),
297 Point.y - dScale * cos(dif + 25.0 * M_DEG_TO_RAD), 1
298 );
299 pLine->Add_Point(Point, 1);
300 pLine->Add_Point(
301 Point.x - dScale * sin(dif - 25.0 * M_DEG_TO_RAD),
302 Point.y - dScale * cos(dif - 25.0 * M_DEG_TO_RAD), 1
303 );
304 }
305
306 if( m_pSlope ) m_pSlope ->Set_Value(x, y, atan(sqrt(A*A + B*B)));
307 if( m_pAspect ) m_pAspect->Set_Value(x, y, dif);
308
309 dir = (int)(dif / M_PI_090);
310
311 dif = dif - dir * M_PI_090;
312
313 switch( m_Method )
314 {
315 case 0:
316 dif = cos(dif) / (cos(dif) + sin(dif));
317 break;
318
319 case 1:
320 dif = dif < M_PI_045 ? 1.0 - tan(dif) / 2.0 : tan(M_PI_090 - dif) / 2.0;
321 break;
322 }
323
324 m_dir.Set_Value(x, y, dir * 2);
325 m_dif.Set_Value(x, y, dif);
326
327 return( true );
328 }
329 }
330
331 return( false );
332 }
333
334
335 ///////////////////////////////////////////////////////////
336 // //
337 // //
338 // //
339 ///////////////////////////////////////////////////////////
340
341 //---------------------------------------------------------
Get_Flow(int x,int y,int Direction)342 inline double CFlow_MassFlux::Get_Flow(int x, int y, int Direction)
343 {
344 if( m_dir.is_InGrid(x, y) )
345 {
346 int i = m_dir.asInt(x, y);
347
348 if( Direction == i )
349 {
350 return( m_dif.asDouble(x, y) );
351 }
352
353 if( Direction == (i + 2) % 8 )
354 {
355 return( 1.0 - m_dif.asDouble(x, y) );
356 }
357 }
358
359 return( 0.0 );
360 }
361
362 //---------------------------------------------------------
Get_Area(int x,int y)363 double CFlow_MassFlux::Get_Area(int x, int y)
364 {
365 if( m_Area.is_NoData(x, y) ) // cell has not been processed yet...
366 {
367 int i, ix, iy;
368 double d;
369
370 m_Area.Set_Value(x, y, 1.0); // add this cell's contribution...
371
372 for(i=0; i<8; i+=2)
373 {
374 ix = Get_xFrom(i, x);
375 iy = Get_yFrom(i, y);
376
377 d = Get_Flow(ix, iy, i);
378
379 if( d > 0.0 ) // which portion drains ith neighbour into this cell???
380 {
381 m_Area.Add_Value(x, y, d * Get_Area(ix, iy)); // then recursive call of this function...
382 }
383 }
384 }
385
386 return( m_Area.asDouble(x, y) ); // return this cell's area...
387 }
388
389
390 ///////////////////////////////////////////////////////////
391 // //
392 // //
393 // //
394 ///////////////////////////////////////////////////////////
395
396 //---------------------------------------------------------
397