1 
2 ///////////////////////////////////////////////////////////
3 //                                                       //
4 //                         SAGA                          //
5 //                                                       //
6 //      System for Automated Geoscientific Analyses      //
7 //                                                       //
8 //                     Tool Library                      //
9 //                    ta_morphometry                     //
10 //                                                       //
11 //-------------------------------------------------------//
12 //                                                       //
13 //                    Ruggedness.cpp                     //
14 //                                                       //
15 //                 Copyright (C) 2010 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 Hamburg                  //
44 //                Germany                                //
45 //                                                       //
46 ///////////////////////////////////////////////////////////
47 
48 //---------------------------------------------------------
49 #include "ruggedness.h"
50 
51 
52 ///////////////////////////////////////////////////////////
53 //														 //
54 //														 //
55 //														 //
56 ///////////////////////////////////////////////////////////
57 
58 //---------------------------------------------------------
CRuggedness_TRI(void)59 CRuggedness_TRI::CRuggedness_TRI(void)
60 {
61 	Set_Name		(_TL("Terrain Ruggedness Index (TRI)"));
62 
63 	Set_Author		("O.Conrad (c) 2010");
64 
65 	Set_Description	(_TW(
66 		"Terrain Ruggedness Index (TRI)"
67 	));
68 
69 	Add_Reference("Riley, S.J., De Gloria, S.D., Elliot, R.", "1999",
70 		"A Terrain Ruggedness that Quantifies Topographic Heterogeneity",
71 		"Intermountain Journal of Science, Vol.5, No.1-4, pp.23-27.",
72 		SG_T("http://download.osgeo.org/qgis/doc/reference-docs/Terrain_Ruggedness_Index.pdf"), SG_T("online")
73 	);
74 
75 	//-----------------------------------------------------
76 	Parameters.Add_Grid(
77 		"", "DEM"		, _TL("Elevation"),
78 		_TL(""),
79 		PARAMETER_INPUT
80 	);
81 
82 	Parameters.Add_Grid(
83 		"", "TRI"		, _TL("Terrain Ruggedness Index (TRI)"),
84 		_TL(""),
85 		PARAMETER_OUTPUT
86 	);
87 
88 	Parameters.Add_Choice(
89 		"", "MODE"		, _TL("Search Mode"),
90 		_TL(""),
91 		CSG_String::Format("%s|%s",
92 			_TL("Square"),
93 			_TL("Circle")
94 		), 1
95 	);
96 
97 	Parameters.Add_Int(
98 		"", "RADIUS"	, _TL("Search Radius"),
99 		_TL("radius in cells"),
100 		1, 1, true
101 	);
102 
103 	m_Cells.Get_Weighting().Set_BandWidth(75.);	// 75%
104 	m_Cells.Get_Weighting().Create_Parameters(Parameters);
105 }
106 
107 
108 ///////////////////////////////////////////////////////////
109 //														 //
110 ///////////////////////////////////////////////////////////
111 
112 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)113 int CRuggedness_TRI::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
114 {
115 	m_Cells.Get_Weighting().Enable_Parameters(*pParameters);
116 
117 	return( CSG_Tool_Grid::On_Parameters_Enable(pParameters, pParameter) );
118 }
119 
120 
121 ///////////////////////////////////////////////////////////
122 //														 //
123 ///////////////////////////////////////////////////////////
124 
125 //---------------------------------------------------------
On_Execute(void)126 bool CRuggedness_TRI::On_Execute(void)
127 {
128 	//-----------------------------------------------------
129 	m_pDEM	= Parameters("DEM")->asGrid();
130 	m_pTRI	= Parameters("TRI")->asGrid();
131 
132 	DataObject_Set_Colors(m_pTRI, 11, SG_COLORS_RED_GREY_BLUE, true);
133 
134 	m_Cells.Get_Weighting().Set_Parameters(Parameters);
135 	m_Cells.Get_Weighting().Set_BandWidth(Parameters("RADIUS")->asInt() * m_Cells.Get_Weighting().Get_BandWidth() / 100.);
136 
137 	if( !m_Cells.Set_Radius(Parameters("RADIUS")->asInt(), Parameters("MODE")->asInt() == 0) )
138 	{
139 		return( false );
140 	}
141 
142 	//-----------------------------------------------------
143 	for(int y=0; y<Get_NY() && Set_Progress(y); y++)
144 	{
145 		#pragma omp parallel for
146 		for(int x=0; x<Get_NX(); x++)
147 		{
148 			Set_Index(x, y);
149 		}
150 	}
151 
152 	//-----------------------------------------------------
153 	m_Cells.Destroy();
154 
155 	return( true );
156 }
157 
158 
159 ///////////////////////////////////////////////////////////
160 //														 //
161 ///////////////////////////////////////////////////////////
162 
163 //---------------------------------------------------------
Set_Index(int x,int y)164 bool CRuggedness_TRI::Set_Index(int x, int y)
165 {
166 	if( m_pDEM->is_InGrid(x, y) )
167 	{
168 		int		i, ix, iy;
169 		double	z, iz, Distance, Weight, n, s;
170 
171 		for(i=0, n=s=0., z=m_pDEM->asDouble(x, y); i<m_Cells.Get_Count(); i++)
172 		{
173 			if( m_Cells.Get_Values(i, ix = x, iy = y, Distance, Weight, true) && Weight > 0. && m_pDEM->is_InGrid(ix, iy) )
174 			{
175 				iz	 = m_pDEM->asDouble(ix, iy);
176 				s	+= SG_Get_Square((z - iz) * Weight);
177 				n	+= Weight;
178 			}
179 		}
180 
181 		//-------------------------------------------------
182 		if( n > 0. )
183 		{
184 			m_pTRI->Set_Value(x, y, sqrt(s / n));
185 
186 			return( true );
187 		}
188 	}
189 
190 	//-----------------------------------------------------
191 	m_pTRI->Set_NoData(x, y);
192 
193 	return( false );
194 }
195 
196 
197 ///////////////////////////////////////////////////////////
198 //														 //
199 //														 //
200 //														 //
201 ///////////////////////////////////////////////////////////
202 
203 //---------------------------------------------------------
CRuggedness_VRM(void)204 CRuggedness_VRM::CRuggedness_VRM(void)
205 {
206 	Set_Name		(_TL("Vector Ruggedness Measure (VRM)"));
207 
208 	Set_Author		("O.Conrad (c) 2010");
209 
210 	Set_Description	(_TW(
211 		"Vector Ruggedness Measure (VRM)"
212 	));
213 
214 	Add_Reference("Sappington, J.M., Longshore, K.M., Thompson, D.B.", "2007",
215 		"Quantifying Landscape Ruggedness for Animal Habitat Analysis: A Case Study Using Bighorn Sheep in the Mojave Desert",
216 		"Journal of Wildlife Management 71(5):1419�1426.",
217 		SG_T("https://wildlife.onlinelibrary.wiley.com/doi/abs/10.2193/2005-723"), SG_T("online")
218 	);
219 
220 	//-----------------------------------------------------
221 	Parameters.Add_Grid(
222 		"", "DEM"		, _TL("Elevation"),
223 		_TL(""),
224 		PARAMETER_INPUT
225 	);
226 
227 	Parameters.Add_Grid(
228 		"", "VRM"		, _TL("Vector Terrain Ruggedness (VRM)"),
229 		_TL(""),
230 		PARAMETER_OUTPUT
231 	);
232 
233 	Parameters.Add_Choice(
234 		"", "MODE"		, _TL("Search Mode"),
235 		_TL(""),
236 		CSG_String::Format("%s|%s",
237 			_TL("Square"),
238 			_TL("Circle")
239 		), 1
240 	);
241 
242 	Parameters.Add_Int(
243 		"", "RADIUS"	, _TL("Search Radius"),
244 		_TL("radius in cells"),
245 		1, 1, true
246 	);
247 
248 	m_Cells.Get_Weighting().Set_BandWidth(75.);	// 75%
249 	m_Cells.Get_Weighting().Create_Parameters(Parameters);
250 }
251 
252 
253 ///////////////////////////////////////////////////////////
254 //														 //
255 ///////////////////////////////////////////////////////////
256 
257 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)258 int CRuggedness_VRM::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
259 {
260 	m_Cells.Get_Weighting().Enable_Parameters(*pParameters);
261 
262 	return( CSG_Tool_Grid::On_Parameters_Enable(pParameters, pParameter) );
263 }
264 
265 
266 ///////////////////////////////////////////////////////////
267 //														 //
268 ///////////////////////////////////////////////////////////
269 
270 //---------------------------------------------------------
On_Execute(void)271 bool CRuggedness_VRM::On_Execute(void)
272 {
273 	//-----------------------------------------------------
274 	m_pDEM	= Parameters("DEM")->asGrid();
275 	m_pVRM	= Parameters("VRM")->asGrid();
276 
277 	DataObject_Set_Colors(m_pVRM, 11, SG_COLORS_RED_GREY_BLUE, true);
278 
279 	m_Cells.Get_Weighting().Set_Parameters(Parameters);
280 	m_Cells.Get_Weighting().Set_BandWidth(Parameters("RADIUS")->asInt() * m_Cells.Get_Weighting().Get_BandWidth() / 100.);
281 
282 	if( !m_Cells.Set_Radius(Parameters("RADIUS")->asInt(), Parameters("MODE")->asInt() == 0) )
283 	{
284 		return( false );
285 	}
286 
287 	//-----------------------------------------------------
288 	int		x, y;
289 
290 	m_X.Create(Get_System(), SG_DATATYPE_Float);
291 	m_Y.Create(Get_System(), SG_DATATYPE_Float);
292 	m_Z.Create(Get_System(), SG_DATATYPE_Float);
293 
294 	for(y=0; y<Get_NY() && Set_Progress(y); y++)
295 	{
296 		#pragma omp parallel for private(x)
297 		for(x=0; x<Get_NX(); x++)
298 		{
299 			double	slope, aspect;
300 
301 			if( m_pDEM->Get_Gradient(x, y, slope, aspect) )
302 			{
303 				m_X.Set_Value(x, y, aspect < 0. ? 0. : sin(slope) * sin(aspect));
304 				m_Y.Set_Value(x, y, aspect < 0. ? 0. : sin(slope) * cos(aspect));
305 				m_Z.Set_Value(x, y, cos(slope));
306 			}
307 			else
308 			{
309 				m_X.Set_NoData(x, y);
310 			}
311 		}
312 	}
313 
314 	//-----------------------------------------------------
315 	for(y=0; y<Get_NY() && Set_Progress(y); y++)
316 	{
317 		#pragma omp parallel for private(x)
318 		for(x=0; x<Get_NX(); x++)
319 		{
320 			Set_Index(x, y);
321 		}
322 	}
323 
324 	//-----------------------------------------------------
325 	m_Cells.Destroy();
326 
327 	m_X.Destroy();
328 	m_Y.Destroy();
329 	m_Z.Destroy();
330 
331 	return( true );
332 }
333 
334 
335 ///////////////////////////////////////////////////////////
336 //														 //
337 ///////////////////////////////////////////////////////////
338 
339 //---------------------------------------------------------
Set_Index(int x,int y)340 bool CRuggedness_VRM::Set_Index(int x, int y)
341 {
342 	if( m_pDEM->is_InGrid(x, y) )
343 	{
344 		int		i, ix, iy;
345 		double	Distance, Weight, n, sx, sy, sz;
346 
347 		for(i=0, n=sx=sy=sz=0.; i<m_Cells.Get_Count(); i++)
348 		{
349 			if( m_Cells.Get_Values(i, ix = x, iy = y, Distance, Weight, true) && Weight > 0. && m_X.is_InGrid(ix, iy) )
350 			{
351 				sx	+= Weight * m_X.asDouble(ix, iy);
352 				sy	+= Weight * m_Y.asDouble(ix, iy);
353 				sz	+= Weight * m_Z.asDouble(ix, iy);
354 				n	+= Weight;
355 			}
356 		}
357 
358 		//-------------------------------------------------
359 		if( n > 0. )
360 		{
361 			m_pVRM->Set_Value(x, y, 1. - sqrt(sx*sx + sy*sy + sz*sz) / n);
362 
363 			return( true );
364 		}
365 	}
366 
367 	//-----------------------------------------------------
368 	m_pVRM->Set_NoData(x, y);
369 
370 	return( false );
371 }
372 
373 
374 ///////////////////////////////////////////////////////////
375 //														 //
376 //														 //
377 //														 //
378 ///////////////////////////////////////////////////////////
379 
380 //---------------------------------------------------------
381