root / org.gvsig.toolbox / trunk / org.gvsig.toolbox / org.gvsig.toolbox.algorithm / src / main / java / es / unex / sextante / rasterize / kriging / KrigingAlgorithm.java @ 59
History | View | Annotate | Download (8.7 KB)
1 |
package es.unex.sextante.rasterize.kriging; |
---|---|
2 |
|
3 |
import java.awt.geom.Point2D; |
4 |
import java.util.Arrays; |
5 |
|
6 |
import Jama.Matrix; |
7 |
import es.unex.sextante.additionalInfo.AdditionalInfoNumericalValue; |
8 |
import es.unex.sextante.closestpts.Point3D; |
9 |
import es.unex.sextante.closestpts.PtAndDistance; |
10 |
import es.unex.sextante.core.Sextante; |
11 |
import es.unex.sextante.dataObjects.IRasterLayer; |
12 |
import es.unex.sextante.exceptions.GeoAlgorithmExecutionException; |
13 |
import es.unex.sextante.exceptions.RepeatedParameterNameException; |
14 |
import es.unex.sextante.rasterWrappers.GridCell; |
15 |
import es.unex.sextante.rasterize.interpolationBase.BaseInterpolationAlgorithm; |
16 |
|
17 |
public class KrigingAlgorithm |
18 |
extends
|
19 |
BaseInterpolationAlgorithm { |
20 |
|
21 |
public static final String SILL = "SILL"; |
22 |
public static final String RANGE = "RANGE"; |
23 |
public static final String NUGGET = "NUGGET"; |
24 |
public static final String MODEL = "MODEL"; |
25 |
public static final String MINPOINTS = "MINPOINTS"; |
26 |
public static final String MAXPOINTS = "MAXPOINTS"; |
27 |
public static final String VARIANCE = "VARIANCE"; |
28 |
public static final String CROSSVALIDATION = "CROSSVALIDATION"; |
29 |
|
30 |
private double m_dNugget; |
31 |
private double m_dScale; |
32 |
private double m_dRange; |
33 |
private double m_dGammas[]; |
34 |
private int m_iMinPoints; |
35 |
private int m_iMaxPoints; |
36 |
private int m_iModel; |
37 |
private boolean m_bCreateVarianceLayer; |
38 |
private double m_dWeights[][]; |
39 |
private Matrix m_Matrix;
|
40 |
private IRasterLayer m_Variance;
|
41 |
|
42 |
|
43 |
@Override
|
44 |
public void defineCharacteristics() { |
45 |
|
46 |
super.defineCharacteristics();
|
47 |
|
48 |
setGroup(Sextante.getText("Rasterization_and_interpolation"));
|
49 |
setName(Sextante.getText("Kriging"));
|
50 |
|
51 |
final String sModels[] = { Sextante.getText("Spherical"), Sextante.getText("Exponential"), Sextante.getText("Gaussian") }; |
52 |
|
53 |
try {
|
54 |
m_Parameters.addNumericalValue(MINPOINTS, Sextante.getText("Mino_number_of_points"),
|
55 |
AdditionalInfoNumericalValue.NUMERICAL_VALUE_INTEGER, 4, 1, Integer.MAX_VALUE); |
56 |
m_Parameters.addNumericalValue(MAXPOINTS, Sextante.getText("Max_number_of_points"),
|
57 |
AdditionalInfoNumericalValue.NUMERICAL_VALUE_INTEGER, 25, 1, Integer.MAX_VALUE); |
58 |
m_Parameters.addSelection(MODEL, Sextante.getText("Model"), sModels);
|
59 |
m_Parameters.addNumericalValue(NUGGET, Sextante.getText("Nugget"), AdditionalInfoNumericalValue.NUMERICAL_VALUE_DOUBLE,
|
60 |
0, 0, Double.MAX_VALUE); |
61 |
m_Parameters.addNumericalValue(SILL, Sextante.getText("Sill"), AdditionalInfoNumericalValue.NUMERICAL_VALUE_DOUBLE, 10, |
62 |
0, Double.MAX_VALUE); |
63 |
m_Parameters.addNumericalValue(RANGE, Sextante.getText("Range"), AdditionalInfoNumericalValue.NUMERICAL_VALUE_DOUBLE,
|
64 |
100, 0, Double.MAX_VALUE); |
65 |
addOutputTable(CROSSVALIDATION, Sextante.getText("Cross_validation"));
|
66 |
addOutputRasterLayer("VARIANCE", Sextante.getText("Variances")); |
67 |
} |
68 |
catch (final RepeatedParameterNameException e) { |
69 |
Sextante.addErrorToLog(e); |
70 |
} |
71 |
|
72 |
} |
73 |
|
74 |
|
75 |
@Override
|
76 |
protected void setValues() throws GeoAlgorithmExecutionException { |
77 |
|
78 |
super.setValues();
|
79 |
|
80 |
m_iMaxPoints = m_Parameters.getParameterValueAsInt(MAXPOINTS); |
81 |
m_iMinPoints = m_Parameters.getParameterValueAsInt(MINPOINTS); |
82 |
m_iModel = m_Parameters.getParameterValueAsInt(MODEL); |
83 |
m_dNugget = m_Parameters.getParameterValueAsDouble(NUGGET); |
84 |
m_dScale = m_Parameters.getParameterValueAsDouble(SILL) - m_dNugget; |
85 |
m_dRange = m_Parameters.getParameterValueAsDouble(RANGE); |
86 |
m_bCreateVarianceLayer = true; //m_Parameters.getParameterValueAsBoolean("VARIANCE"); |
87 |
m_dWeights = new double[m_iMaxPoints + 1][m_iMaxPoints + 1]; |
88 |
m_Matrix = new Matrix(m_dWeights);
|
89 |
m_dGammas = new double[m_iMaxPoints + 1]; |
90 |
if (m_bCreateVarianceLayer) {
|
91 |
m_Variance = getNewRasterLayer(VARIANCE, m_Layer.getName() + Sextante.getText("[variances]"),
|
92 |
IRasterLayer.RASTER_DATA_TYPE_DOUBLE); |
93 |
} |
94 |
|
95 |
|
96 |
} |
97 |
|
98 |
|
99 |
@Override
|
100 |
protected double getValueAt(final int x, |
101 |
final int y) { |
102 |
|
103 |
final Point2D pt = m_AnalysisExtent.getWorldCoordsFromGridCoords(new GridCell(x, y, 0)); |
104 |
final PtAndDistance[] nearestPoints = m_SearchEngine.getClosestPoints(pt.getX(), pt.getY(), m_dDistance); |
105 |
|
106 |
Arrays.sort(nearestPoints);
|
107 |
|
108 |
final int n = Math.min(nearestPoints.length, m_iMaxPoints); |
109 |
m_NearestPoints = new PtAndDistance[n];
|
110 |
System.arraycopy(nearestPoints, 0, m_NearestPoints, 0, n); |
111 |
|
112 |
final double dValue = interpolate(pt.getX(), pt.getY()); |
113 |
|
114 |
return dValue;
|
115 |
|
116 |
} |
117 |
|
118 |
|
119 |
@Override
|
120 |
protected double getValueAt(final double x, |
121 |
final double y) { |
122 |
|
123 |
try {
|
124 |
final PtAndDistance[] nearestPoints = m_SearchEngine.getClosestPoints(x, y, m_dDistance); |
125 |
Arrays.sort(nearestPoints);
|
126 |
final int n = Math.min(nearestPoints.length, m_iMaxPoints) - 1; |
127 |
m_NearestPoints = new PtAndDistance[n];
|
128 |
System.arraycopy(nearestPoints, 1, m_NearestPoints, 0, n); |
129 |
|
130 |
return interpolate(x, y);
|
131 |
} |
132 |
catch (final Exception e) { |
133 |
return NO_DATA;
|
134 |
} |
135 |
|
136 |
} |
137 |
|
138 |
|
139 |
@Override
|
140 |
protected double interpolate(final double x, |
141 |
final double y) { |
142 |
|
143 |
int i, j, nPoints;
|
144 |
double dLambda, dValue, dVariance;
|
145 |
Point3D pt; |
146 |
|
147 |
if ((nPoints = getWeights(x, y)) >= m_iMinPoints) {
|
148 |
for (i = 0; i < nPoints; i++) { |
149 |
m_dGammas[i] = getWeight(m_NearestPoints[i].getDist()); |
150 |
} |
151 |
|
152 |
m_dGammas[nPoints] = 1.0;
|
153 |
|
154 |
for (i = 0, dValue = 0.0, dVariance = 0.0; i < nPoints; i++) { |
155 |
pt = m_NearestPoints[i].getPt(); |
156 |
for (j = 0, dLambda = 0.0; j <= nPoints; j++) { |
157 |
dLambda += m_dWeights[i][j] * m_dGammas[j]; |
158 |
} |
159 |
|
160 |
dValue += dLambda * pt.getZ(); |
161 |
|
162 |
if (m_bCreateVarianceLayer) {
|
163 |
dVariance += dLambda * m_dGammas[i]; |
164 |
} |
165 |
} |
166 |
|
167 |
if (m_bCreateVarianceLayer) {
|
168 |
final GridCell cell = m_AnalysisExtent.getGridCoordsFromWorldCoords(x, y);
|
169 |
m_Variance.setCellValue(cell.getX(), cell.getY(), dVariance); |
170 |
} |
171 |
|
172 |
return dValue;
|
173 |
} |
174 |
|
175 |
if (m_bCreateVarianceLayer) {
|
176 |
final GridCell cell = m_AnalysisExtent.getGridCoordsFromWorldCoords(x, y);
|
177 |
m_Variance.setNoData(cell.getX(), cell.getY()); |
178 |
} |
179 |
|
180 |
return NO_DATA;
|
181 |
} |
182 |
|
183 |
|
184 |
private double getWeight(double d) { |
185 |
|
186 |
if (d == 0.0) { |
187 |
d = 0.0001;
|
188 |
} |
189 |
|
190 |
switch (m_iModel) {
|
191 |
case 0: // Spherical Model |
192 |
if (d >= m_dRange) {
|
193 |
d = m_dNugget + m_dScale; |
194 |
} |
195 |
else {
|
196 |
d = m_dNugget + m_dScale * (3 * d / (2 * m_dRange) - d * d * d / (2 * m_dRange * m_dRange * m_dRange)); |
197 |
} |
198 |
break;
|
199 |
|
200 |
case 1: // Exponential Model |
201 |
d = m_dNugget + m_dScale * (1 - Math.exp(-3 * d / m_dRange)); |
202 |
break;
|
203 |
|
204 |
case 2: // Gaussian Model |
205 |
d = 1 - Math.exp(-3 * d / (m_dRange * m_dRange)); |
206 |
d = m_dNugget + m_dScale * d * d; |
207 |
break;
|
208 |
} |
209 |
|
210 |
return (d);
|
211 |
} |
212 |
|
213 |
|
214 |
private int getWeights(final double x, |
215 |
final double y) { |
216 |
|
217 |
int i, j, n;
|
218 |
double dx, dy;
|
219 |
Point3D pt, pt2; |
220 |
|
221 |
m_dWeights = m_Matrix.getArray(); |
222 |
|
223 |
if ((n = Math.min(m_NearestPoints.length, m_iMaxPoints)) >= m_iMinPoints) { |
224 |
//n = Math.min(n, m_iMaxPoints);
|
225 |
for (i = 0; i < n; i++) { |
226 |
pt = m_NearestPoints[i].getPt(); |
227 |
m_dWeights[i][i] = 0.0;
|
228 |
m_dWeights[i][n] = m_dWeights[n][i] = 1.0;
|
229 |
for (j = i + 1; j < n; j++) { |
230 |
pt2 = m_NearestPoints[j].getPt(); |
231 |
dx = pt.getX() - pt2.getX(); |
232 |
dy = pt.getY() - pt2.getY(); |
233 |
m_dWeights[i][j] = m_dWeights[j][i] = getWeight(Math.sqrt(dx * dx + dy * dy));
|
234 |
} |
235 |
} |
236 |
|
237 |
m_dWeights[n][n] = 0.0;
|
238 |
|
239 |
final Matrix subMatrix = m_Matrix.getMatrix(0, n, 0, n); |
240 |
|
241 |
try {
|
242 |
final Matrix inverse = subMatrix.inverse();
|
243 |
m_Matrix.setMatrix(0, n, 0, n, inverse); |
244 |
} |
245 |
catch (final RuntimeException e) { |
246 |
return 0; |
247 |
} |
248 |
|
249 |
//m_dWeights = inverse.getArray();
|
250 |
|
251 |
} |
252 |
|
253 |
return n;
|
254 |
|
255 |
|
256 |
} |
257 |
|
258 |
|
259 |
} |