1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Tool Library //
9 // Grid_Gridding //
10 // //
11 //-------------------------------------------------------//
12 // //
13 // Interpolation.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 "Interpolation.h"
52
53
54 ///////////////////////////////////////////////////////////
55 // //
56 // //
57 // //
58 ///////////////////////////////////////////////////////////
59
60 //---------------------------------------------------------
CInterpolation(bool bCrossValidation,bool bMultiThreading)61 CInterpolation::CInterpolation(bool bCrossValidation, bool bMultiThreading)
62 {
63 m_bMultiThreading = bMultiThreading;
64
65 //-----------------------------------------------------
66 Parameters.Add_Shapes("",
67 "POINTS" , _TL("Points"),
68 _TL(""),
69 PARAMETER_INPUT, SHAPE_TYPE_Point
70 );
71
72 Parameters.Add_Table_Field("POINTS",
73 "FIELD" , _TL("Attribute"),
74 _TL("")
75 );
76
77 //-----------------------------------------------------
78 if( bCrossValidation )
79 {
80 Parameters.Add_Choice("",
81 "CV_METHOD" , _TL("Cross Validation"),
82 _TL(""),
83 CSG_String::Format("%s|%s|%s|%s",
84 _TL("none"),
85 _TL("leave one out"),
86 _TL("2-fold"),
87 _TL("k-fold")
88 ), 0
89 );
90
91 Parameters.Add_Table("CV_METHOD",
92 "CV_SUMMARY" , _TL("Cross Validation Summary"),
93 _TL(""),
94 PARAMETER_OUTPUT_OPTIONAL
95 );
96
97 Parameters.Add_Shapes("CV_METHOD",
98 "CV_RESIDUALS" , _TL("Cross Validation Residuals"),
99 _TL(""),
100 PARAMETER_OUTPUT_OPTIONAL, SHAPE_TYPE_Point
101 );
102
103 Parameters.Add_Int("CV_METHOD",
104 "CV_SAMPLES" , _TL("Cross Validation Subsamples"),
105 _TL("number of subsamples for k-fold cross validation"),
106 10, 2, true
107 );
108 }
109
110 //-----------------------------------------------------
111 m_Grid_Target.Create(&Parameters, true, "", "TARGET_");
112 }
113
114
115 ///////////////////////////////////////////////////////////
116 // //
117 ///////////////////////////////////////////////////////////
118
119 //---------------------------------------------------------
On_Parameter_Changed(CSG_Parameters * pParameters,CSG_Parameter * pParameter)120 int CInterpolation::On_Parameter_Changed(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
121 {
122 if( pParameter->Cmp_Identifier("POINTS") )
123 {
124 m_Grid_Target.Set_User_Defined(pParameters, pParameter->asShapes());
125 }
126
127 m_Grid_Target.On_Parameter_Changed(pParameters, pParameter);
128
129 return( CSG_Tool::On_Parameter_Changed(pParameters, pParameter) );
130 }
131
132 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)133 int CInterpolation::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
134 {
135 if( pParameter->Cmp_Identifier("CV_METHOD") )
136 {
137 pParameters->Set_Enabled("CV_SUMMARY" , pParameter->asInt() != 0); // !none
138 pParameters->Set_Enabled("CV_RESIDUALS", pParameter->asInt() == 1); // leave one out
139 pParameters->Set_Enabled("CV_SAMPLES" , pParameter->asInt() == 3); // k-fold
140 }
141
142 m_Grid_Target.On_Parameters_Enable(pParameters, pParameter);
143
144 return( CSG_Tool::On_Parameters_Enable(pParameters, pParameter) );
145 }
146
147
148 ///////////////////////////////////////////////////////////
149 // //
150 ///////////////////////////////////////////////////////////
151
152 //---------------------------------------------------------
On_Execute(void)153 bool CInterpolation::On_Execute(void)
154 {
155 //-----------------------------------------------------
156 m_pPoints = Parameters("POINTS")->asShapes();
157 m_zField = Parameters("FIELD" )->asInt ();
158
159 if( m_pPoints->Get_Count() < 2 )
160 {
161 Error_Set(_TL("not enough points"));
162
163 return( false );
164 }
165
166 //-----------------------------------------------------
167 m_pGrid = m_Grid_Target.Get_Grid();
168
169 if( m_pGrid == NULL )
170 {
171 return( false );
172 }
173
174 m_pGrid->Fmt_Name("%s.%s [%s]", m_pPoints->Get_Name(), Parameters("FIELD")->asString(), Get_Name().c_str());
175
176 //-----------------------------------------------------
177 if( !Interpolate() )
178 {
179 return( false );
180 }
181
182 _Get_Cross_Validation();
183
184 //-----------------------------------------------------
185 return( true );
186 }
187
188
189 ///////////////////////////////////////////////////////////
190 // //
191 ///////////////////////////////////////////////////////////
192
193 //---------------------------------------------------------
_Interpolate(int x,int y)194 inline void CInterpolation::_Interpolate(int x, int y)
195 {
196 double z;
197
198 if( Get_Value(m_pGrid->Get_System().Get_Grid_to_World(x, y), z) )
199 {
200 m_pGrid->Set_Value(x, y, z);
201 }
202 else
203 {
204 m_pGrid->Set_NoData(x, y);
205 }
206 }
207
208 //---------------------------------------------------------
Interpolate(void)209 bool CInterpolation::Interpolate(void)
210 {
211 if( !On_Initialize() )
212 {
213 return( false );
214 }
215
216 //-----------------------------------------------------
217 for(int y=0; y<m_pGrid->Get_NY() && Set_Progress(y, m_pGrid->Get_NY()); y++)
218 {
219 if( m_bMultiThreading )
220 {
221 #pragma omp parallel for
222 for(int x=0; x<m_pGrid->Get_NX(); x++)
223 {
224 _Interpolate(x, y);
225 }
226 }
227 else
228 {
229 for(int x=0; x<m_pGrid->Get_NX(); x++)
230 {
231 _Interpolate(x, y);
232 }
233 }
234 }
235
236 //-----------------------------------------------------
237 On_Finalize();
238
239 return( true );
240 }
241
242
243 ///////////////////////////////////////////////////////////
244 // //
245 ///////////////////////////////////////////////////////////
246
247 //---------------------------------------------------------
_Get_Cross_Validation(void)248 bool CInterpolation::_Get_Cross_Validation(void)
249 {
250 if( !Parameters("CV_METHOD") )
251 {
252 return( true );
253 }
254
255 //-----------------------------------------------------
256 int nSubSets;
257
258 switch( Parameters("CV_METHOD")->asInt() )
259 {
260 default: // none
261 return( true );
262
263 case 1: // leave one out (LOOVC)
264 nSubSets = 1;
265 break;
266
267 case 2: // 2-fold
268 nSubSets = 2;
269 break;
270
271 case 3: // k-fold
272 nSubSets = Parameters("CV_SAMPLES")->asInt();
273 break;
274 }
275
276 //-----------------------------------------------------
277 Process_Set_Text(_TL("Cross Validation"));
278
279 CSG_Simple_Statistics SFull, SR, SE;
280
281 CSG_Shapes *pFull = m_pPoints;
282
283 int i, nSamples = 0;
284
285 for(i=0; i<pFull->Get_Count(); i++)
286 {
287 CSG_Shape *pPoint = pFull->Get_Shape(i);
288
289 if( !pPoint->is_NoData(m_zField) )
290 {
291 SFull += pPoint->asDouble(m_zField);
292 }
293 }
294
295 //-----------------------------------------------------
296 // leave-one-out cross validation (LOOCV)
297
298 if( nSubSets <= 1 || nSubSets > pFull->Get_Count() / 2 )
299 {
300 CSG_Shapes Sample(*pFull), *pResiduals; m_pPoints = &Sample;
301
302 if( (pResiduals = Parameters("CV_RESIDUALS")->asShapes()) != NULL )
303 {
304 pResiduals->Create(SHAPE_TYPE_Point, CSG_String::Format("%s [%s, %s]", m_pPoints->Get_Name(), Get_Name().c_str(), _TL("Residuals")));
305 pResiduals->Add_Field(pFull->Get_Field_Name(m_zField), SG_DATATYPE_Double);
306 pResiduals->Add_Field("PREDICTED", SG_DATATYPE_Double);
307 pResiduals->Add_Field("RESIDUALS" , SG_DATATYPE_Double);
308 }
309
310 for(i=pFull->Get_Count()-1; i>=0 && Set_Progress(pFull->Get_Count()-1-i, pFull->Get_Count()); i--)
311 {
312 CSG_Shape *pPoint = pFull->Get_Shape(i);
313
314 Sample.Del_Shape(i);
315
316 double z;
317
318 if( !pPoint->is_NoData(m_zField) && On_Initialize() && Get_Value(pPoint->Get_Point(0), z) )
319 {
320 nSamples++;
321
322 SE += SG_Get_Square(z - pPoint->asDouble(m_zField));
323 SR += SG_Get_Square(z - (SFull.Get_Sum() - pPoint->asDouble(m_zField)) / Sample.Get_Count());
324
325 if( pResiduals )
326 {
327 CSG_Shape *pResidual = pResiduals->Add_Shape();
328
329 pResidual->Add_Point(pPoint->Get_Point(0));
330 pResidual->Set_Value(0, pPoint->asDouble(m_zField));
331 pResidual->Set_Value(1, z);
332 pResidual->Set_Value(2, pPoint->asDouble(m_zField) - z);
333 }
334 }
335
336 Sample.Add_Shape(pPoint);
337 }
338 }
339
340 //-----------------------------------------------------
341 // k-fold cross validation
342
343 else
344 {
345 CSG_Array_Int SubSet(pFull->Get_Count());
346
347 for(i=0; i<pFull->Get_Count(); i++)
348 {
349 SubSet[i] = i % nSubSets;
350 }
351
352 //-------------------------------------------------
353 for(int iSubSet=0; iSubSet<nSubSets && Process_Get_Okay(); iSubSet++)
354 {
355 CSG_Simple_Statistics iSFull;
356
357 CSG_Shapes Sample[2];
358
359 Sample[0].Create(SHAPE_TYPE_Point, SG_T(""), pFull); m_pPoints = &Sample[0];
360 Sample[1].Create(SHAPE_TYPE_Point, SG_T(""), pFull);
361
362 for(i=0; i<pFull->Get_Count(); i++)
363 {
364 CSG_Shape *pPoint = pFull->Get_Shape(i);
365
366 if( !pPoint->is_NoData(m_zField) )
367 {
368 if( SubSet[i] == iSubSet )
369 {
370 Sample[1].Add_Shape(pPoint);
371 }
372 else
373 {
374 Sample[0].Add_Shape(pPoint);
375
376 iSFull += pPoint->asDouble(m_zField);
377 }
378 }
379 }
380
381 //---------------------------------------------
382 if( On_Initialize() )
383 {
384 nSamples++;
385
386 for(i=0; i<Sample[1].Get_Count(); i++)
387 {
388 CSG_Shape *pPoint = Sample[1].Get_Shape(i);
389
390 double z;
391
392 if( Get_Value(pPoint->Get_Point(0), z) )
393 {
394 SE += SG_Get_Square(z - pPoint->asDouble(m_zField));
395 SR += SG_Get_Square(z - iSFull.Get_Mean());
396 }
397 }
398 }
399 }
400 }
401
402 //-----------------------------------------------------
403 if( nSamples == 0 )
404 {
405 return( false );
406 }
407
408 //-----------------------------------------------------
409 CSG_Table *pSummary = Parameters("CV_SUMMARY")->asTable();
410
411 if( pSummary )
412 {
413 pSummary->Destroy();
414 pSummary->Set_Name(_TL("Cross Validation"));
415
416 pSummary->Add_Field(_TL("Parameter"), SG_DATATYPE_String);
417 pSummary->Add_Field(_TL("Value" ), SG_DATATYPE_Double);
418
419 #define CV_ADD_SUMMARY(name, value) { CSG_Table_Record *pR = pSummary->Add_Record();\
420 pR->Set_Value(0, name);\
421 pR->Set_Value(1, value);\
422 }
423
424 CV_ADD_SUMMARY("SAMPLES", nSamples);
425 CV_ADD_SUMMARY("MSE" , SE.Get_Mean());
426 CV_ADD_SUMMARY("RMSE" , sqrt(SE.Get_Mean()));
427 CV_ADD_SUMMARY("NRMSE" , sqrt(SE.Get_Mean()) / SFull.Get_Range() * 100.0);
428 CV_ADD_SUMMARY("R2" , SR.Get_Sum() / (SR.Get_Sum() + SE.Get_Sum()) * 100.0);
429 }
430
431 //-----------------------------------------------------
432 CSG_String Message;
433
434 Message += CSG_String::Format("\n%s:\n" , _TL("Cross Validation"));
435 Message += CSG_String::Format("\t%s:\t%s\n" , _TL("Method" ), Parameters("CV_METHOD")->asString());
436 Message += CSG_String::Format("\t%s:\t%d\n" , _TL("Samples"), nSamples);
437 Message += CSG_String::Format("\t%s:\t%f\n" , _TL("MSE" ), SE.Get_Mean());
438 Message += CSG_String::Format("\t%s:\t%f\n" , _TL("RMSE" ), sqrt(SE.Get_Mean()));
439 Message += CSG_String::Format("\t%s:\t%.2f\n", _TL("NRMSE" ), sqrt(SE.Get_Mean()) / SFull.Get_Range() * 100.0);
440 Message += CSG_String::Format("\t%s:\t%.2f\n", _TL("R2" ), SR.Get_Sum() / (SR.Get_Sum() + SE.Get_Sum()) * 100.0);
441
442 Message_Add(Message, false);
443
444 //-----------------------------------------------------
445 return( true );
446 }
447
448
449 ///////////////////////////////////////////////////////////
450 // //
451 // //
452 // //
453 ///////////////////////////////////////////////////////////
454
455 //---------------------------------------------------------
456