1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Tool Library //
9 // ta_morphometry //
10 // //
11 //-------------------------------------------------------//
12 // //
13 // Morphometry.cpp //
14 // //
15 // Copyright (C) 2003 by //
16 // Olaf Conrad //
17 // //
18 //-------------------------------------------------------//
19 // //
20 // This file is part of 'SAGA - System for Automated //
21 // Geoscientific Analyses'. SAGA is free software; you //
22 // can redistribute it and/or modify it under the terms //
23 // of the GNU General Public License as published by the //
24 // Free Software Foundation, either version 2 of the //
25 // License, or (at your option) any later version. //
26 // //
27 // SAGA is distributed in the hope that it will be //
28 // useful, but WITHOUT ANY WARRANTY; without even the //
29 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
30 // PARTICULAR PURPOSE. See the GNU General Public //
31 // License for more details. //
32 // //
33 // You should have received a copy of the GNU General //
34 // Public License along with this program; if not, see //
35 // <http://www.gnu.org/licenses/>. //
36 // //
37 //-------------------------------------------------------//
38 // //
39 // e-mail: oconrad@saga-gis.org //
40 // //
41 // contact: Olaf Conrad //
42 // Institute of Geography //
43 // University of Goettingen //
44 // Goldschmidtstr. 5 //
45 // 37077 Goettingen //
46 // Germany //
47 // //
48 ///////////////////////////////////////////////////////////
49
50 //---------------------------------------------------------
51 #include "Morphometry.h"
52
53
54 ///////////////////////////////////////////////////////////
55 // //
56 // //
57 // //
58 ///////////////////////////////////////////////////////////
59
60 //---------------------------------------------------------
CMorphometry(void)61 CMorphometry::CMorphometry(void)
62 {
63 Set_Name (_TL("Slope, Aspect, Curvature"));
64
65 Set_Author ("O.Conrad (c) 2001");
66
67 Set_Description (_TW(
68 "Calculates the local morphometric terrain parameters slope, aspect and if supported "
69 "by the chosen method also the curvature. Besides tangential curvature also its "
70 "horizontal and vertical components (i.e. plan and profile curvature) can be calculated."
71 ));
72
73 Add_Reference("Travis, M.R., Elsner, G.H., Iverson, W.D., Johnson, C.G.", "1975",
74 "VIEWIT: computation of seen areas, slope, and aspect for land-use planning",
75 "USDA F.S. Gen. Tech. Rep. PSW-11/1975, 70p. Berkeley, California, U.S.A."
76 );
77
78 Add_Reference("Tarboton, D.G.", "1997",
79 "A new method for the determination of flow directions and upslope areas in grid digital elevation models",
80 "Water Resources Research, Vol.33, No.2, p.309-319."
81 );
82
83 Add_Reference("Horn, B. K.", "1981",
84 "Hill shading and the relectance map",
85 "Proceedings of the IEEE, v. 69, no. 1, p.14-47."
86 );
87
88 Add_Reference("Beasley, D.B., Huggins, L.F.", "1982",
89 "ANSWERS: User's manual",
90 "U.S. EPA-905/9-82-001, Chicago, IL. 54pp."
91 );
92
93 Add_Reference("Costa-Cabral, M., Burges, S.J.", "1994",
94 "Digital Elevation Model Networks (DEMON): a model of flow over hillslopes for computation of contributing and dispersal areas",
95 "Water Resources Research, v. 30, no. 6, p.1681-1692."
96 );
97
98 Add_Reference("Evans, I.S.", "1979",
99 "An integrated system of terrain analysis and slope mapping",
100 "Final report on grant DA-ERO-591-73-G0040, University of Durham, England."
101 );
102
103 Add_Reference("Bauer, J., Rohdenburg, H., Bork, H.-R.", "1985",
104 "Ein Digitales Reliefmodell als Vorraussetzung fuer ein deterministisches Modell der Wasser- und Stoff-Fluesse",
105 "In: Bork, H.-R., Rohdenburg, H. [Eds.]: Parameteraufbereitung fuer deterministische Gebietswassermodelle, Grundlagenarbeiten zur Analyse von Agrar-Oekosystemen, Landschaftsgenese und Landschaftsoekologie, H.10, p.1-15."
106 );
107
108 Add_Reference("Heerdegen, R.G., Beran, M.A.", "1982",
109 "Quantifying source areas through land surface curvature",
110 "Journal of Hydrology, Vol.57."
111 );
112
113 Add_Reference("Olaya, V.", "2006",
114 "Basic Land-Surface Parameters",
115 "In: Hengl, T., Reuter, H.I. [Eds.]: Geomorphometry: Concepts, Software, Applications. Developments in Soil Science, Elsevier, Vol.33, 141-169."
116 );
117
118 Add_Reference("Zevenbergen, L.W., Thorne, C.R.", "1987",
119 "Quantitative analysis of land surface topography",
120 "Earth Surface Processes and Landforms, 12: 47-56."
121 );
122
123 Add_Reference("Haralick, R.M.", "1983",
124 "Ridge and valley detection on digital images",
125 "Computer Vision, Graphics and Image Processing, Vol.22, No.1, p.28-38."
126 );
127
128 Add_Reference("Florinsky, I.V.", "2009",
129 "Computation of the third?order partial derivatives from a digital elevation model",
130 "International Journal of Geographical Information Science, 23(2), p.213-231.",
131 SG_T("https://doi.org/10.1080/13658810802527499"), SG_T("doi: 10.1080/13658810802527499")
132 );
133
134 //-----------------------------------------------------
135 Parameters.Add_Grid(
136 "", "ELEVATION" , _TL("Elevation"),
137 _TL(""),
138 PARAMETER_INPUT
139 );
140
141 //-----------------------------------------------------
142 Parameters.Add_Grid(
143 "", "SLOPE" , _TL("Slope"),
144 _TL(""),
145 PARAMETER_OUTPUT
146 );
147
148 Parameters.Add_Grid(
149 "", "ASPECT" , _TL("Aspect"),
150 _TL(""),
151 PARAMETER_OUTPUT
152 );
153
154 Parameters.Add_Grid(
155 "", "C_GENE" , _TL("General Curvature"),
156 _TL(""),
157 PARAMETER_OUTPUT_OPTIONAL
158 );
159
160 Parameters.Add_Grid(
161 "", "C_PROF" , _TL("Profile Curvature"),
162 _TL(""),
163 PARAMETER_OUTPUT_OPTIONAL
164 );
165
166 Parameters.Add_Grid(
167 "", "C_PLAN" , _TL("Plan Curvature"),
168 _TL(""),
169 PARAMETER_OUTPUT_OPTIONAL
170 );
171
172 Parameters.Add_Grid(
173 "", "C_TANG" , _TL("Tangential Curvature"),
174 _TL(""),
175 PARAMETER_OUTPUT_OPTIONAL
176 );
177
178 Parameters.Add_Grid(
179 "", "C_LONG" , _TL("Longitudinal Curvature"),
180 _TL("Zevenbergen & Thorne (1987) refer to this as profile curvature"),
181 PARAMETER_OUTPUT_OPTIONAL
182 );
183
184 Parameters.Add_Grid(
185 "", "C_CROS" , _TL("Cross-Sectional Curvature"),
186 _TL("Zevenbergen & Thorne (1987) refer to this as plan curvature"),
187 PARAMETER_OUTPUT_OPTIONAL
188 );
189
190 Parameters.Add_Grid(
191 "", "C_MINI" , _TL("Minimal Curvature"),
192 _TL(""),
193 PARAMETER_OUTPUT_OPTIONAL
194 );
195
196 Parameters.Add_Grid(
197 "", "C_MAXI" , _TL("Maximal Curvature"),
198 _TL(""),
199 PARAMETER_OUTPUT_OPTIONAL
200 );
201
202 Parameters.Add_Grid(
203 "", "C_TOTA" , _TL("Total Curvature"),
204 _TL(""),
205 PARAMETER_OUTPUT_OPTIONAL
206 );
207
208 Parameters.Add_Grid(
209 "", "C_ROTO" , _TL("Flow Line Curvature"),
210 _TL(""),
211 PARAMETER_OUTPUT_OPTIONAL
212 );
213
214 //-----------------------------------------------------
215 Parameters.Add_Choice(
216 "", "METHOD" , _TL("Method"),
217 _TL(""),
218 CSG_String::Format("%s|%s|%s|%s|%s|%s|%s|%s|%s",
219 _TL("maximum slope (Travis et al. 1975)"),
220 _TL("maximum triangle slope (Tarboton 1997)"),
221 _TL("least squares fitted plane (Horn 1981, Costa-Cabral & Burgess 1996)"),
222 _TL("6 parameter 2nd order polynom (Evans 1979)"),
223 _TL("6 parameter 2nd order polynom (Heerdegen & Beran 1982)"),
224 _TL("6 parameter 2nd order polynom (Bauer, Rohdenburg, Bork 1985)"),
225 _TL("9 parameter 2nd order polynom (Zevenbergen & Thorne 1987)"),
226 _TL("10 parameter 3rd order polynom (Haralick 1983)"),
227 _TL("10 parameter 3rd order polynom (Florinsky 2009)")
228 ), 6
229 );
230
231 Parameters.Add_Choice(
232 "SLOPE" , "UNIT_SLOPE" , _TL("Unit"),
233 _TL(""),
234 CSG_String::Format("%s|%s|%s",
235 _TL("radians"),
236 _TL("degree"),
237 _TL("percent rise")
238 ), 0
239 );
240
241 Parameters.Add_Choice(
242 "ASPECT", "UNIT_ASPECT" , _TL("Unit"),
243 _TL(""),
244 CSG_String::Format("%s|%s",
245 _TL("radians"),
246 _TL("degree")
247 ), 0
248 );
249 }
250
251
252 ///////////////////////////////////////////////////////////
253 // //
254 ///////////////////////////////////////////////////////////
255
256 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)257 int CMorphometry::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
258 {
259 if( pParameter->Cmp_Identifier("METHOD") )
260 {
261 bool bOn;
262
263 bOn = pParameter->asInt() >= 3 || pParameter->asInt() == 0;
264 pParameters->Set_Enabled("C_GENE", bOn);
265 pParameters->Set_Enabled("C_PROF", bOn);
266 pParameters->Set_Enabled("C_PLAN", bOn);
267
268 bOn = pParameter->asInt() >= 3;
269 pParameters->Set_Enabled("C_TANG", bOn);
270 pParameters->Set_Enabled("C_LONG", bOn);
271 pParameters->Set_Enabled("C_CROS", bOn);
272 pParameters->Set_Enabled("C_MINI", bOn);
273 pParameters->Set_Enabled("C_MAXI", bOn);
274 pParameters->Set_Enabled("C_TOTA", bOn);
275 }
276
277 return( CSG_Tool::On_Parameters_Enable(pParameters, pParameter) );
278 }
279
280
281 ///////////////////////////////////////////////////////////
282 // //
283 ///////////////////////////////////////////////////////////
284
285 //---------------------------------------------------------
On_Execute(void)286 bool CMorphometry::On_Execute(void)
287 {
288 m_pDTM = Parameters("ELEVATION")->asGrid();
289
290 m_pSlope = Parameters("SLOPE" )->asGrid();
291 m_pAspect = Parameters("ASPECT" )->asGrid();
292
293 m_pC_Gene = Parameters("C_GENE" )->asGrid();
294 m_pC_Prof = Parameters("C_PROF" )->asGrid();
295 m_pC_Plan = Parameters("C_PLAN" )->asGrid();
296 m_pC_Tang = Parameters("C_TANG" )->asGrid();
297 m_pC_Long = Parameters("C_LONG" )->asGrid();
298 m_pC_Cros = Parameters("C_CROS" )->asGrid();
299 m_pC_Mini = Parameters("C_MINI" )->asGrid();
300 m_pC_Maxi = Parameters("C_MAXI" )->asGrid();
301 m_pC_Tota = Parameters("C_TOTA" )->asGrid();
302 m_pC_Roto = Parameters("C_ROTO" )->asGrid();
303
304 int Method = Parameters("METHOD" )->asInt ();
305
306 if( Method == 0 )
307 {
308 m_pC_Tang = m_pC_Long = m_pC_Cros = m_pC_Mini = m_pC_Maxi = m_pC_Tota = m_pC_Roto = NULL;
309 }
310 else if( Method < 3 )
311 {
312 m_pC_Gene = m_pC_Prof = m_pC_Plan =
313 m_pC_Tang = m_pC_Long = m_pC_Cros = m_pC_Mini = m_pC_Maxi = m_pC_Tota = m_pC_Roto = NULL;
314 }
315
316 //-----------------------------------------------------
317 DataObject_Set_Colors(m_pSlope , 11, SG_COLORS_RED_GREEN , true);
318 DataObject_Set_Colors(m_pAspect, 11, SG_COLORS_ASPECT_3 , false);
319 DataObject_Set_Colors(m_pC_Gene, 11, SG_COLORS_RED_GREY_BLUE, true);
320 DataObject_Set_Colors(m_pC_Prof, 11, SG_COLORS_RED_GREY_BLUE, true);
321 DataObject_Set_Colors(m_pC_Plan, 11, SG_COLORS_RED_GREY_BLUE, true);
322 DataObject_Set_Colors(m_pC_Tang, 11, SG_COLORS_RED_GREY_BLUE, true);
323 DataObject_Set_Colors(m_pC_Long, 11, SG_COLORS_RED_GREY_BLUE, true);
324 DataObject_Set_Colors(m_pC_Cros, 11, SG_COLORS_RED_GREY_BLUE, true);
325 DataObject_Set_Colors(m_pC_Mini, 11, SG_COLORS_RED_GREY_BLUE, true);
326 DataObject_Set_Colors(m_pC_Maxi, 11, SG_COLORS_RED_GREY_BLUE, true);
327 DataObject_Set_Colors(m_pC_Tota, 11, SG_COLORS_YELLOW_RED , false);
328 DataObject_Set_Colors(m_pC_Roto, 11, SG_COLORS_RED_GREY_BLUE, true);
329
330 //-----------------------------------------------------
331 m_Unit_Slope = Parameters("UNIT_SLOPE" )->asInt();
332
333 if( m_Unit_Slope == 0 )
334 {
335 m_pSlope->Set_Unit(_TL("Radians"));
336 }
337 else if( m_Unit_Slope == 1 )
338 {
339 m_pSlope->Set_Unit(_TL("Degree"));
340 }
341 else // if( m_Unit_Slope == 2 )
342 {
343 m_pSlope->Set_Unit(_TL("Percent"));
344 }
345
346 //-----------------------------------------------------
347 m_Unit_Aspect = Parameters("UNIT_ASPECT")->asInt();
348
349 if( m_Unit_Aspect == 0 )
350 {
351 m_pAspect->Set_Unit(_TL("Radians"));
352 }
353 else // if( m_Unit_Aspect == 1 )
354 {
355 m_pAspect->Set_Unit(_TL("Degree"));
356 }
357
358 //-----------------------------------------------------
359 for(int y=0; y<Get_NY() && Set_Progress(y); y++)
360 {
361 #pragma omp parallel for
362 for(int x=0; x<Get_NX(); x++)
363 {
364 if( m_pDTM->is_NoData(x, y) )
365 {
366 Set_NoData(x, y);
367 }
368 else switch( Method )
369 {
370 case 0: Set_MaximumSlope(x, y); break;
371 case 1: Set_Tarboton (x, y); break;
372 case 2: Set_LeastSquare (x, y); break;
373 case 3: Set_Evans (x, y); break;
374 case 4: Set_Heerdegen (x, y); break;
375 case 5: Set_BRM (x, y); break;
376 default: Set_Zevenbergen (x, y); break;
377 case 7: Set_Haralick (x, y); break;
378 case 8: Set_Florinsky (x, y); break;
379 }
380 }
381 }
382
383 return( true );
384 }
385
386
387 ///////////////////////////////////////////////////////////
388 // //
389 // //
390 // //
391 ///////////////////////////////////////////////////////////
392
393 //---------------------------------------------------------
394 // Indexing of the Submatrix:
395 //
396 // +-------+ +-------+ +-------+
397 // | 7 0 1 | | 2 5 8 | | 8 5 2 |
398 // | 6 * 2 | => | 1 4 7 | or | 7 4 1 |
399 // | 5 4 3 | | 0 3 6 | | 6 3 0 |
400 // +-------+ +-------+ +-------+
401 //
402 //---------------------------------------------------------
Get_SubMatrix3x3(int x,int y,double Z[9],int Orientation)403 inline void CMorphometry::Get_SubMatrix3x3(int x, int y, double Z[9], int Orientation)
404 {
405 static const int Indexes[][8] =
406 {
407 { 5, 8, 7, 6, 3, 0, 1, 2 },
408 { 5, 2, 1, 0, 3, 6, 7, 8 }
409 };
410
411 int *Index = (int *)Indexes[Orientation];
412
413 double z = m_pDTM->asDouble(x, y);
414
415 Z[4] = 0.;
416
417 for(int i=0; i<8; i++)
418 {
419 int ix = Get_xTo(i, x);
420 int iy = Get_yTo(i, y);
421
422 if( m_pDTM->is_InGrid(ix, iy) )
423 {
424 Z[Index[i]] = m_pDTM->asDouble(ix, iy) - z;
425 }
426 else
427 {
428 ix = Get_xTo(i + 4, x);
429 iy = Get_yTo(i + 4, y);
430
431 if( m_pDTM->is_InGrid(ix, iy) )
432 {
433 Z[Index[i]] = z - m_pDTM->asDouble(ix, iy);
434 }
435 else
436 {
437 Z[Index[i]] = 0.;
438 }
439 }
440 }
441 }
442
443 //---------------------------------------------------------
Get_SubMatrix5x5(int x,int y,double Z[25],int Orientation)444 inline void CMorphometry::Get_SubMatrix5x5(int x, int y, double Z[25], int Orientation)
445 {
446 double z = m_pDTM->asDouble(x,y);
447
448 if( Orientation == 0 )
449 {
450 for(int i=0, iy=y-2; iy<=y+2; iy++)
451 {
452 int jy = iy < 0 ? 0 : (iy >= Get_NY() ? Get_NY() - 1 : iy);
453
454 for(int ix=x-2; ix<=x+2; ix++, i++)
455 {
456 int jx = ix < 0 ? 0 : (ix >= Get_NX() ? Get_NX() - 1 : ix);
457
458 Z[i] = m_pDTM->is_InGrid(jx, jy) ? m_pDTM->asDouble(jx, jy) - z : 0.;
459 }
460 }
461 }
462 else
463 {
464 for(int i=0, iy=y+2; iy>=y-2; iy--)
465 {
466 int jy = iy < 0 ? 0 : (iy >= Get_NY() ? Get_NY() - 1 : iy);
467
468 for(int ix=x-2; ix<=x+2; ix++, i++)
469 {
470 int jx = ix < 0 ? 0 : (ix >= Get_NX() ? Get_NX() - 1 : ix);
471
472 Z[i] = m_pDTM->is_InGrid(jx, jy) ? m_pDTM->asDouble(jx, jy) - z : 0.;
473 }
474 }
475 }
476 }
477
478
479 ///////////////////////////////////////////////////////////
480 // //
481 ///////////////////////////////////////////////////////////
482
483 //---------------------------------------------------------
484 #define SET_NODATA(grid) if( grid ) grid->Set_NoData(x, y);
485 #define SET_VALUE(grid, value) if( grid ) grid->Set_Value(x, y, value);
486
487 //---------------------------------------------------------
Set_NoData(int x,int y)488 inline void CMorphometry::Set_NoData(int x, int y)
489 {
490 SET_NODATA(m_pSlope )
491 SET_NODATA(m_pAspect)
492 SET_NODATA(m_pC_Gene)
493 SET_NODATA(m_pC_Prof)
494 SET_NODATA(m_pC_Plan)
495 SET_NODATA(m_pC_Tang)
496 SET_NODATA(m_pC_Long)
497 SET_NODATA(m_pC_Cros)
498 SET_NODATA(m_pC_Mini)
499 SET_NODATA(m_pC_Maxi)
500 SET_NODATA(m_pC_Tota)
501 SET_NODATA(m_pC_Roto)
502 }
503
504 //---------------------------------------------------------
Set_Gradient(int x,int y,double Slope,double Aspect)505 inline void CMorphometry::Set_Gradient(int x, int y, double Slope, double Aspect)
506 {
507 if ( m_Unit_Slope == 2 )
508 {
509 SET_VALUE(m_pSlope, 100. * Slope);
510 }
511 else if( m_Unit_Slope == 1 )
512 {
513 SET_VALUE(m_pSlope, atan(Slope) * M_RAD_TO_DEG);
514 }
515 else
516 {
517 SET_VALUE(m_pSlope, atan(Slope));
518 }
519
520 //-----------------------------------------------------
521 if( m_Unit_Aspect == 1 && Aspect >= 0. )
522 {
523 SET_VALUE(m_pAspect, Aspect * M_RAD_TO_DEG);
524 }
525 else
526 {
527 SET_VALUE(m_pAspect, Aspect);
528 }
529 }
530
531
532 ///////////////////////////////////////////////////////////
533 // //
534 ///////////////////////////////////////////////////////////
535
536 //---------------------------------------------------------
Set_From_Polynom(int x,int y,double r,double t,double s,double p,double q)537 inline void CMorphometry::Set_From_Polynom(int x, int y, double r, double t, double s, double p, double q)
538 {
539 //-----------------------------------------------------
540 double p2_q2 = p*p + q*q;
541
542 Set_Gradient(x, y, sqrt(p2_q2),
543 p != 0. ? M_PI_180 + atan2(q, p)
544 : q > 0. ? M_PI_270
545 : q < 0. ? M_PI_090
546 : m_pAspect ? m_pAspect->Get_NoData_Value() : -1
547 );
548
549 //-----------------------------------------------------
550 if( p2_q2 )
551 {
552 double spq = s * p * q, p2 = p*p, q2 = q*q;
553
554 SET_VALUE(m_pC_Gene, -2 * (r + t));
555 SET_VALUE(m_pC_Prof, -(r * p2 + t * q2 + 2 * spq) / (p2_q2 * pow(1 + p2_q2, 1.5)));
556 SET_VALUE(m_pC_Plan, -(t * p2 + r * q2 - 2 * spq) / ( pow( p2_q2, 1.5)));
557 SET_VALUE(m_pC_Tang, -(t * p2 + r * q2 - 2 * spq) / (p2_q2 * pow(1 + p2_q2, 0.5)));
558 SET_VALUE(m_pC_Long, -2 * (r * p2 + t * q2 + spq) / (p2_q2 ));
559 SET_VALUE(m_pC_Cros, -2 * (t * p2 + r * q2 - spq) / (p2_q2 ));
560 SET_VALUE(m_pC_Mini, -r/2 - t/2 - sqrt(0.5 * (r - t)*(r - t) + s*s));
561 SET_VALUE(m_pC_Maxi, -r/2 - t/2 + sqrt(0.5 * (r - t)*(r - t) + s*s));
562 SET_VALUE(m_pC_Tota, r*r + 2 * s*s + t*t);
563 SET_VALUE(m_pC_Roto, (p2 - q2) * s - p * q * (r - t)); // rotor
564 // SET_VALUE(m_pC_Gaus, (r * t - 2 * s*s) / (1 + p2_q2)); // total gaussian
565 }
566 }
567
568
569 ///////////////////////////////////////////////////////////
570 // //
571 // The Methods //
572 // //
573 ///////////////////////////////////////////////////////////
574
575 //---------------------------------------------------------
576 // Maximum Slope (Travis et al., 1975, Peucker & Douglas, 1975))
577 //
578 // Travis, M.R., Elsner, G.H., Iverson, W.D., and Johnson, C.G. 1975:
579 // VIEWIT: computation of seen areas, slope, and aspect for land-use planning.
580 // USDA F.S. Gen. Tech. Rep. PSW-11/1975, 70p. Berkeley, California, U.S.A.
581 //
582 //---------------------------------------------------------
Set_MaximumSlope(int x,int y)583 void CMorphometry::Set_MaximumSlope(int x, int y)
584 {
585 int i, ix, iy, j, Aspect;
586 double z, Z[8], Slope, Curv, hCurv, a, b;
587
588 //-----------------------------------------------------
589 z = m_pDTM->asDouble(x, y);
590 Slope = Curv = 0.;
591 Aspect = -1;
592
593 for(i=0; i<8; i++)
594 {
595 if( !m_pDTM->is_InGrid(ix = Get_xTo(i, x), iy = Get_yTo(i, y)) )
596 {
597 Z[i] = 0.;
598 }
599 else
600 {
601 Z[i] = (z - m_pDTM->asDouble(ix, iy)) / Get_Length(i);
602 Curv += Z[i];
603
604 if( Z[i] > Slope )
605 {
606 Aspect = i;
607 Slope = Z[i];
608 }
609 }
610 }
611
612 Set_Gradient(x, y, Slope, Aspect * M_PI_045);
613
614 //-------------------------------------------------
615 if( Aspect < 0. )
616 {
617 SET_NODATA(m_pAspect);
618
619 SET_NODATA(m_pC_Gene);
620 SET_NODATA(m_pC_Prof);
621 SET_NODATA(m_pC_Plan);
622 }
623 else
624 {
625 //---------------------------------------------
626 // Let's now estimate the plan curvature...
627
628 for(i=Aspect+1, j=0, a=0.; i<Aspect+8; i++, j++)
629 {
630 if( Z[i % 8] < 0. )
631 {
632 a = j + Z[(i - 1) % 8] / (Z[(i - 1) % 8] - Z[i % 8]);
633 break;
634 }
635 }
636
637 if( a != 0. )
638 {
639 for(i=Aspect+7, j=0, b=0.; i>Aspect; i--, j++)
640 {
641 if( Z[i % 8] < 0. )
642 {
643 b = j + Z[(i + 1) % 8] / (Z[(i + 1) % 8] - Z[i % 8]);
644 break;
645 }
646 }
647
648 hCurv = 45. * (a + b) - 180.;
649 }
650 else
651 {
652 hCurv = 180.;
653 }
654
655 //---------------------------------------------
656 SET_VALUE(m_pC_Gene, Curv);
657 SET_VALUE(m_pC_Prof, Z[Aspect] + Z[(Aspect + 4) % 8]);
658 SET_VALUE(m_pC_Plan, hCurv);
659 }
660 }
661
662 //---------------------------------------------------------
663 // Maximum Triangle Slope
664 //
665 // Tarboton, D.G. (1997):
666 // 'A new method for the determination of flow directions and upslope areas in grid digital elevation models',
667 // Water Resources Research, Vol.33, No.2, p.309-319
668 //
669 //---------------------------------------------------------
Set_Tarboton(int x,int y)670 void CMorphometry::Set_Tarboton(int x, int y)
671 {
672 int i, ix, iy, j;
673 double z, Z[8], iSlope, iAspect, Slope, Aspect, G, H;
674
675 //-----------------------------------------------------
676 z = m_pDTM->asDouble(x, y);
677
678 for(i=0; i<8; i++)
679 {
680 ix = Get_xTo(i, x);
681 iy = Get_yTo(i, y);
682
683 if( m_pDTM->is_InGrid(ix, iy) )
684 {
685 Z[i] = m_pDTM->asDouble(ix, iy);
686 }
687 else
688 {
689 ix = Get_xTo(i + 4, x);
690 iy = Get_yTo(i + 4, y);
691
692 if( m_pDTM->is_InGrid(ix, iy) )
693 {
694 Z[i] = z - (m_pDTM->asDouble(ix, iy) - z);
695 }
696 else
697 {
698 Z[i] = z;
699 }
700 }
701 }
702
703 //---------------------------------------------
704 Slope = 0.;
705 Aspect = -1.;
706
707 for(i=0, j=1; i<8; i++, j=(j+1)%8)
708 {
709 if( i % 2 ) // i => diagonal
710 {
711 G = (z - Z[j]) / Get_Cellsize();
712 H = (Z[j] - Z[i]) / Get_Cellsize();
713 }
714 else // i => orthogonal
715 {
716 G = (z - Z[i]) / Get_Cellsize();
717 H = (Z[i] - Z[j]) / Get_Cellsize();
718 }
719
720 if( H < 0. )
721 {
722 iAspect = 0.;
723 iSlope = G;
724 }
725 else if( H > G )
726 {
727 iAspect = M_PI_045;
728 iSlope = (z - Z[i % 2 ? i : j]) / (sqrt(2.) * Get_Cellsize());
729 }
730 else
731 {
732 iAspect = atan(H / G);
733 iSlope = sqrt(G*G + H*H);
734 }
735
736 if( iSlope > Slope )
737 {
738 Aspect = i * M_PI_045 + (i % 2 ? M_PI_045 - iAspect : iAspect);
739 Slope = iSlope;
740 }
741 }
742
743 //---------------------------------------------
744 if( Aspect < 0. )
745 {
746 Set_NoData(x, y);
747 }
748 else
749 {
750 Set_Gradient(x, y, Slope, Aspect);
751 }
752 }
753
754 //---------------------------------------------------------
755 // Least Squares or Best Fit Plane (Horn 1981, Beasley & Huggins 1982, Costa-Cabral & Burgess 1994)
756 //
757 // Horn, B. K. (1981):
758 // Hill shading and the relectance map.
759 // Proceedings of the IEEE, v. 69, no. 1, p 14-47.
760 //
761 // Beasley, D.B. and Huggins, L.F. 1982:
762 // ANSWERS: User�s manual.
763 // U.S. EPA-905/9-82-001, Chicago, IL. 54pp.
764 //
765 // Costa-Cabral, M., and Burges, S.J., 1994:
766 // Digital Elevation Model Networks (DEMON): a model of flow over hillslopes for computation of contributing and dispersal areas
767 // Water Resources Research, v. 30, no. 6, p. 1681-1692.
768 //
769 //---------------------------------------------------------
Set_LeastSquare(int x,int y)770 void CMorphometry::Set_LeastSquare(int x, int y)
771 {
772 double Z[9], a, b;
773
774 Get_SubMatrix3x3(x, y, Z);
775
776 a = ((Z[2] + 2 * Z[5] + Z[8]) - (Z[0] + 2 * Z[3] + Z[6])) / (8 * Get_Cellsize());
777 b = ((Z[6] + 2 * Z[7] + Z[8]) - (Z[0] + 2 * Z[1] + Z[2])) / (8 * Get_Cellsize());
778
779 Set_Gradient(x, y, sqrt(a*a + b*b),
780 a != 0. ? M_PI_180 + atan2(b, a)
781 : b > 0. ? M_PI_270
782 : b < 0. ? M_PI_090
783 : m_pAspect ? m_pAspect->Get_NoData_Value() : -1
784 );
785 }
786
787 //---------------------------------------------------------
788 // Quadratic Function Approximation (Heerdegen & Beran, 1984)
789 //
790 // Evans, I.S. (1979):
791 // An integrated system of terrain analysis and slope mapping.
792 // Final report on grant DA-ERO-591-73-G0040. University of Durham, England.
793 //
794 //---------------------------------------------------------
795 // f(z) = Ax^2 + By^2 + Cxy + Dx + Ey + F
796 //
797 //---------------------------------------------------------
Set_Evans(int x,int y)798 void CMorphometry::Set_Evans(int x, int y)
799 {
800 double Z[9], A, B, C, D, E;
801
802 Get_SubMatrix3x3(x, y, Z, 1);
803
804 A = (Z[0] + Z[2] + Z[3] + Z[5] + Z[6] + Z[8] - 2 * (Z[1] + Z[4] + Z[7])) / (6 * Get_Cellarea());
805 B = (Z[0] + Z[1] + Z[2] + Z[6] + Z[7] + Z[8] - 2 * (Z[3] + Z[4] + Z[5])) / (6 * Get_Cellarea());
806 C = (Z[2] + Z[6] - Z[0] - Z[8]) / (4 * Get_Cellarea());
807 D = (Z[2] + Z[5] + Z[8] - Z[0] - Z[3] - Z[6]) / (6 * Get_Cellsize());
808 E = (Z[0] + Z[1] + Z[2] - Z[6] - Z[7] - Z[8]) / (6 * Get_Cellsize());
809
810 Set_From_Polynom(x, y, A, B, C, D, E);
811 }
812
813 //---------------------------------------------------------
814 // Quadratic Function Approximation (Heerdegen & Beran, 1984)
815 //
816 // Heerdegen, R.G. / Beran, M.A. (1982):
817 // Quantifying source areas through land surface curvature.
818 // Journal of Hydrology, Vol.57
819 //
820 //---------------------------------------------------------
821 // f(z) = Ax^2 + By^2 + Cxy + Dx + Ey + F
822 //
823 //---------------------------------------------------------
Set_Heerdegen(int x,int y)824 void CMorphometry::Set_Heerdegen(int x, int y)
825 {
826 double Z[9], A, B, C, D, E, a, b;
827
828 Get_SubMatrix3x3(x, y, Z);
829
830 a = Z[0] + Z[2] + Z[3] + Z[5] + Z[6] + Z[8];
831 b = Z[0] + Z[1] + Z[2] + Z[6] + Z[7] + Z[8];
832
833 A = (0.3 * a - 0.2 * b) / ( Get_Cellarea());
834 B = (0.3 * b - 0.2 * a) / ( Get_Cellarea());
835 C = ( Z[0] - Z[2] - Z[6] + Z[8]) / (4 * Get_Cellarea());
836 D = (-Z[0] + Z[2] - Z[3] + Z[5] - Z[6] + Z[8]) / (6 * Get_Cellsize());
837 E = (-Z[0] - Z[1] - Z[2] + Z[6] + Z[7] + Z[8]) / (6 * Get_Cellsize());
838
839 Set_From_Polynom(x, y, A, B, C, D, E);
840 }
841
842 //---------------------------------------------------------
843 // Quadratic Function Approximation (Bauer, Rohdenburg & Bork, 1985)
844 //
845 // Bauer, J. / Rohdenburg, H. / Bork, H.-R., (1985):
846 // 'Ein Digitales Reliefmodell als Vorraussetzung fuer ein deterministisches Modell der Wasser- und Stoff-Fluesse',
847 // Landschaftsgenese und Landschaftsoekologie, H.10, Parameteraufbereitung fuer deterministische Gebiets-Wassermodelle,
848 // Grundlagenarbeiten zu Analyse von Agrar-Oekosystemen, (Eds.: Bork, H.-R. / Rohdenburg, H.), p.1-15
849 //
850 //---------------------------------------------------------
851 // f(z) = Ax^2 + By^2 + Cxy + Dx + Ey + F
852 //
853 //---------------------------------------------------------
Set_BRM(int x,int y)854 void CMorphometry::Set_BRM(int x, int y)
855 {
856 double Z[9], A, B, C, D, E;
857
858 Get_SubMatrix3x3(x, y, Z);
859
860 A = ( (Z[0] + Z[2] + Z[3] + Z[5] + Z[6] + Z[8]) - 2 * (Z[1] + Z[4] + Z[7]) ) / ( Get_Cellarea());
861 B = ( (Z[0] + Z[6] + Z[1] + Z[7] + Z[2] + Z[8]) - 2 * (Z[3] + Z[4] + Z[5]) ) / ( Get_Cellarea());
862 C = ( Z[8] + Z[0] - Z[7] ) / (4 * Get_Cellarea());
863 D = ( (Z[2] - Z[0]) + (Z[5] - Z[3]) + (Z[8]-Z[6]) ) / (6 * Get_Cellsize());
864 E = ( (Z[6] - Z[0]) + (Z[7] - Z[1]) + (Z[8]-Z[2]) ) / (6 * Get_Cellsize());
865
866 Set_From_Polynom(x, y, A, B, C, D, E);
867 }
868
869 //---------------------------------------------------------
870 // Quadratic Function Approximation (Zevenbergen und Thorne, 1986)
871 //
872 // Zevenbergen, L.W. and C.R. Thorne. 1987:
873 // Quantitative analysis of land surface topography
874 // Earth Surface Processes and Landforms, 12: 47-56.
875 //
876 //---------------------------------------------------------
877 // f(z) = Ax^2y^2 + Bx^2y + Cxy^2 + Dx^2 + Ey^2 + Fxy + Gx + Hy + I
878 //
879 //---------------------------------------------------------
Set_Zevenbergen(int x,int y)880 void CMorphometry::Set_Zevenbergen(int x, int y)
881 {
882 double Z[9], D, E, F, G, H;
883
884 Get_SubMatrix3x3(x, y, Z);
885
886 D = ((Z[3] + Z[5]) / 2. - Z[4]) / ( Get_Cellarea());
887 E = ((Z[1] + Z[7]) / 2. - Z[4]) / ( Get_Cellarea());
888 F = (Z[0] - Z[2] - Z[6] + Z[8]) / (4 * Get_Cellarea());
889 G = (Z[5] - Z[3]) / (2 * Get_Cellsize());
890 H = (Z[7] - Z[1]) / (2 * Get_Cellsize());
891
892 Set_From_Polynom(x, y, D, E, F, G, H);
893 }
894
895 //---------------------------------------------------------
896 // Cubic Function Approximation (Haralick, 1991)
897 //
898 // R.M. Haralick (1983):
899 // 'Ridge and Valley Detection on digital images',
900 // Computer Vision, Graphics and Image Processing, Vol.22, No.1, p.28-38
901 //
902 //---------------------------------------------------------
903 // f(z) = Ax^3 + By^3 + Cx^2y + Dxy^2 + Ex^2 + Fy^2 + Gxy + Hx + Iy + J
904 //
905 //---------------------------------------------------------
Set_Haralick(int x,int y)906 void CMorphometry::Set_Haralick(int x, int y)
907 {
908 //-----------------------------------------------------
909 // Matrices for Finite Difference solution...
910
911 const int Mtrx[][5][5] = {
912 { { 31,- 5,-17,- 5, 31}, {-44,-62,-68,-62,-44}, { 0, 0, 0, 0, 0}, { 44, 62, 68, 62, 44}, {-31, 5, 17, 5,-31} },
913 { { 31,-44, 0, 44,-31}, {- 5,-62, 0, 62, 5}, {-17,-68, 0, 68, 17}, {- 5,-62, 0, 62, 5}, { 31,-44, 0, 44,-31} },
914 { { 2, 2, 2, 2, 2}, {- 1,- 1,- 1,- 1,- 1}, {- 2,- 2,- 2,- 2,- 2}, {- 1,- 1,- 1,- 1,- 1}, { 2, 2, 2, 2, 2} },
915 { { 4, 2, 0,- 2,- 4}, { 2, 1, 0,- 1,- 2}, { 0, 0, 0, 0, 0}, {- 2,- 1, 0, 1, 2}, {- 4,- 2, 0, 2, 4} },
916 { { 2,- 1,- 2,- 1, 2}, { 2,- 1,- 2,- 1, 2}, { 2,- 1,- 2,- 1, 2}, { 2,- 1,- 2,- 1, 2}, { 2,- 1,- 2,- 1, 2} }, };
917
918 //-----------------------------------------------------
919 double Z[25], k[5];
920
921 Get_SubMatrix5x5(x, y, Z);
922
923 for(int i=0; i<5; i++)
924 {
925 k[i] = 0.;
926
927 for(int iy=0, n=0; iy<5; iy++)
928 {
929 for(int ix=0; ix<5; ix++, n++)
930 {
931 k[i] += Z[n] * Mtrx[i][ix][iy];
932 }
933 }
934 }
935
936 k[0] /= 420. * Get_Cellsize();
937 k[1] /= 420. * Get_Cellsize();
938 k[2] /= 70. * Get_Cellarea();
939 k[3] /= 100. * Get_Cellarea();
940 k[4] /= 70. * Get_Cellarea();
941
942 Set_From_Polynom(x, y, k[4], k[2], k[3], k[1], k[0]);
943 }
944
945 //---------------------------------------------------------
946 // Cubic Function Approximation (Florinsky, 2009)
947 //
948 // I.V. Florinsky (2009):
949 // 'Computation of the third-order partial derivatives from a digital elevation model',
950 // International Journal of Geographical Information Science, 23(2), p.213-231.
951 //
952 //---------------------------------------------------------
953 // f(z) = a/6x^3 + d/6y^3 + b/2x^2y + c/2xy^2 + r/2x^2 + t/2y^2 + sxy + px + qy + u
954 //
955 //---------------------------------------------------------
Set_Florinsky(int x,int y)956 void CMorphometry::Set_Florinsky(int x, int y)
957 {
958 double z[26], r, t, s, p, q;
959
960 Get_SubMatrix5x5(x, y, z + 1, 1);
961
962 r = ( (z[ 1] + z[ 5] + z[ 6] + z[10] + z[11] + z[15] + z[16] + z[20] + z[21] + z[25]) * 2.
963 - (z[ 3] + z[ 8] + z[13] + z[18] + z[23]) * 2.
964 - z[ 2] - z[ 4] - z[ 7] - z[ 9] - z[12] - z[14] - z[17] - z[19] - z[22] - z[24]
965 ) / ( 35. * Get_Cellarea());
966
967 t = ( (z[ 1] + z[ 2] + z[ 3] + z[ 4] + z[ 5] + z[21] + z[22] + z[23] + z[24] + z[25]) * 2.
968 - (z[11] + z[12] + z[13] + z[14] + z[15]) * 2.
969 - z[ 6] - z[ 7] - z[ 8] - z[ 9] - z[10] - z[16] - z[17] - z[18] - z[19] - z[20]
970 ) / ( 35. * Get_Cellarea());
971
972 s = ( z[ 9] + z[17] - z[ 7] - z[19]
973 + (z[ 5] + z[21] - z[ 1] - z[25]) * 4.
974 + (z[ 4] + z[10] + z[16] + z[22] - z[ 2] - z[ 6] - z[20] - z[24]) * 2.
975 ) / (100. * Get_Cellarea());
976
977 p = ( (z[ 4] + z[24] - z[ 2] - z[22]) * 44.
978 + (z[ 1] + z[21] - z[ 5] - z[25] + (z[ 9] + z[19] - z[ 7] - z[17]) * 2.) * 31.
979 + (z[15] - z[11] +(z[14] - z[12]) * 4.) * 17.
980 + (z[10] + z[20] - z[ 6] - z[16]) * 5.
981 ) / (420. * Get_Cellsize());
982
983 q = ( (z[ 6] + z[10] - z[16] - z[20]) * 44.
984 + (z[21] + z[25] - z[ 1] - z[ 5] + (z[ 7] + z[ 9] - z[17] - z[19]) * 2.) * 31.
985 + (z[ 3] - z[23] +(z[ 8] - z[18]) * 4.) * 17.
986 + (z[ 2] + z[ 4] - z[22] - z[24]) * 5.
987 ) / (420. * Get_Cellsize());
988
989 Set_From_Polynom(x, y, r / 2., t / 2., s, q, p);
990 }
991
992
993 ///////////////////////////////////////////////////////////
994 // //
995 // //
996 // //
997 ///////////////////////////////////////////////////////////
998
999 //---------------------------------------------------------
1000