1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Tool Library //
9 // statistics_kriging //
10 // //
11 //-------------------------------------------------------//
12 // //
13 // kriging_simple.cpp //
14 // //
15 // Olaf Conrad (C) 2015 //
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 "kriging_simple.h"
49
50
51 ///////////////////////////////////////////////////////////
52 // //
53 // //
54 // //
55 ///////////////////////////////////////////////////////////
56
57 //---------------------------------------------------------
CKriging_Simple(void)58 CKriging_Simple::CKriging_Simple(void)
59 {
60 //-----------------------------------------------------
61 Set_Name (_TL("Simple Kriging"));
62
63 Set_Author ("O.Conrad (c) 2015");
64
65 Set_Description (_TW(
66 "Simple Kriging for grid interpolation from irregular sample points."
67 ));
68
69 //-----------------------------------------------------
70 }
71
72
73 ///////////////////////////////////////////////////////////
74 // //
75 ///////////////////////////////////////////////////////////
76
77 //---------------------------------------------------------
Get_Weights(const CSG_Matrix & Points,CSG_Matrix & W)78 bool CKriging_Simple::Get_Weights(const CSG_Matrix &Points, CSG_Matrix &W)
79 {
80 int n = Points.Get_NRows();
81
82 if( n < 1 || !W.Create(n, n) )
83 {
84 return( false );
85 }
86
87 for(int i=0; i<n; i++)
88 {
89 W[i][i] = 0.; // diagonal...
90
91 for(int j=i+1; j<n; j++)
92 {
93 W[i][j] = W[j][i] = Get_Weight(Points[i], Points[j]);
94 }
95 }
96
97 return( W.Set_Inverse(m_Search.is_Okay(), n) );
98 }
99
100
101 ///////////////////////////////////////////////////////////
102 // //
103 ///////////////////////////////////////////////////////////
104
105 //---------------------------------------------------------
Get_Value(double x,double y,double & v,double & e)106 bool CKriging_Simple::Get_Value(double x, double y, double &v, double &e)
107 {
108 CSG_Matrix __Points, __W; double **P, **W; int i, n = 0;
109
110 if( !m_Search.is_Okay() ) // global
111 {
112 n = m_Points.Get_NRows();
113 P = m_Points.Get_Data ();
114 W = m_W .Get_Data ();
115 }
116 else if( Get_Points(x, y, __Points) && Get_Weights(__Points, __W) ) // local
117 {
118 n = __Points.Get_NRows();
119 P = __Points.Get_Data ();
120 W = __W .Get_Data ();
121 }
122
123 if( n < 1 )
124 {
125 return( false );
126 }
127
128 //-----------------------------------------------------
129 CSG_Vector G(n);
130
131 for(i=0; i<n; i++)
132 {
133 G[i] = Get_Weight(x, y, P[i][0], P[i][1]);
134 }
135
136 for(i=0, v=0., e=0.; i<n; i++)
137 {
138 double Lambda = 0.;
139
140 for(int j=0; j<n; j++)
141 {
142 Lambda += W[i][j] * G[j];
143 }
144
145 v += Lambda * P[i][2];
146 e += Lambda * G[i];
147 }
148
149 //-----------------------------------------------------
150 return( true );
151 }
152
153
154 ///////////////////////////////////////////////////////////
155 // //
156 // //
157 // //
158 ///////////////////////////////////////////////////////////
159
160 //---------------------------------------------------------
161