1 
2 ///////////////////////////////////////////////////////////
3 //                                                       //
4 //                         SAGA                          //
5 //                                                       //
6 //      System for Automated Geoscientific Analyses      //
7 //                                                       //
8 //                     Tool Library:                     //
9 //                      grids_tools                      //
10 //                                                       //
11 //-------------------------------------------------------//
12 //                                                       //
13 //                   gridding3d_idw.cpp                  //
14 //                                                       //
15 //                  Olaf Conrad (C) 2019                 //
16 //                                                       //
17 //-------------------------------------------------------//
18 //                                                       //
19 // This file is part of 'SAGA - System for Automated     //
20 // Geoscientific Analyses'. SAGA is free software; you   //
21 // can redistribute it and/or modify it under the terms  //
22 // of the GNU General Public License as published by the //
23 // Free Software Foundation, either version 2 of the     //
24 // License, or (at your option) any later version.       //
25 //                                                       //
26 // SAGA is distributed in the hope that it will be       //
27 // useful, but WITHOUT ANY WARRANTY; without even the    //
28 // implied warranty of MERCHANTABILITY or FITNESS FOR A  //
29 // PARTICULAR PURPOSE. See the GNU General Public        //
30 // License for more details.                             //
31 //                                                       //
32 // You should have received a copy of the GNU General    //
33 // Public License along with this program; if not, see   //
34 // <http://www.gnu.org/licenses/>.                       //
35 //                                                       //
36 //-------------------------------------------------------//
37 //                                                       //
38 //    e-mail:     oconrad@saga-gis.org                   //
39 //                                                       //
40 //    contact:    Olaf Conrad                            //
41 //                Institute of Geography                 //
42 //                University of Hamburg                  //
43 //                Germany                                //
44 //                                                       //
45 ///////////////////////////////////////////////////////////
46 
47 //---------------------------------------------------------
48 #include "gridding3d_idw.h"
49 
50 
51 ///////////////////////////////////////////////////////////
52 //														 //
53 //														 //
54 //														 //
55 ///////////////////////////////////////////////////////////
56 
57 //---------------------------------------------------------
CGridding3D_IDW(void)58 CGridding3D_IDW::CGridding3D_IDW(void)
59 {
60 	//-----------------------------------------------------
61 	Set_Name		(_TL("Inverse Distance Weighted (3D)"));
62 
63 	Set_Author		("O.Conrad (c) 2019");
64 
65 	Set_Description	(_TW(
66 		"Inverse distance weighted interpolation for 3-dimensional data points. "
67 		"Output will be a grid collection with evenly spaced Z-levels "
68 		"representing the 3rd dimension. "
69 	));
70 
71 	//-----------------------------------------------------
72 	Parameters.Add_Shapes("",
73 		"POINTS"	, _TL("Points"),
74 		_TL(""),
75 		PARAMETER_INPUT, SHAPE_TYPE_Point
76 	);
77 
78 	Parameters.Add_Table_Field("POINTS",
79 		"Z_FIELD"	, _TL("Z"),
80 		_TL("")
81 	);
82 
83 	Parameters.Add_Double("POINTS",
84 		"Z_SCALE"	, _TL("Z Factor"),
85 		_TL(""),
86 		1.
87 	);
88 
89 	Parameters.Add_Table_Field("POINTS",
90 		"V_FIELD"	, _TL("Value"),
91 		_TL("")
92 	);
93 
94 	//-----------------------------------------------------
95 	m_Grid_Target.Create(&Parameters, false, "", "TARGET_");
96 
97 	m_Grid_Target.Add_Grids("GRIDS", _TL("Grid Collection"), false, true);
98 
99 	//-----------------------------------------------------
100 	m_Search_Options.Create(&Parameters, "NODE_SEARCH", 1);
101 
102 	//-----------------------------------------------------
103 	m_Weighting.Set_Weighting (SG_DISTWGHT_IDW);
104 	m_Weighting.Set_IDW_Offset(false);
105 	m_Weighting.Set_IDW_Power (2.);
106 
107 	m_Weighting.Create_Parameters(Parameters);
108 }
109 
110 
111 ///////////////////////////////////////////////////////////
112 //														 //
113 ///////////////////////////////////////////////////////////
114 
115 //---------------------------------------------------------
On_Parameter_Changed(CSG_Parameters * pParameters,CSG_Parameter * pParameter)116 int CGridding3D_IDW::On_Parameter_Changed(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
117 {
118 	if( pParameter->Cmp_Identifier("POINTS") )
119 	{
120 		m_Grid_Target.Set_User_Defined(pParameters, pParameter->asShapes());
121 	}
122 
123 	if( pParameter->Cmp_Identifier("POINTS") || pParameter->Cmp_Identifier("Z_FIELD") )
124 	{
125 		CSG_Shapes	*pPoints = (*pParameters)("POINTS")->asShapes();
126 
127 		if( pPoints )
128 		{
129 			int	zField	= pPoints->Get_Vertex_Type() == SG_VERTEX_TYPE_XY ? (*pParameters)("Z_FIELD")->asInt() : -1;
130 
131 			m_Grid_Target.Set_User_Defined_ZLevels(pParameters,
132 				zField < 0 ? pPoints->Get_ZMin() : pPoints->Get_Minimum(zField),
133 				zField < 0 ? pPoints->Get_ZMax() : pPoints->Get_Maximum(zField), 10
134 			);
135 		}
136 	}
137 
138 	m_Grid_Target.On_Parameter_Changed(pParameters, pParameter);
139 
140 	return( CSG_Tool::On_Parameter_Changed(pParameters, pParameter) );
141 }
142 
143 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)144 int CGridding3D_IDW::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
145 {
146 	if( pParameter->Cmp_Identifier("POINTS") )
147 	{
148 		pParameters->Set_Enabled("Z_FIELD", pParameter->asShapes() && pParameter->asShapes()->Get_Vertex_Type() == SG_VERTEX_TYPE_XY);
149 	}
150 
151 	m_Grid_Target   .On_Parameters_Enable(pParameters, pParameter);
152 
153 	m_Search_Options.On_Parameters_Enable(pParameters, pParameter);
154 
155 	m_Weighting.Enable_Parameters(*pParameters);
156 
157 	return( CSG_Tool::On_Parameters_Enable(pParameters, pParameter) );
158 }
159 
160 
161 ///////////////////////////////////////////////////////////
162 //														 //
163 ///////////////////////////////////////////////////////////
164 
165 //---------------------------------------------------------
On_Execute(void)166 bool CGridding3D_IDW::On_Execute(void)
167 {
168 	CSG_Grids	*pGrids	= m_Grid_Target.Get_Grids("GRIDS");
169 
170 	if( pGrids == NULL )
171 	{
172 		return( false );
173 	}
174 
175 	pGrids->Fmt_Name("%s.%s [%s]",
176 		Parameters("POINTS")->asShapes()->Get_Name(),
177 		Parameters("V_FIELD")->asString(), Get_Name().c_str()
178 	);
179 
180 	//-----------------------------------------------------
181 	if( !Initialize() )
182 	{
183 		Finalize();
184 
185 		return( false );
186 	}
187 
188 	double	zScale	= Parameters("Z_SCALE")->asDouble();
189 
190 	//-----------------------------------------------------
191 	for(int x=0; x<pGrids->Get_NX() && Set_Progress(x, pGrids->Get_NX()); x++)
192 	{
193 		#pragma omp parallel for
194 		for(int y=0; y<pGrids->Get_NY(); y++)
195 		{
196 			double	c[3], v;
197 
198 			c[0]	= pGrids->Get_XMin() + x * pGrids->Get_Cellsize();
199 			c[1]	= pGrids->Get_YMin() + y * pGrids->Get_Cellsize();
200 
201 			for(int z=0; z<pGrids->Get_NZ(); z++)
202 			{
203 				c[2]	= pGrids->Get_Z(z) * zScale;
204 
205 				if( Get_Value(c, v) )
206 				{
207 					pGrids->Set_Value(x, y, z, v);
208 				}
209 				else
210 				{
211 					pGrids->Set_NoData(x, y, z);
212 				}
213 			}
214 		}
215 	}
216 
217 	//-----------------------------------------------------
218 	Finalize();
219 
220 	return( true );
221 }
222 
223 
224 ///////////////////////////////////////////////////////////
225 //														 //
226 ///////////////////////////////////////////////////////////
227 
228 //---------------------------------------------------------
Initialize(void)229 bool CGridding3D_IDW::Initialize(void)
230 {
231 	CSG_Shapes	*pPoints	= Parameters("POINTS")->asShapes();
232 
233 	int	Field	= Parameters("V_FIELD")->asInt();
234 
235 	int	zField	= pPoints->Get_Vertex_Type() == SG_VERTEX_TYPE_XY ? Parameters("Z_FIELD")->asInt() : -1;
236 
237 	double	zScale	= Parameters("Z_SCALE")->asDouble();
238 
239 	if( zScale == 0. )
240 	{
241 		Error_Set(_TL("Z factor is zero! Please use 2D instead of 3D interpolation."));
242 
243 		return( false );
244 	}
245 
246 	//-----------------------------------------------------
247 	m_Points.Create(4, pPoints->Get_Count());
248 
249 	int	n	= 0;
250 
251 	for(int i=0; i<pPoints->Get_Count(); i++)
252 	{
253 		CSG_Shape	*pPoint	= pPoints->Get_Shape(i);
254 
255 		if( !pPoint->is_NoData(Field) )
256 		{
257 			m_Points[n][0]	= pPoint->Get_Point(0).x;
258 			m_Points[n][1]	= pPoint->Get_Point(0).y;
259 			m_Points[n][2]	= zScale * (zField < 0 ? pPoint->Get_Z(0) : pPoint->asDouble(zField));
260 			m_Points[n][3]	= pPoint->asDouble(Field);
261 
262 			n++;
263 		}
264 	}
265 
266 	if( n < 1 )
267 	{
268 		Error_Set(_TL("no valid points in data set"));
269 
270 		return( false );
271 	}
272 
273 	m_Points.Set_Rows(n);	// resize if there are no-data values
274 
275 	//-----------------------------------------------------
276 	if( !m_Search_Options.Do_Use_All(true) && !m_Search.Create(m_Points) )
277 	{
278 		Error_Set(_TL("failed to initialize search engine"));
279 
280 		return( false );
281 	}
282 
283 	return( true );
284 }
285 
286 //---------------------------------------------------------
Finalize(void)287 bool CGridding3D_IDW::Finalize(void)
288 {
289 	m_Search.Destroy();
290 	m_Points.Destroy();
291 
292 	return( true );
293 }
294 
295 
296 ///////////////////////////////////////////////////////////
297 //														 //
298 ///////////////////////////////////////////////////////////
299 
300 //---------------------------------------------------------
Get_Value(double Coordinate[3],double & Value)301 bool CGridding3D_IDW::Get_Value(double Coordinate[3], double &Value)
302 {
303 	CSG_Array_Int	Index;	CSG_Vector	Distance;
304 
305 	if( m_Search.is_Okay() )
306 	{
307 		if( m_Search.Get_Nearest_Points(Coordinate,
308 			m_Search_Options.Get_Max_Points(),
309 			m_Search_Options.Get_Radius(), Index, Distance
310 		) < m_Search_Options.Get_Min_Points() )
311 		{
312 			return( false );
313 		}
314 	}
315 
316 	CSG_Simple_Statistics	s;
317 
318 	int	nPoints	= m_Search.is_Okay() ? (int)Index.Get_Size() : m_Points.Get_NRows();
319 
320 	for(int i=0; i<nPoints; i++)
321 	{
322 		double	*p	= m_Points[m_Search.is_Okay() ? Index[i] : i];
323 		double	 d	= m_Search.is_Okay() ? Distance[i] : Get_Distance(Coordinate, p);
324 
325 		if( d > 0. )
326 		{
327 			s.Add_Value(p[3], m_Weighting.Get_Weight(d));
328 		}
329 		else	// d == 0! there is a point at the requested coordinate!
330 		{
331 			s.Create();	s	+= p[3];
332 
333 			for(++i; i<nPoints; i++)	// is there more than one point?!
334 			{
335 				p	= m_Points[m_Search.is_Okay() ? Index[i] : i];
336 				d	= m_Search.is_Okay() ? Distance[i] : Get_Distance(Coordinate, p);
337 
338 				if( d <= 0. )
339 				{
340 					s	+= p[3];
341 				}
342 			}
343 		}
344 	}
345 
346 	Value	= s.Get_Mean();
347 
348 	return( true );
349 }
350 
351 //---------------------------------------------------------
Get_Distance(double Coordinate[3],double Point[3])352 inline double CGridding3D_IDW::Get_Distance(double Coordinate[3], double Point[3])
353 {
354 	double	dx	= Coordinate[0] - Point[0];
355 	double	dy	= Coordinate[1] - Point[1];
356 	double	dz	= Coordinate[2] - Point[2];
357 
358 	return( sqrt(dx*dx + dy*dy + dz*dz) );
359 }
360 
361 //---------------------------------------------------------
is_Identical(double Coordinate[3],double Point[3])362 inline bool CGridding3D_IDW::is_Identical(double Coordinate[3], double Point[3])
363 {
364 	return( Coordinate[0] == Point[0]
365 		&&  Coordinate[1] == Point[1]
366 		&&  Coordinate[2] == Point[2]
367 	);
368 }
369 
370 
371 ///////////////////////////////////////////////////////////
372 //														 //
373 //														 //
374 //														 //
375 ///////////////////////////////////////////////////////////
376 
377 //---------------------------------------------------------
378