1 /**********************************************************
2 * Version $Id$
3 *********************************************************/
4
5 ///////////////////////////////////////////////////////////
6 // //
7 // SAGA //
8 // //
9 // System for Automated Geoscientific Analyses //
10 // //
11 // Tool Library //
12 // ta_channels //
13 // //
14 //-------------------------------------------------------//
15 // //
16 // ChannelNetwork_Distance.cpp //
17 // //
18 // Copyright (C) 2003 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 Goettingen //
47 // Goldschmidtstr. 5 //
48 // 37077 Goettingen //
49 // Germany //
50 // //
51 ///////////////////////////////////////////////////////////
52
53 //---------------------------------------------------------
54
55
56 ///////////////////////////////////////////////////////////
57 // //
58 // //
59 // //
60 ///////////////////////////////////////////////////////////
61
62 //---------------------------------------------------------
63 #include "ChannelNetwork_Distance.h"
64
65
66 ///////////////////////////////////////////////////////////
67 // //
68 // //
69 // //
70 ///////////////////////////////////////////////////////////
71
72 //---------------------------------------------------------
CChannelNetwork_Distance(void)73 CChannelNetwork_Distance::CChannelNetwork_Distance(void)
74 {
75 //-----------------------------------------------------
76 Set_Name (_TL("Overland Flow Distance to Channel Network"));
77
78 Set_Author ("O.Conrad (c) 2001-14");
79
80 Set_Description (_TW(
81 "This tool calculates overland flow distances to a channel network "
82 "based on gridded digital elevation data and channel network information. "
83 "The flow algorithm may be either Deterministic 8 (O'Callaghan & Mark 1984) "
84 "or Multiple Flow Direction (Freeman 1991). Sediment Delivery Rates (SDR) "
85 "according to Ali & De Boer (2010) can be computed optionally. "
86 ));
87
88 Add_Reference(
89 "Ali, K. F., De Boer, D. H.", "2010",
90 "Spatially distributed erosion and sediment yield modeling in the upper Indus River basin",
91 "Water Resources Research, 46(8), W08504. doi:10.1029/2009WR008762"
92 );
93
94 Add_Reference(
95 "Freeman, G.T.", "1991",
96 "Calculating catchment area with divergent flow based on a regular grid",
97 "Computers and Geosciences, 17:413-22."
98 );
99
100 Add_Reference(
101 "O'Callaghan, J.F., Mark, D.M.", "1984",
102 "The extraction of drainage networks from digital elevation data",
103 "Computer Vision, Graphics and Image Processing, 28:323-344."
104 );
105
106 Add_Reference(
107 "Nobre, A.D., Cuartas, L.A., Hodnett, M., Renno, C.D., Rodrigues, G., Silveira, A., Waterloo, M., Saleska S.", "2011",
108 "Height Above the Nearest Drainage - a hydrologically relevant new terrain model",
109 "Journal of Hydrology, Vol. 404, Issues 1-2, pp. 13-29, ISSN 0022-1694, 10.1016/j.jhydrol.2011.03.051.",
110 SG_T("http://www.sciencedirect.com/science/article/pii/S0022169411002599")
111 );
112
113 //-----------------------------------------------------
114 Parameters.Add_Grid("",
115 "ELEVATION" , _TL("Elevation"),
116 _TL("A grid that contains elevation data."),
117 PARAMETER_INPUT
118 );
119
120 Parameters.Add_Grid("",
121 "CHANNELS" , _TL("Channel Network"),
122 _TW("A grid providing information about the channel network. It is assumed that no-data cells are not part "
123 "of the channel network. Vice versa all others cells are recognised as channel network members."),
124 PARAMETER_INPUT
125 );
126
127 Parameters.Add_Grid("",
128 "ROUTE" , _TL("Preferred Routing"),
129 _TL("Downhill flow is bound to preferred routing cells, where these are not no-data. Helps to model e.g. small ditches, that are not well represented in the elevation data."),
130 PARAMETER_INPUT_OPTIONAL
131 );
132
133 //-----------------------------------------------------
134 Parameters.Add_Grid("",
135 "DISTANCE" , _TL("Overland Flow Distance"),
136 _TW("The overland flow distance in map units. "
137 "It is assumed that the (vertical) elevation data use the same units as the (horizontal) grid coordinates."),
138 PARAMETER_OUTPUT
139 );
140
141 Parameters.Add_Grid("",
142 "DISTVERT" , _TL("Vertical Overland Flow Distance"),
143 _TL("This is the vertical component of the overland flow"),
144 PARAMETER_OUTPUT
145 );
146
147 Parameters.Add_Grid("",
148 "DISTHORZ" , _TL("Horizontal Overland Flow Distance"),
149 _TL("This is the horizontal component of the overland flow"),
150 PARAMETER_OUTPUT
151 );
152
153 //-----------------------------------------------------
154 Parameters.Add_Grid("",
155 "TIME" , _TL("Flow Travel Time"),
156 _TL("flow travel time to channel expressed in hours based on Manning's Equation"),
157 PARAMETER_OUTPUT_OPTIONAL
158 );
159
160 Parameters.Add_Grid("",
161 "SDR" , _TL("Sediment Yield Delivery Ratio"),
162 _TL("This is the horizontal component of the overland flow"),
163 PARAMETER_OUTPUT_OPTIONAL
164 );
165
166 //-----------------------------------------------------
167 Parameters.Add_Grid("",
168 "FIELDS" , _TL("Fields"),
169 _TL("If set, output is given about the number of fields a flow path visits downhill. For D8 only."),
170 PARAMETER_INPUT_OPTIONAL
171 );
172
173 Parameters.Add_Grid("",
174 "PASSES" , _TL("Fields Visited"),
175 _TL("Number of fields a flow path visits downhill starting at a cell. For D8 only."),
176 PARAMETER_OUTPUT, true, SG_DATATYPE_Short
177 );
178
179 //-----------------------------------------------------
180 Parameters.Add_Choice("",
181 "METHOD" , _TL("Flow Algorithm"),
182 _TL("Choose a flow routing algorithm that shall be used for the overland flow distance calculation:\n- D8\n- MFD"),
183 CSG_String::Format("%s|%s|",
184 _TL("D8"),
185 _TL("MFD")
186 ), 1
187 );
188
189 Parameters.Add_Bool("",
190 "BOUNDARY" , _TL("Boundary Cells"),
191 _TL("Take cells at the boundary of the DEM as channel."),
192 false
193 );
194
195 Parameters.Add_Double("SDR",
196 "FLOW_B" , _TL("Beta"),
197 _TL("catchment specific parameter for sediment delivery ratio calculation"),
198 1.0, 0.0, true
199 );
200
201 Parameters.Add_Grid_or_Const("",
202 "FLOW_K" , _TL("Manning-Strickler Coefficient"),
203 _TL("Manning-Strickler coefficient for flow travel time estimation (reciprocal of Manning's Roughness Coefficient)"),
204 20.0, 0.0, true
205 );
206
207 Parameters.Add_Grid_or_Const("",
208 "FLOW_R" , _TL("Flow Depth"),
209 _TL("flow depth [m] for flow travel time estimation"),
210 0.05, 0.0, true
211 );
212 }
213
214
215 ///////////////////////////////////////////////////////////
216 // //
217 ///////////////////////////////////////////////////////////
218
219 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)220 int CChannelNetwork_Distance::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
221 {
222 if( pParameter->Cmp_Identifier("METHOD") )
223 {
224 pParameters->Set_Enabled("FIELDS", pParameter->asInt() == 0);
225 pParameters->Set_Enabled("PASSES", pParameter->asInt() == 0 && (*pParameters)("FIELDS")->asPointer() != NULL);
226 }
227
228 if( pParameter->Cmp_Identifier("FIELDS") )
229 {
230 pParameters->Set_Enabled("PASSES", pParameter->is_Enabled() && pParameter->asPointer() != NULL);
231 }
232
233 if( pParameter->Cmp_Identifier("TIME") )
234 {
235 pParameters->Set_Enabled("FLOW_K", pParameter->asDataObject() != NULL);
236 pParameters->Set_Enabled("FLOW_R", pParameter->asDataObject() != NULL);
237 pParameters->Set_Enabled("SDR" , pParameter->asDataObject() != NULL);
238 }
239
240 if( pParameter->Cmp_Identifier("SDR") )
241 {
242 pParameters->Set_Enabled("FLOW_B", pParameter->asDataObject() != NULL);
243 }
244
245 return( CSG_Tool_Grid::On_Parameters_Enable(pParameters, pParameter) );
246 }
247
248
249 ///////////////////////////////////////////////////////////
250 // //
251 ///////////////////////////////////////////////////////////
252
253 //---------------------------------------------------------
On_Execute(void)254 bool CChannelNetwork_Distance::On_Execute(void)
255 {
256 //-----------------------------------------------------
257 m_pDEM = Parameters("ELEVATION")->asGrid();
258 m_pChannels = Parameters("CHANNELS" )->asGrid();
259 m_pRoute = Parameters("ROUTE" )->asGrid();
260
261 m_pDistance = Parameters("DISTANCE" )->asGrid();
262 m_pDistVert = Parameters("DISTVERT" )->asGrid();
263 m_pDistHorz = Parameters("DISTHORZ" )->asGrid();
264
265 m_pTime = Parameters("TIME" )->asGrid();
266 m_pSDR = Parameters("SDR" )->asGrid();
267
268 m_Flow_B = Parameters("FLOW_B" )->asDouble();
269 m_Flow_K = Parameters("FLOW_K" )->asDouble();
270 m_Flow_R = Parameters("FLOW_R" )->asDouble();
271 m_pFlow_K = Parameters("FLOW_K" )->asGrid();
272 m_pFlow_R = Parameters("FLOW_R" )->asGrid();
273
274 m_pFields = Parameters("FIELDS" )->asGrid();
275 m_pPasses = Parameters("PASSES" )->asGrid();
276
277 //-----------------------------------------------------
278 m_pDistance->Assign_NoData();
279 m_pDistVert->Assign_NoData();
280 m_pDistHorz->Assign_NoData();
281
282 if( m_pTime )
283 {
284 m_pTime->Assign_NoData();
285 m_pTime->Set_Unit("hours");
286 }
287
288 if( m_pSDR )
289 {
290 m_pSDR ->Assign_NoData();
291 }
292
293 if( m_pFields && m_pPasses )
294 {
295 m_pPasses->Set_NoData_Value(-1.0);
296 m_pPasses->Assign_NoData();
297 }
298
299 if( !m_pDEM->Set_Index() )
300 {
301 Error_Set(_TL("index creation failed"));
302
303 return( false );
304 }
305
306 //-----------------------------------------------------
307 int Method = Parameters("METHOD")->asInt();
308 bool bBoundary = Parameters("BOUNDARY")->asBool();
309
310 for(sLong n=0; n<Get_NCells() && Set_Progress_NCells(n); n++)
311 {
312 int x, y;
313
314 if( m_pDEM->Get_Sorted(n, x, y, false) ) // ascending, only valid dem cells
315 {
316 if( is_Channel(x, y, bBoundary) )
317 {
318 m_pDistance->Set_Value(x, y, 0.0);
319 m_pDistVert->Set_Value(x, y, 0.0);
320 m_pDistHorz->Set_Value(x, y, 0.0);
321
322 if( m_pTime ) m_pTime ->Set_Value(x, y, 0.0);
323 if( m_pFields ) m_pPasses->Set_Value(x, y, 0.0);
324 }
325 else // if( m_pDEM->is_InGrid(x, y) )
326 {
327 switch( Method )
328 {
329 default: Set_D8 (x, y); break;
330 case 1: Set_MFD(x, y); break;
331 }
332
333 if( m_pSDR && !m_pTime->is_NoData(x, y) )
334 {
335 m_pSDR->Set_Value(x, y, exp(-m_Flow_B * m_pTime->asDouble(x, y)));
336 }
337 }
338 }
339 }
340
341 //-----------------------------------------------------
342 return( true );
343 }
344
345
346 ///////////////////////////////////////////////////////////
347 // //
348 ///////////////////////////////////////////////////////////
349
350 //---------------------------------------------------------
is_Channel(int x,int y,bool bBoundary)351 inline bool CChannelNetwork_Distance::is_Channel(int x, int y, bool bBoundary)
352 {
353 if( !m_pChannels->is_NoData(x, y) )
354 {
355 return( true );
356 }
357
358 if( bBoundary )
359 {
360 for(int i=0; i<8; i++)
361 {
362 if( !m_pDEM->is_InGrid(Get_xTo(i, x), Get_yTo(i, y)) )
363 {
364 return( true );
365 }
366 }
367 }
368
369 return( false );
370 }
371
372
373 ///////////////////////////////////////////////////////////
374 // //
375 ///////////////////////////////////////////////////////////
376
377 //---------------------------------------------------------
Get_Travel_Time(int x,int y,int Direction)378 inline double CChannelNetwork_Distance::Get_Travel_Time(int x, int y, int Direction)
379 {
380 int ix = Get_xTo(Direction, x), iy = Get_yTo(Direction, y);
381
382 double v, k, R, dz, dx;
383
384 dz = m_pDEM->is_InGrid(ix, iy) ? m_pDEM->asDouble(x, y) - m_pDEM->asDouble(ix, iy) : 0.1;
385 dx = Get_Length(Direction);
386
387 // k = m_pFlow_K && !m_pFlow_K->is_NoData(x, y) ? m_pFlow_K->asDouble(x, y) : m_Flow_K;
388 k = !m_pFlow_K || (m_pFlow_K->is_NoData(x, y) && m_pFlow_K->is_NoData(ix, iy)) ? m_Flow_K
389 : m_pFlow_K->is_NoData( x, y) ? m_pFlow_K->asDouble(ix, iy)
390 : m_pFlow_K->is_NoData(ix, iy) ? m_pFlow_K->asDouble( x, y)
391 : (m_pFlow_K->asDouble(x, y) + m_pFlow_K->asDouble(ix, iy)) / 2.0;
392
393 // R = m_pFlow_R && !m_pFlow_R->is_NoData(x, y) ? m_pFlow_R->asDouble(x, y) : m_Flow_R;
394 R = !m_pFlow_R || (m_pFlow_R->is_NoData(x, y) && m_pFlow_R->is_NoData(ix, iy)) ? m_Flow_R
395 : m_pFlow_R->is_NoData( x, y) ? m_pFlow_R->asDouble(ix, iy)
396 : m_pFlow_R->is_NoData(ix, iy) ? m_pFlow_R->asDouble( x, y)
397 : (m_pFlow_R->asDouble(x, y) + m_pFlow_R->asDouble(ix, iy)) / 2.0;
398
399 v = k * pow(R, 2.0 / 3.0) * sqrt(dz / dx); // [m / s], simplified Manning equation
400
401 return( dx / (v * 3600.0) ); // return travel time in hours
402 }
403
404
405 ///////////////////////////////////////////////////////////
406 // //
407 ///////////////////////////////////////////////////////////
408
409 //---------------------------------------------------------
Get_D8(int x,int y,int & Direction)410 bool CChannelNetwork_Distance::Get_D8(int x, int y, int &Direction)
411 {
412 double dz, dzMax = 0.0, z = m_pDEM->asDouble(x, y);
413
414 Direction = -1;
415
416 if( m_pRoute )
417 {
418 for(int i=0; i<8; i++)
419 {
420 int ix = Get_xTo(i, x), iy = Get_yTo(i, y);
421
422 if( m_pDEM->is_InGrid(ix, iy) && !m_pRoute->is_NoData(ix, iy) && (dz = (z - m_pDEM->asDouble(ix, iy)) / Get_Length(i)) > dzMax )
423 {
424 Direction = i; dzMax = dz;
425 }
426 }
427
428 if( Direction >= 0 )
429 {
430 return( true );
431 }
432 }
433
434 for(int i=0; i<8; i++)
435 {
436 int ix = Get_xTo(i, x), iy = Get_yTo(i, y);
437
438 if( m_pDEM->is_InGrid(ix, iy) && !m_pDistance->is_NoData(ix, iy) && (dz = (z - m_pDEM->asDouble(ix, iy)) / Get_Length(i)) > dzMax )
439 {
440 Direction = i; dzMax = dz;
441 }
442 }
443
444 return( Direction >= 0 );
445 }
446
447 //---------------------------------------------------------
Set_D8(int x,int y)448 bool CChannelNetwork_Distance::Set_D8(int x, int y)
449 {
450 int Direction;
451
452 if( !Get_D8(x, y, Direction) )
453 {
454 return( false );
455 }
456
457 int ix = Get_xTo(Direction, x), iy = Get_yTo(Direction, y);
458
459 double dz = m_pDEM->asDouble(x, y) - m_pDEM->asDouble(ix, iy);
460 double dx = Get_Length(Direction);
461
462 m_pDistance->Set_Value(x, y, m_pDistance->asDouble(ix, iy) + sqrt(dz*dz + dx*dx));
463 m_pDistVert->Set_Value(x, y, m_pDistVert->asDouble(ix, iy) + dz);
464 m_pDistHorz->Set_Value(x, y, m_pDistHorz->asDouble(ix, iy) + dx);
465
466 if( m_pTime )
467 {
468 double dt = Get_Travel_Time(x, y, Direction);
469
470 m_pTime->Set_Value(x, y, m_pTime->asDouble(ix, iy) + dt);
471 }
472
473 if( m_pFields )
474 {
475 int Crossed = m_pFields->asDouble(ix, iy) == m_pFields->asDouble(x, y) ? 0 : 1;
476
477 m_pPasses->Set_Value(x, y, m_pPasses->asInt(ix, iy) + Crossed);
478 }
479
480 return( true );
481 }
482
483
484 ///////////////////////////////////////////////////////////
485 // //
486 ///////////////////////////////////////////////////////////
487
488 //---------------------------------------------------------
Get_MFD(int x,int y,CSG_Vector & Flow)489 bool CChannelNetwork_Distance::Get_MFD(int x, int y, CSG_Vector &Flow)
490 {
491 const double MFD_Convergence = 1.1;
492
493 double dz, Sum = 0.0, z = m_pDEM->asDouble(x, y);
494
495 if( m_pRoute )
496 {
497 for(int i=0; i<8; i++)
498 {
499 int ix = Get_xTo(i, x), iy = Get_yTo(i, y);
500
501 if( m_pDEM->is_InGrid(ix, iy) && !m_pRoute->is_NoData(ix, iy) && (dz = z - m_pDEM->asDouble(ix, iy)) > 0.0 )
502 {
503 Sum += (Flow[i] = pow(dz / Get_Length(i), MFD_Convergence));
504 }
505 }
506
507 if( Sum > 0.0 )
508 {
509 Flow *= 1. / Sum;
510
511 return( true );
512 }
513 }
514
515 for(int i=0; i<8; i++)
516 {
517 int ix = Get_xTo(i, x), iy = Get_yTo(i, y);
518
519 if( m_pDEM->is_InGrid(ix, iy) && !m_pDistance->is_NoData(ix, iy) && (dz = z - m_pDEM->asDouble(ix, iy)) > 0.0 )
520 {
521 Sum += (Flow[i] = pow(dz / Get_Length(i), MFD_Convergence));
522 }
523 }
524
525 if( Sum > 0.0 )
526 {
527 Flow *= 1. / Sum;
528
529 return( true );
530 }
531
532 return( false );
533 }
534
535 //---------------------------------------------------------
Set_MFD(int x,int y)536 bool CChannelNetwork_Distance::Set_MFD(int x, int y)
537 {
538 CSG_Vector Flow(8);
539
540 if( !Get_MFD(x, y, Flow) )
541 {
542 return( false );
543 }
544
545 double Distance = 0.0, DistVert = 0.0, DistHorz = 0.0, Time = 0.0, SDR = 0.0;
546
547 double z = m_pDEM->asDouble(x, y);
548
549 for(int i=0; i<8; i++)
550 {
551 if( Flow[i] > 0.0 )
552 {
553 int ix = Get_xTo(i, x), iy = Get_yTo(i, y);
554
555 double dz = z - m_pDEM->asDouble(ix, iy);
556 double dx = Get_Length(i);
557
558 Distance += Flow[i] * (m_pDistance->asDouble(ix, iy) + sqrt(dz*dz + dx*dx));
559 DistVert += Flow[i] * (m_pDistVert->asDouble(ix, iy) + dz);
560 DistHorz += Flow[i] * (m_pDistHorz->asDouble(ix, iy) + dx);
561
562 if( m_pTime )
563 {
564 double dt = Get_Travel_Time(x, y, i);
565
566 Time += Flow[i] * (m_pTime->asDouble(ix, iy) + dt);
567 }
568 }
569 }
570
571 if( Distance > 0.0 )
572 {
573 m_pDistance->Set_Value(x, y, Distance);
574 m_pDistVert->Set_Value(x, y, DistVert);
575 m_pDistHorz->Set_Value(x, y, DistHorz);
576
577 if( m_pTime )
578 {
579 m_pTime->Set_Value(x, y, Time);
580 }
581 }
582
583 return( true );
584 }
585
586
587 ///////////////////////////////////////////////////////////
588 // //
589 // //
590 // //
591 ///////////////////////////////////////////////////////////
592
593 //---------------------------------------------------------
594