1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Tool Library //
9 // statistics_points //
10 // //
11 //-------------------------------------------------------//
12 // //
13 // GSPoints_Variogram_Surface.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 "GSPoints_Variogram_Surface.h"
50
51
52 ///////////////////////////////////////////////////////////
53 // //
54 // //
55 // //
56 ///////////////////////////////////////////////////////////
57
58 //---------------------------------------------------------
CGSPoints_Variogram_Surface(void)59 CGSPoints_Variogram_Surface::CGSPoints_Variogram_Surface(void)
60 {
61 Set_Name (_TL("Variogram Surface"));
62
63 Set_Author ("O.Conrad (c) 2010");
64
65 Set_Description(
66 _TL("Calculates a variogram surface.")
67 );
68
69 //-----------------------------------------------------
70 Parameters.Add_Shapes("",
71 "POINTS" , _TL("Points"),
72 _TL(""),
73 PARAMETER_INPUT, SHAPE_TYPE_Point
74 );
75
76 Parameters.Add_Table_Field("POINTS",
77 "FIELD" , _TL("Attribute"),
78 _TL("")
79 );
80
81 Parameters.Add_Int("",
82 "DISTCOUNT" , _TL("Number of Distance Classes"),
83 _TL(""),
84 10, 1, true
85 );
86
87 Parameters.Add_Int("",
88 "NSKIP" , _TL("Skip Number"),
89 _TL(""),
90 1, 1, true
91 );
92
93 Parameters.Add_Grid_Output("", "COUNT" , _TL("Number of Pairs" ), _TL(""));
94 Parameters.Add_Grid_Output("", "VARIANCE" , _TL("Variogram Surface" ), _TL(""));
95 Parameters.Add_Grid_Output("", "COVARIANCE", _TL("Covariance Surface"), _TL(""));
96 }
97
98
99 ///////////////////////////////////////////////////////////
100 // //
101 ///////////////////////////////////////////////////////////
102
103 //---------------------------------------------------------
On_Execute(void)104 bool CGSPoints_Variogram_Surface::On_Execute(void)
105 {
106 CSG_Shapes *pPoints = Parameters("POINTS")->asShapes();
107
108 int Field = Parameters("FIELD" )->asInt();
109 int nSkip = Parameters("NSKIP" )->asInt();
110 int nDistances = Parameters("DISTCOUNT")->asInt();
111
112 double Lag = pPoints->Get_Extent().Get_XRange() < pPoints->Get_Extent().Get_YRange()
113 ? pPoints->Get_Extent().Get_XRange() / nDistances
114 : pPoints->Get_Extent().Get_YRange() / nDistances;
115
116 int nx = 1 + (int)(pPoints->Get_Extent().Get_XRange() / Lag);
117 int ny = 1 + (int)(pPoints->Get_Extent().Get_YRange() / Lag);
118
119 CSG_Grid *pCount = SG_Create_Grid(SG_DATATYPE_Int , 1 + 2 * nx, 1 + 2 * ny, Lag, -nx * Lag, -ny * Lag);
120 CSG_Grid *pVariance = SG_Create_Grid(SG_DATATYPE_Float, 1 + 2 * nx, 1 + 2 * ny, Lag, -nx * Lag, -ny * Lag);
121 CSG_Grid *pCovariance = SG_Create_Grid(SG_DATATYPE_Float, 1 + 2 * nx, 1 + 2 * ny, Lag, -nx * Lag, -ny * Lag);
122
123 pCount ->Fmt_Name("%s [%s]" , pPoints->Get_Name(), _TL("Count" ));
124 pVariance ->Fmt_Name("%s [%s: %s]", pPoints->Get_Name(), _TL("Variogram Surface" ), pPoints->Get_Field_Name(Field));
125 pCovariance ->Fmt_Name("%s [%s: %s]", pPoints->Get_Name(), _TL("Covariance Surface"), pPoints->Get_Field_Name(Field));
126
127 Parameters("COUNT" )->Set_Value(pCount );
128 Parameters("VARIANCE" )->Set_Value(pVariance );
129 Parameters("COVARIANCE")->Set_Value(pCovariance);
130
131 //-----------------------------------------------------
132 for(int i=0, n=0; i<pPoints->Get_Count() && Set_Progress(n, 0.5 * SG_Get_Square(pPoints->Get_Count() / nSkip)); i+=nSkip)
133 {
134 CSG_Shape *pPoint = pPoints->Get_Shape(i);
135
136 if( !pPoint->is_NoData(Field) )
137 {
138 TSG_Point pi = pPoint->Get_Point(0);
139 double zi = pPoint->asDouble(Field);
140
141 for(int j=i+nSkip; j<pPoints->Get_Count(); j+=nSkip, n++)
142 {
143 pPoint = pPoints->Get_Shape(j);
144
145 if( !pPoint->is_NoData(Field) )
146 {
147 TSG_Point pj = pPoint->Get_Point(0);
148 double zj = pPoint->asDouble(Field);
149
150 double v = SG_Get_Square(zi - zj);
151 double c = (zi - pPoints->Get_Mean(Field)) * (zj - pPoints->Get_Mean(Field));
152
153 pj.x = (pi.x - pj.x) / Lag;
154 pj.y = (pi.y - pj.y) / Lag;
155
156 int x = (int)(pj.x + (pj.x > 0. ? 0.5 : -0.5));
157 int y = (int)(pj.y + (pj.y > 0. ? 0.5 : -0.5));
158
159 pCount ->Add_Value(nx + x, ny + y, 1);
160 pCount ->Add_Value(nx - x, ny - y, 1);
161 pVariance ->Add_Value(nx + x, ny + y, v);
162 pVariance ->Add_Value(nx - x, ny - y, v);
163 pCovariance->Add_Value(nx + x, ny + y, c);
164 pCovariance->Add_Value(nx - x, ny - y, c);
165 }
166 }
167 }
168 }
169
170 //-----------------------------------------------------
171 for(sLong iCell=0; iCell<pCount->Get_NCells(); iCell++)
172 {
173 if( pCount->asInt(iCell) > 0 )
174 {
175 pVariance ->Mul_Value(iCell, 0.5 / pCount->asDouble(iCell));
176 pCovariance->Mul_Value(iCell, 1.0 / pCount->asDouble(iCell));
177 }
178 else
179 {
180 pVariance ->Set_NoData(iCell);
181 pCovariance->Set_NoData(iCell);
182 }
183 }
184
185 //-----------------------------------------------------
186 return( true );
187 }
188
189
190 ///////////////////////////////////////////////////////////
191 // //
192 // //
193 // //
194 ///////////////////////////////////////////////////////////
195
196 //---------------------------------------------------------
197