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