+ * This is useful so that there's no need to convert the world point into the
+ * grid point by using the mathematical transform, which can lead to rounding
+ * errors which in turn invite all kinds of trouble (see
+ * https://github.com/moovida/jgrasstools/issues/12, for example)
+ *
+ * @author Rafael Almeida (@rafaelalmeida)
+ * @since 0.7.9
+ *
+ */
+public class OmsGridPoint {
+ private GridGeometry2D gridGeometry;
+ private GridCoordinates2D gridCoordinates;
+ private Coordinate worldCoordinates;
+
+ public GridCoordinates2D getGridCoordinates() {
+ return gridCoordinates;
+ }
+
+ public void setGridCoordinates(GridCoordinates2D gridCoordinates) {
+ this.gridCoordinates = gridCoordinates;
+ }
+
+ public Coordinate getWorldCoordinates() {
+ return worldCoordinates;
+ }
+
+ public void setWorldCoordinates(Coordinate worldCoordinates) {
+ this.worldCoordinates = worldCoordinates;
+ }
+
+ public GridGeometry2D getGridGeometry() {
+ return gridGeometry;
+ }
+
+ public void setGridGeometry(GridGeometry2D gridGeometry) {
+ this.gridGeometry = gridGeometry;
+ }
+
+ @Override
+ public String toString() {
+ return "OmsGridPoint [gridCoordinates=" + gridCoordinates
+ + ", worldCoordinates=" + worldCoordinates + "]";
+ }
+
+}
diff --git a/hortonmachine/src/main/java/org/jgrasstools/hortonmachine/modules/statistics/kriging/OmsKriging.java b/hortonmachine/src/main/java/org/jgrasstools/hortonmachine/modules/statistics/kriging/OmsKriging.java
index 8b836c2fd..f97e53373 100644
--- a/hortonmachine/src/main/java/org/jgrasstools/hortonmachine/modules/statistics/kriging/OmsKriging.java
+++ b/hortonmachine/src/main/java/org/jgrasstools/hortonmachine/modules/statistics/kriging/OmsKriging.java
@@ -51,6 +51,7 @@
import java.util.Iterator;
import java.util.LinkedHashMap;
import java.util.List;
+import java.util.Map;
import java.util.Set;
import javax.media.jai.iterator.RandomIterFactory;
@@ -67,13 +68,13 @@
import oms3.annotations.Out;
import oms3.annotations.Status;
+import org.geotools.coverage.grid.GridCoordinates2D;
import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.grid.GridGeometry2D;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.geotools.feature.FeatureCollection;
import org.geotools.feature.FeatureIterator;
import org.geotools.feature.SchemaException;
-import org.geotools.geometry.DirectPosition2D;
import org.jgrasstools.gears.libs.exceptions.ModelsIllegalargumentException;
import org.jgrasstools.gears.libs.exceptions.ModelsRuntimeException;
import org.jgrasstools.gears.libs.modules.JGTModel;
@@ -84,9 +85,7 @@
import org.jgrasstools.gears.utils.math.matrixes.LinearSystem;
import org.jgrasstools.hortonmachine.i18n.HortonMessageHandler;
import org.opengis.feature.simple.SimpleFeature;
-import org.opengis.geometry.DirectPosition;
import org.opengis.geometry.MismatchedDimensionException;
-import org.opengis.referencing.operation.MathTransform;
import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Geometry;
@@ -129,17 +128,19 @@ public class OmsKriging extends JGTModel {
public String fPointZ = null;
/**
- * Define the calculation mode. It can be 0 or 1.
- *
- *
When mode == 1, the values are in a regular grid (the coordinates
- * are stored in a {@link GridCoverage2D), so parameter gridToInterpolate
- * must be set, and the calculated values will be in the outGrid field.
- */
+ * Define the calculation mode. It can be 0 or 1.
+ *
+ *
+ * When mode == 0, the values to calculate are in a non-regular grid (the
+ * coordinates are stored in a {@link FeatureCollection}, so parameters
+ * inInterpolate and fInterpolateid must be set, and the calculated values
+ * will be in the outData field.
+ *
+ *
+ * When mode == 1, the values are in a regular grid (the coordinates are
+ * stored in a {@link GridCoverage2D}, so parameter gridToInterpolate must
+ * be set, and the calculated values will be in the outGrid field.
+ */
@Description(OMSKRIGING_pMode_DESCRIPTION)
@In
public int pMode = 0;
@@ -168,9 +169,21 @@ public class OmsKriging extends JGTModel {
@In
public boolean doLogarithmic = false;
+ /**
+ * The grid of points to interpolate, used when pMode == 1.
+ *
+ * @see #pMode
+ */
@Description(OMSKRIGING_inInterpolationGrid_DESCRIPTION)
@In
public GridGeometry2D inInterpolationGrid = null;
+
+ /**
+ * The internal coverage name GeoTools will put in the generated coverage
+ * when pMode == 1 ({@linkplain #outGrid} field).
+ */
+ @In
+ public String pInterpolatedGridName = "InterpolatedValuesRaster";
public int defaultVariogramMode = 0;
@@ -216,6 +229,12 @@ public class OmsKriging extends JGTModel {
private double west;
private double xres;
private double yres;
+
+ /**
+ * The map of point id -> {@linkplain OmsGridPoint} extracted from the
+ * provided inInterpolationGrid when pMode == 1.
+ */
+ private Map pointIdToGridPoint;
/**
* Executing ordinary kriging.
@@ -357,30 +376,21 @@ public void process() throws Exception {
}
}
}
- LinkedHashMap pointsToInterpolateId2Coordinates = null;
- // vecchio int numPointToInterpolate = getNumPoint(inInterpolate);
- int numPointToInterpolate = 0;
-
- /*
- * if the isLogarithmic is true then execute the model with log value.
- */
- // vecchio double[] result = new double[numPointToInterpolate];
if (pMode == 0) {
- pointsToInterpolateId2Coordinates = getCoordinate(numPointToInterpolate, inInterpolate, fInterpolateid);
+ getCoordinate(inInterpolate, fInterpolateid);
} else if (pMode == 1) {
- pointsToInterpolateId2Coordinates = getCoordinate(inInterpolationGrid);
- numPointToInterpolate = pointsToInterpolateId2Coordinates.size();
+ getCoordinate(inInterpolationGrid);
} else {
throw new ModelsIllegalargumentException("The parameter pMode can only be 0 or 1.", this, pm);
}
- Set pointsToInterpolateIdSet = pointsToInterpolateId2Coordinates.keySet();
+ Set pointsToInterpolateIdSet = pointIdToGridPoint.keySet();
Iterator idIterator = pointsToInterpolateIdSet.iterator();
int j = 0;
- // vecchio int[] idArray = new int[inInterpolate.size()];
- int[] idArray = new int[pointsToInterpolateId2Coordinates.size()];
- double[] result = new double[pointsToInterpolateId2Coordinates.size()];
+
+ int[] idArray = new int[pointIdToGridPoint.size()];
+ double[] result = new double[pointIdToGridPoint.size()];
if (n1 != 0) {
if (doLogarithmic) {
for( int i = 0; i < nStaz; i++ ) {
@@ -408,7 +418,7 @@ public void process() throws Exception {
double sum = 0.;
int id = idIterator.next();
idArray[j] = id;
- Coordinate coordinate = (Coordinate) pointsToInterpolateId2Coordinates.get(id);
+ Coordinate coordinate = (Coordinate) pointIdToGridPoint.get(id).getWorldCoordinates();
xStation[n1] = coordinate.x;
yStation[n1] = coordinate.y;
zStation[n1] = coordinate.z;
@@ -460,7 +470,7 @@ public void process() throws Exception {
if (pMode == 0) {
storeResult(result, idArray);
} else {
- storeResult(result, pointsToInterpolateId2Coordinates);
+ storeResult(result);
}
} else {
pm.errorMessage("No rain for this time step");
@@ -475,7 +485,7 @@ public void process() throws Exception {
if (pMode == 0) {
storeResult(result, idArray);
} else {
- storeResult(result, pointsToInterpolateId2Coordinates);
+ storeResult(result);
}
}
}
@@ -523,55 +533,44 @@ private void verifyInput() {
}
/**
- * Store the result in a HashMap (if the mode is 0 or 1)
- *
- * @param result2
- * the result of the model
- * @param id
- * the associated id of the calculating points.
- * @throws SchemaException
- * @throws SchemaException
- */
- private void storeResult( double[] result2, int[] id ) throws SchemaException {
+ * Store the result in a HashMap (if pMode == 0)
+ *
+ * @param interpolatedValues
+ * the result of the model
+ * @param id
+ * the associated id of the calculated points.
+ */
+ private void storeResult( double[] interpolatedValues, int[] id ) {
outData = new HashMap();
- for( int i = 0; i < result2.length; i++ ) {
- outData.put(id[i], new double[]{checkResultValue(result2[i])});
+ for( int i = 0; i < interpolatedValues.length; i++ ) {
+ outData.put(id[i], new double[]{checkResultValue(interpolatedValues[i])});
}
}
- private void storeResult( double[] interpolatedValues, HashMap interpolatedCoordinatesMap )
+ private void storeResult( double[] interpolatedValues )
throws MismatchedDimensionException, Exception {
WritableRandomIter outIter = RandomIterFactory.createWritable(outWR, null);
- Set pointsToInterpolateIdSett = interpolatedCoordinatesMap.keySet();
- Iterator idIterator = pointsToInterpolateIdSett.iterator();
+ Set pointsToInterpolateIdSet = pointIdToGridPoint.keySet();
+ Iterator idIterator = pointsToInterpolateIdSet.iterator();
int c = 0;
- MathTransform transf = inInterpolationGrid.getCRSToGrid2D();
-
- final DirectPosition gridPoint = new DirectPosition2D();
-
+
while( idIterator.hasNext() ) {
+ // Advance iterator
int id = idIterator.next();
- Coordinate coordinate = (Coordinate) interpolatedCoordinatesMap.get(id);
-
- DirectPosition point = new DirectPosition2D(inInterpolationGrid.getCoordinateReferenceSystem(), coordinate.x,
- coordinate.y);
- transf.transform(point, gridPoint);
-
- double[] gridCoord = gridPoint.getCoordinate();
- int x = (int) gridCoord[0];
- int y = (int) gridCoord[1];
-
- outIter.setSample(x, y, 0, checkResultValue(interpolatedValues[c]));
+
+ // Get grid coordinates and set sample
+ GridCoordinates2D coords = pointIdToGridPoint.get(id).getGridCoordinates();
+ outIter.setSample(coords.x, coords.y, 0, checkResultValue(interpolatedValues[c]));
c++;
-
}
- RegionMap regionMap = CoverageUtilities.gridGeometry2RegionParamsMap(inInterpolationGrid);
-
- outGrid = CoverageUtilities
- .buildCoverage("gridded", outWR, regionMap, inInterpolationGrid.getCoordinateReferenceSystem());
+ // Build the output coverage grid
+ RegionMap regionMap = CoverageUtilities
+ .gridGeometry2RegionParamsMap(inInterpolationGrid);
+ outGrid = CoverageUtilities.buildCoverage(pInterpolatedGridName, outWR,
+ regionMap, inInterpolationGrid.getCoordinateReferenceSystem());
}
@@ -584,6 +583,8 @@ private double checkResultValue( double resultValue ) {
private LinkedHashMap getCoordinate( GridGeometry2D grid ) {
LinkedHashMap out = new LinkedHashMap();
+ pointIdToGridPoint = new LinkedHashMap();
+
int count = 0;
RegionMap regionMap = CoverageUtilities.gridGeometry2RegionParamsMap(grid);
cols = regionMap.getCols();
@@ -600,11 +601,25 @@ private LinkedHashMap getCoordinate( GridGeometry2D grid )
for( int i = 0; i < cols; i++ ) {
easting = easting + xres;
for( int j = 0; j < rows; j++ ) {
+ // Calculate world coordinate
+ Coordinate worldCoordinate = new Coordinate();
northing = northing + yres;
- Coordinate coordinate = new Coordinate();
- coordinate.x = west + i * xres;
- coordinate.y = south + j * yres;
- out.put(count, coordinate);
+ worldCoordinate.x = west + i * xres;
+ worldCoordinate.y = south + j * yres;
+
+ // Pre-calculate grid coordinate and put it on the map
+ GridCoordinates2D gridCoordinates = new GridCoordinates2D(i, rows - j - 1);
+
+ // Create the OmsGridPoint object
+ OmsGridPoint gridPoint = new OmsGridPoint();
+ gridPoint.setGridGeometry(inInterpolationGrid);
+ gridPoint.setWorldCoordinates(worldCoordinate);
+ gridPoint.setGridCoordinates(gridCoordinates);
+
+ // Add the point to the map
+ pointIdToGridPoint.put(count, gridPoint);
+
+ // Advance count
count++;
}
}
@@ -616,14 +631,13 @@ private LinkedHashMap getCoordinate( GridGeometry2D grid )
* Extract the coordinate of a FeatureCollection in a HashMap with an ID as
* a key.
*
- * @param nStaz
* @param collection
* @throws Exception
* if a fiel of elevation isn't the same of the collection
*/
- private LinkedHashMap getCoordinate( int nStaz, SimpleFeatureCollection collection, String idField )
+ private void getCoordinate( SimpleFeatureCollection collection, String idField )
throws Exception {
- LinkedHashMap id2CoordinatesMap = new LinkedHashMap();
+ pointIdToGridPoint = new LinkedHashMap();
FeatureIterator iterator = collection.features();
Coordinate coordinate = null;
try {
@@ -641,13 +655,16 @@ private LinkedHashMap getCoordinate( int nStaz, SimpleFeatu
}
}
coordinate.z = z;
- id2CoordinatesMap.put(name, coordinate);
+
+ // Create the point
+ OmsGridPoint point = new OmsGridPoint();
+ point.setWorldCoordinates(coordinate);
+
+ pointIdToGridPoint.put(name, point);
}
} finally {
iterator.close();
}
-
- return id2CoordinatesMap;
}
/**
diff --git a/hortonmachine/src/test/java/org/jgrasstools/hortonmachine/models/hm/TestKriging.java b/hortonmachine/src/test/java/org/jgrasstools/hortonmachine/models/hm/TestKriging.java
index 21be0882c..5cdc5f60e 100644
--- a/hortonmachine/src/test/java/org/jgrasstools/hortonmachine/models/hm/TestKriging.java
+++ b/hortonmachine/src/test/java/org/jgrasstools/hortonmachine/models/hm/TestKriging.java
@@ -1,677 +1,205 @@
-//package org.jgrasstools.hortonmachine.models.hm;
-//
-//import java.awt.geom.Point2D;
-//import java.io.File;
-//import java.net.URL;
-//import java.util.HashMap;
-//import java.util.Iterator;
-//import java.util.Set;
-//
-//import org.geotools.coverage.grid.GridCoverage2D;
-//import org.geotools.coverage.grid.GridGeometry2D;
-//import org.geotools.data.simple.SimpleFeatureCollection;
-//import org.geotools.filter.text.cql2.CQL;
-//import org.geotools.referencing.CRS;
-//import org.jgrasstools.gears.io.rasterreader.OmsRasterReader;
-//import org.jgrasstools.gears.io.rasterwriter.OmsRasterWriter;
-//import org.jgrasstools.gears.io.shapefile.OmsShapefileFeatureReader;
-//import org.jgrasstools.gears.io.timedependent.OmsTimeSeriesIteratorReader;
-//import org.jgrasstools.gears.io.timedependent.OmsTimeSeriesIteratorWriter;
-//import org.jgrasstools.gears.utils.coverage.CoverageUtilities;
-//import org.jgrasstools.hortonmachine.modules.statistics.kriging.OmsKriging;
-//import org.jgrasstools.hortonmachine.utils.HMTestCase;
-//import org.jgrasstools.hortonmachine.utils.HMTestMaps;
-//import org.opengis.feature.simple.SimpleFeature;
-//import org.opengis.filter.Filter;
-//import org.opengis.referencing.crs.CoordinateReferenceSystem;
-//
-//import com.vividsolutions.jts.geom.Coordinate;
-//import com.vividsolutions.jts.geom.Geometry;
-//
-///**
-// * Test the kriging model.
-// *
-// * @author daniele andreis
-// *
-// */
-//public class TestKriging extends HMTestCase {
-//
-// private File stazioniFile;
-// private File puntiFile;
-// private File krigingRainFile;
-// private String interpolatedRainPath;
-// private File krigingRain2File;
-// private File krigingRain3File;
-// private File krigingRain4File;
-// private File stazioniGridFile;
-//
-// @Override
-// protected void setUp() throws Exception {
-//
-// URL stazioniUrl = this.getClass().getClassLoader().getResource("rainstations.shp");
-// stazioniFile = new File(stazioniUrl.toURI());
-//
-// URL puntiUrl = this.getClass().getClassLoader().getResource("basins_passirio_width0.shp");
-// puntiFile = new File(puntiUrl.toURI());
-//
-// URL krigingRainUrl = this.getClass().getClassLoader().getResource("rain_test.csv");
-// krigingRainFile = new File(krigingRainUrl.toURI());
-//
-// URL krigingRain2Url = this.getClass().getClassLoader().getResource("rain_test2A.csv");
-// krigingRain2File = new File(krigingRain2Url.toURI());
-//
-// URL krigingRain3Url = this.getClass().getClassLoader().getResource("rain_test3A.csv");
-// krigingRain3File = new File(krigingRain3Url.toURI());
-//
-// URL stazioniGridUrl = this.getClass().getClassLoader().getResource("rainstationgrid.shp");
-// stazioniGridFile = new File(stazioniGridUrl.toURI());
-//
-// URL krigingRain4Url = this.getClass().getClassLoader().getResource("rain_test_grid.csv");
-// krigingRain4File = new File(krigingRain4Url.toURI());
-//
-// File interpolatedRainFile = new File(krigingRainFile.getParentFile(), "kriging_interpolated.csv");
-// interpolatedRainPath = interpolatedRainFile.getAbsolutePath();
-// // interpolatedRainPath = interpolatedRainPath.replaceFirst("target",
-// // "src" + File.separator + File.separator + "test");
-// interpolatedRainPath = interpolatedRainPath.replaceFirst("target", "src" + File.separator + "test");
-//
-// interpolatedRainPath = interpolatedRainPath.replaceFirst("test-classes", "resources");
-//
-// super.setUp();
-// }
-//
-// @SuppressWarnings("nls")
-// public void testKriging() throws Exception {
-// //
-// String stationIdField = "ID_PUNTI_M";
-//
-// OmsShapefileFeatureReader stationsReader = new OmsShapefileFeatureReader();
-// stationsReader.file = stazioniGridFile.getAbsolutePath();
-// stationsReader.readFeatureCollection();
-// SimpleFeatureCollection stationsFC = stationsReader.geodata;
-//
-// // OmsShapefileFeatureReader interpolatedPointsReader = new OmsShapefileFeatureReader();
-// // interpolatedPointsReader.file = puntiFile.getAbsolutePath();
-// // interpolatedPointsReader.readFeatureCollection();
-//
-// OmsTimeSeriesIteratorReader reader = new OmsTimeSeriesIteratorReader();
-// reader.file = krigingRain4File.getAbsolutePath();
-// reader.idfield = "ID";
-// reader.tStart = "2000-01-01 00:00";
-// reader.tTimestep = 60;
-// // reader.tEnd = "2000-01-01 00:00";
-// reader.fileNovalue = "-9999";
-// reader.initProcess();
-//
-// OmsKriging kriging = new OmsKriging();
-// kriging.pm = pm;
-//
-// GridGeometry2D gridGeometry2D = CoverageUtilities.gridGeometryFromRegionValues(5204514.51713, 5141634.51713,
-// 686136.82243, 601576.82243, 2114, 1572, HMTestMaps.getCrs());
-// kriging.inInterpolationGrid = gridGeometry2D;
-//
-// kriging.inStations = stationsFC;
-// kriging.fStationsid = stationIdField;
-//
-// // kriging.inInterpolate = interpolatedPointsFC;
-// kriging.fInterpolateid = "netnum";
-//
-// // it doesn't execute the model with log value.
-// kriging.doLogarithmic = false;
-// /*
-// * Set up the model in order to use the variogram with an explicit
-// * integral scale and variance.
-// */
-// // kriging.pVariance = 3.5;
-// // kriging.pIntegralscale = new double[]{10000, 10000, 100};
-// kriging.defaultVariogramMode = 1;
-// kriging.pA = 123537.0;
-// kriging.pNug = 0.0;
-// kriging.pS = 1.678383;
-// /*
-// * Set up the model in order to run with a FeatureCollection as point to
-// * interpolated. In this case only 2D.
-// */
-// kriging.pMode = 1;
-// kriging.pSemivariogramType = 1;
-//
-// // OmsTimeSeriesIteratorWriter writer = new OmsTimeSeriesIteratorWriter();
-// // writer.file = interpolatedRainPath;
-// //
-// // writer.tStart = reader.tStart;
-// // writer.tTimestep = reader.tTimestep;
-// while( reader.doProcess ) {
-// reader.nextRecord();
-// HashMap id2ValueMap = reader.outFolder;
-// kriging.inData = id2ValueMap;
-// kriging.executeKriging();
-// /*
-// * Extract the result.
-// */
-//
-// double[] values = id2ValueMap.get(1331);
-// Filter filter = CQL.toFilter(stationIdField + " = 1331");
-// SimpleFeatureCollection subCollection = stationsFC.subCollection(filter);
-// assertTrue(subCollection.size() == 1);
-//
-// SimpleFeature station = subCollection.features().next();
-// Geometry geometry = (Geometry) station.getDefaultGeometry();
-// Coordinate stationCoordinate = geometry.getCoordinate();
-//
-// GridCoverage2D krigingRaster = kriging.outGrid;
-// double[] expected = krigingRaster.evaluate(new Point2D.Double(stationCoordinate.x, stationCoordinate.y),
-// (double[]) null);
-//
-// assertEquals(expected[0], values[0], 0.01);
-//
-// // HashMap result = kriging.outFolder;
-// // Set pointsToInterpolateResult = result.keySet();
-// // Iterator iteratorTest = pointsToInterpolateResult
-// // .iterator();
-//
-// int iii = 0;
-// // while (iteratorTest.hasNext() && iii<12) {
-// // double expected;
-// // if (j == 0) {
-// // expected = 0.3390869;
-// // } else if (j == 1) {
-// // expected = 0.2556174;
-// // } else if (j == 2) {
-// // expected = 0.2428944;
-// // } else if (j == 3) {
-// // expected = 0.2613782;
-// // } else if (j == 4) {
-// // expected = 0.3112850;
-// // } else if (j == 5) {
-// // expected = 0.2983679;
-// // } else if (j == 6) {
-// // expected = 0.3470377;
-// // } else if (j == 7) {
-// // expected = 0.3874065;
-// // } else if (j == 8) {
-// // expected = 0.2820323;
-// // } else if (j == 9) {
-// // expected = 0.1945515;
-// // } else if (j == 10) {
-// // expected = 0.1698022;
-// // } else if (j == 11) {
-// // expected = 0.2405134;
-// // } else if (j == 12) {
-// // expected = 0.2829313;
-// // } else {
-// // expected = 1.0;
-// // }
-// //
-//
-// //
-// // int id = iteratorTest.next();
-// // double[] actual = result.get(id);
-// // iii+=1;
-// //
-// // //assertEquals(expected, actual[0], 0.001);
-// // j=j+1;
-// // }
-// // iii=0;
-// // j=0;
-// // writer.inData = result;
-// // writer.writeNextLine();
-//
-// }
-//
-// reader.close();
-// // writer.close();
-// }
-// //
-// //
-//
-// // ///////////////////////////////////////////////////////////////////////////////////////////
-// // /////////////////////////////////TEST 1
-// // PASSA////////////////////////////////////////////////////
-// // /////////////////////////////////////////////////////////////////////////////////////////
-// public void testKriging1() throws Exception {
-//
-// OmsShapefileFeatureReader stationsReader = new OmsShapefileFeatureReader();
-// stationsReader.file = stazioniFile.getAbsolutePath();
-// stationsReader.readFeatureCollection();
-// SimpleFeatureCollection stationsFC = stationsReader.geodata;
-//
-// OmsShapefileFeatureReader interpolatedPointsReader = new OmsShapefileFeatureReader();
-// interpolatedPointsReader.file = puntiFile.getAbsolutePath();
-// interpolatedPointsReader.readFeatureCollection();
-// SimpleFeatureCollection interpolatedPointsFC = interpolatedPointsReader.geodata;
-//
-// OmsTimeSeriesIteratorReader reader = new OmsTimeSeriesIteratorReader();
-// reader.file = krigingRainFile.getAbsolutePath();
-// reader.idfield = "ID";
-// reader.tStart = "2000-01-01 00:00";
-// reader.tTimestep = 60;
-// // reader.tEnd = "2000-01-01 00:00";
-// reader.fileNovalue = "-9999";
-//
-// reader.initProcess();
-//
-// OmsKriging kriging = new OmsKriging();
-// kriging.pm = pm;
-//
-// kriging.inStations = stationsFC;
-// kriging.fStationsid = "ID_PUNTI_M";
-//
-// kriging.inInterpolate = interpolatedPointsFC;
-// kriging.fInterpolateid = "netnum";
-//
-// // it doesn't execute the model with log value.
-// kriging.doLogarithmic = false;
-// /*
-// * Set up the model in order to use the variogram with an explicit
-// * integral scale and variance.
-// */
-// // kriging.pVariance = 3.5;
-// // kriging.pIntegralscale = new double[]{10000, 10000, 100};
-// kriging.defaultVariogramMode = 1;
-// kriging.pA = 123537.0;
-// kriging.pNug = 0.0;
-// kriging.pS = 1.678383;
-// /*
-// * Set up the model in order to run with a FeatureCollection as point to
-// * interpolated. In this case only 2D.
-// */
-// kriging.pMode = 0;
-// kriging.pSemivariogramType = 1;
-//
-// OmsTimeSeriesIteratorWriter writer = new OmsTimeSeriesIteratorWriter();
-// writer.file = interpolatedRainPath;
-//
-// writer.tStart = reader.tStart;
-// writer.tTimestep = reader.tTimestep;
-// int j = 0;
-// while( reader.doProcess ) {
-// reader.nextRecord();
-// HashMap id2ValueMap = reader.outFolder;
-// kriging.inData = id2ValueMap;
-// kriging.executeKriging();
-// /*
-// * Extract the result.
-// */
-// HashMap result = kriging.outFolder;
-// Set pointsToInterpolateResult = result.keySet();
-// Iterator iteratorTest = pointsToInterpolateResult.iterator();
-//
-// int iii = 0;
-//
-// while( iteratorTest.hasNext() && iii < 12 ) {
-// double expected;
-// if (j == 0) {
-// expected = 0.3390869;
-// } else if (j == 1) {
-// expected = 0.2556174;
-// } else if (j == 2) {
-// expected = 0.2428944;
-// } else if (j == 3) {
-// expected = 0.2613782;
-// } else if (j == 4) {
-// expected = 0.3112850;
-// } else if (j == 5) {
-// expected = 0.2983679;
-// } else if (j == 6) {
-// expected = 0.3470377;
-// } else if (j == 7) {
-// expected = 0.3874065;
-// } else if (j == 8) {
-// expected = 0.2820323;
-// } else if (j == 9) {
-// expected = 0.1945515;
-// } else if (j == 10) {
-// expected = 0.1698022;
-// } else if (j == 11) {
-// expected = 0.2405134;
-// } else if (j == 12) {
-// expected = 0.2829313;
-// } else {
-// expected = 1.0;
-// }
-//
-// int id = iteratorTest.next();
-// double[] actual = result.get(id);
-// iii += 1;
-//
-// assertEquals(expected, actual[0], 0.001);
-// j = j + 1;
-// }
-// iii = 0;
-// j = 0;
-// writer.inData = result;
-// writer.writeNextLine();
-//
-// }
-//
-// reader.close();
-// writer.close();
-// }
-//
-// // ///////////////////////////////////////////////////////////////////////////////////////////
-// // /////////////////////////////////FINE TEST
-// // 1PASSA////////////////////////////////////////////////////
-// // /////////////////////////////////////////////////////////////////////////////////////////
-// //
-// //
-// //
-// // ///////////////////////////////////////////////////////////////////////////////////////////
-// // /////////////////////////////////TEST 2
-// // PASSA////////////////////////////////////////////////////
-// // /////////////////////////////////////////////////////////////////////////////////////////
-//
-// /**
-// * Run the kriging models.
-// *
-// *
-// * This is the case which all the station have the same value.
-// *
-// * @throws Exception
-// * @throws Exception
-// */
-// public void testKriging2() throws Exception {
-// OmsShapefileFeatureReader stationsReader = new OmsShapefileFeatureReader();
-// stationsReader.file = stazioniFile.getAbsolutePath();
-// stationsReader.readFeatureCollection();
-// SimpleFeatureCollection stationsFC = stationsReader.geodata;
-//
-// OmsShapefileFeatureReader interpolatedPointsReader = new OmsShapefileFeatureReader();
-// interpolatedPointsReader.file = puntiFile.getAbsolutePath();
-// interpolatedPointsReader.readFeatureCollection();
-// SimpleFeatureCollection interpolatedPointsFC = interpolatedPointsReader.geodata;
-//
-// OmsTimeSeriesIteratorReader reader = new OmsTimeSeriesIteratorReader();
-// reader.file = krigingRain2File.getAbsolutePath();
-// reader.idfield = "ID";
-// reader.tStart = "2000-01-01 00:00";
-// reader.tTimestep = 60;
-// // reader.tEnd = "2000-01-01 00:00";
-// reader.fileNovalue = "-9999";
-//
-// reader.initProcess();
-//
-// OmsKriging kriging = new OmsKriging();
-// kriging.pm = pm;
-//
-// kriging.inStations = stationsFC;
-// kriging.fStationsid = "ID_PUNTI_M";
-//
-// kriging.inInterpolate = interpolatedPointsFC;
-// kriging.fInterpolateid = "netnum";
-//
-// // it doesn't execute the model with log value.
-// kriging.doLogarithmic = false;
-// /*
-// * Set up the model in order to use the variogram with an explicit integral scale and
-// variance.
-// */
-// kriging.pVariance = 0.5;
-// kriging.pIntegralscale = new double[]{10000, 10000, 100};
-// /*
-// * Set up the model in order to run with a FeatureCollection as point to interpolated. In this
-// case only 2D.
-// */
-// kriging.pMode = 0;
-//
-// OmsTimeSeriesIteratorWriter writer = new OmsTimeSeriesIteratorWriter();
-// writer.file = interpolatedRainPath;
-//
-// writer.tStart = reader.tStart;
-// writer.tTimestep = reader.tTimestep;
-//
-// while( reader.doProcess ) {
-// reader.nextRecord();
-// HashMap id2ValueMap = reader.outFolder;
-// kriging.inData = id2ValueMap;
-// kriging.executeKriging();
-// /*
-// * Extract the result.
-// */
-// HashMap result = kriging.outFolder;
-// Set pointsToInterpolateResult = result.keySet();
-// Iterator iterator = pointsToInterpolateResult.iterator();
-// while( iterator.hasNext() ) {
-// int id = iterator.next();
-// double[] actual = result.get(id);
-// assertEquals(1.0, actual[0], 0);
-// }
-// writer.inData = result;
-// writer.writeNextLine();
-// }
-//
-// reader.close();
-// writer.close();
-// }
-// // /////////////////////////////////////////////////////////////////////////////////////////
-// // ///////////////////////////////FINE TEST 2
-// // PASSA////////////////////////////////////////////////////
-// // ///////////////////////////////////////////////////////////////////////////////////////
-//
-// // /////////////////////////////////////////////////////////////////////////////////////////
-// // /////////////////////////////// TEST 3
-// // PASSA////////////////////////////////////////////////////
-// // ///////////////////////////////////////////////////////////////////////////////////////
-// // /**
-// // * Run the kriging models.
-// // *
-// // *
-// // * This is the case that defaultMode=0.
-// // *
-// // * @throws Exception
-// // * @throws Exception
-// // */
-// public void testKriging4() throws Exception {
-// OmsShapefileFeatureReader stationsReader = new OmsShapefileFeatureReader();
-// stationsReader.file = stazioniFile.getAbsolutePath();
-// stationsReader.readFeatureCollection();
-// SimpleFeatureCollection stationsFC = stationsReader.geodata;
-//
-// OmsShapefileFeatureReader interpolatedPointsReader = new OmsShapefileFeatureReader();
-// interpolatedPointsReader.file = puntiFile.getAbsolutePath();
-// interpolatedPointsReader.readFeatureCollection();
-// SimpleFeatureCollection interpolatedPointsFC = interpolatedPointsReader.geodata;
-//
-// OmsTimeSeriesIteratorReader reader = new OmsTimeSeriesIteratorReader();
-// reader.file = krigingRainFile.getAbsolutePath();
-// reader.idfield = "ID";
-// reader.tStart = "2000-01-01 00:00";
-// reader.tTimestep = 60;
-// // reader.tEnd = "2000-01-01 00:00";
-// reader.fileNovalue = "-9999";
-//
-// reader.initProcess();
-//
-// OmsKriging kriging = new OmsKriging();
-// kriging.pm = pm;
-//
-// kriging.inStations = stationsFC;
-// kriging.fStationsid = "ID_PUNTI_M";
-//
-// kriging.inInterpolate = interpolatedPointsFC;
-// kriging.fInterpolateid = "netnum";
-//
-// // it doesn't execute the model with log value.
-// kriging.doLogarithmic = false;
-// /*
-// * Set up the model in order to use the variogram with an explicit integral scale and
-// variance.
-// */
-// kriging.pVariance = 3.5;
-// kriging.pIntegralscale = new double[]{10000, 10000, 100};
-// /*
-// * Set up the model in order to run with a FeatureCollection as point to interpolated. In this
-// case only 2D.
-// */
-// kriging.pMode = 0;
-//
-// kriging.doIncludezero = false;
-// OmsTimeSeriesIteratorWriter writer = new OmsTimeSeriesIteratorWriter();
-// writer.file = interpolatedRainPath;
-//
-// writer.tStart = reader.tStart;
-// writer.tTimestep = reader.tTimestep;
-//
-// while( reader.doProcess ) {
-// reader.nextRecord();
-// HashMap id2ValueMap = reader.outFolder;
-// kriging.inData = id2ValueMap;
-// kriging.executeKriging();
-// /*
-// * Extract the result.
-// */
-// HashMap result = kriging.outFolder;
-// double[][] test = HMTestMaps.outKriging4;
-// for( int i = 0; i < test.length; i++ ) {
-// double actual = result.get((int) test[i][0])[0];
-// double expected = test[i][1];
-// assertEquals(expected, actual, 0.01);
-// }
-//
-// writer.inData = result;
-// writer.writeNextLine();
-// }
-//
-// reader.close();
-// writer.close();
-// }
-// // /////////////////////////////////////////////////////////////////////////////////////////
-// // ///////////////////////////////FINE TEST 3
-// // PASSA////////////////////////////////////////////////////
-// // ///////////////////////////////////////////////////////////////////////////////////////
-//
-// // /////////////////////////////////////////////////////////////////////////////////////////
-// // /////////////////////////////// TEST 4
-// // PASSA////////////////////////////////////////////////////
-// // ///////////////////////////////////////////////////////////////////////////////////////
-// /**
-// * Run the kriging models.
-// *
-// *
-// * This is the case which there is only one station.
-// *
-// * @throws Exception
-// * @throws Exception
-// */
-// public void testKriging5() throws Exception {
-// OmsShapefileFeatureReader stationsReader = new OmsShapefileFeatureReader();
-// stationsReader.file = stazioniFile.getAbsolutePath();
-// stationsReader.readFeatureCollection();
-// SimpleFeatureCollection stationsFC = stationsReader.geodata;
-//
-// OmsShapefileFeatureReader interpolatedPointsReader = new OmsShapefileFeatureReader();
-// interpolatedPointsReader.file = puntiFile.getAbsolutePath();
-// interpolatedPointsReader.readFeatureCollection();
-// SimpleFeatureCollection interpolatedPointsFC = interpolatedPointsReader.geodata;
-//
-// OmsTimeSeriesIteratorReader reader = new OmsTimeSeriesIteratorReader();
-// reader.file = krigingRain3File.getAbsolutePath();
-// reader.idfield = "ID";
-// reader.tStart = "2000-01-01 00:00";
-// reader.tTimestep = 60;
-// // reader.tEnd = "2000-01-01 00:00";
-// reader.fileNovalue = "-9999";
-//
-// reader.initProcess();
-//
-// OmsKriging kriging = new OmsKriging();
-// kriging.pm = pm;
-//
-// kriging.inStations = stationsFC;
-// kriging.fStationsid = "ID_PUNTI_M";
-//
-// kriging.inInterpolate = interpolatedPointsFC;
-// kriging.fInterpolateid = "netnum";
-//
-// // it doesn't execute the model with log value.
-// kriging.doLogarithmic = false;
-// /*
-// * Set up the model in order to use the variogram with an explicit integral scale and
-// variance.
-// */
-// kriging.pVariance = 0.5;
-// kriging.pIntegralscale = new double[]{10000, 10000, 100};
-// /*
-// * Set up the model in order to run with a FeatureCollection as point to interpolated. In this
-// case only 2D.
-// */
-// kriging.pMode = 0;
-//
-// OmsTimeSeriesIteratorWriter writer = new OmsTimeSeriesIteratorWriter();
-// writer.file = interpolatedRainPath;
-//
-// writer.tStart = reader.tStart;
-// writer.tTimestep = reader.tTimestep;
-// int j = 0;
-// while( reader.doProcess ) {
-// reader.nextRecord();
-// HashMap id2ValueMap = reader.outFolder;
-// kriging.inData = id2ValueMap;
-// kriging.executeKriging();
-// /*
-// * Extract the result.
-// */
-// HashMap result = kriging.outFolder;
-// Set pointsToInterpolateResult = result.keySet();
-// Iterator iteratorTest = pointsToInterpolateResult.iterator();
-// double expected;
-// if (j == 0) {
-// expected = 10.0;
-// } else if (j == 1) {
-// expected = 15;
-// } else if (j == 2) {
-// expected = 1;
-// } else if (j == 3) {
-// expected = 2;
-// } else if (j == 4) {
-// expected = 2;
-// } else if (j == 5) {
-// expected = 0;
-// } else if (j == 6) {
-// expected = 0;
-// } else if (j == 7) {
-// expected = 23;
-// } else if (j == 8) {
-// expected = 50;
-// } else if (j == 9) {
-// expected = 70;
-// } else if (j == 10) {
-// expected = 30;
-// } else if (j == 11) {
-// expected = 10;
-// } else if (j == 12) {
-// expected = 2;
-// } else {
-// expected = 1.0;
-// }
-//
-// while( iteratorTest.hasNext() ) {
-// int id = iteratorTest.next();
-// double[] actual = result.get(id);
-//
-// assertEquals(expected, actual[0], 0);
-// }
-// writer.inData = result;
-// writer.writeNextLine();
-// j++;
-// }
-//
-// reader.close();
-// writer.close();
-// }
-// // ///////////////////////////////////////////////////////////////////////////////////////
-// // /////////////////////////////FINE TEST 4
-// // PASSA////////////////////////////////////////////////////
-// // /////////////////////////////////////////////////////////////////////////////////////
-// @Override
-// protected void tearDown() throws Exception {
-// File remove = new File(interpolatedRainPath);
-// if (remove.exists()) {
-// if (!remove.delete()) {
-// remove.deleteOnExit();
-// }
-// }
-//
-// super.tearDown();
-// }
-//
-//}
+package org.jgrasstools.hortonmachine.models.hm;
+
+import static org.junit.Assert.assertEquals;
+import static org.junit.Assert.assertTrue;
+
+import java.awt.geom.Point2D;
+import java.awt.image.ComponentSampleModel;
+import java.awt.image.DataBuffer;
+import java.awt.image.SampleModel;
+import java.awt.image.WritableRaster;
+import java.io.ByteArrayOutputStream;
+import java.io.File;
+import java.io.FileOutputStream;
+import java.io.IOException;
+import java.net.URL;
+import java.util.HashMap;
+
+import javax.media.jai.RasterFactory;
+
+import org.apache.commons.io.IOUtils;
+import org.geotools.coverage.CoverageFactoryFinder;
+import org.geotools.coverage.grid.GridCoordinates2D;
+import org.geotools.coverage.grid.GridCoverage2D;
+import org.geotools.coverage.grid.GridCoverageFactory;
+import org.geotools.coverage.grid.GridGeometry2D;
+import org.geotools.data.simple.SimpleFeatureCollection;
+import org.geotools.filter.text.cql2.CQL;
+import org.geotools.gce.geotiff.GeoTiffWriter;
+import org.jgrasstools.gears.io.shapefile.OmsShapefileFeatureReader;
+import org.jgrasstools.gears.io.timedependent.OmsTimeSeriesIteratorReader;
+import org.jgrasstools.gears.libs.monitor.DummyProgressMonitor;
+import org.jgrasstools.gears.utils.coverage.CoverageUtilities;
+import org.jgrasstools.hortonmachine.modules.statistics.kriging.OmsKriging;
+import org.jgrasstools.hortonmachine.utils.HMTestMaps;
+import org.junit.Test;
+import org.opengis.feature.simple.SimpleFeature;
+import org.opengis.filter.Filter;
+
+import com.vividsolutions.jts.geom.Coordinate;
+import com.vividsolutions.jts.geom.Geometry;
+
+/**
+ * Test the kriging model.
+ *
+ * @author Daniele Andreis
+ * @author Rafael Almeida
+ *
+ */
+public class TestKriging {
+
+ /**
+ * If enabled, will write test results into GeoTIFF files for manual
+ * inspection and ease of development, beyond automated testing.
+ */
+ private final static boolean ENABLE_TEST_GEOTIFF_WRITING = true;
+
+ @Test
+ public void testKriging() throws Exception {
+ // Load test case shapefile
+ URL stazioniGridUrl = this.getClass().getClassLoader().getResource("rainstationgrid.shp");
+ File stazioniGridFile = new File(stazioniGridUrl.toURI());
+
+ // Read the features
+ String stationIdField = "ID_PUNTI_M";
+ OmsShapefileFeatureReader stationsReader = new OmsShapefileFeatureReader();
+ stationsReader.file = stazioniGridFile.getAbsolutePath();
+ stationsReader.readFeatureCollection();
+ SimpleFeatureCollection stationsFC = stationsReader.geodata;
+
+ // Load the time series data for this test case
+ URL krigingRain4Url = this.getClass().getClassLoader().getResource("rain_test_grid.csv");
+ File krigingRain4File = new File(krigingRain4Url.toURI());
+
+ // Setup the time series reader
+ OmsTimeSeriesIteratorReader reader = new OmsTimeSeriesIteratorReader();
+ reader.file = krigingRain4File.getAbsolutePath();
+ reader.idfield = "ID";
+ reader.tStart = "2000-01-01 00:00";
+ reader.tTimestep = 60;
+ reader.fileNovalue = "-9999";
+ reader.initProcess();
+
+ // Create the kriging handler
+ OmsKriging kriging = new OmsKriging();
+ kriging.pm = new DummyProgressMonitor();
+
+ // Defines the grids of points to be interpolated
+ kriging.pMode = 1;
+ GridGeometry2D gridGeometry2D = CoverageUtilities.gridGeometryFromRegionValues(5204514.51713, 5141634.51713,
+ 686136.82243, 601576.82243, 2114, 1572, HMTestMaps.getCrs());
+ kriging.inInterpolationGrid = gridGeometry2D;
+
+ // Set up the station data
+ kriging.inStations = stationsFC;
+ kriging.fStationsid = stationIdField;
+ kriging.fInterpolateid = "netnum";
+
+ // Disable logarithmic mode (which would do the kriging with the log of
+ // data)
+ kriging.doLogarithmic = false;
+
+ // Set up the exponential variogram
+ kriging.defaultVariogramMode = 1;
+ kriging.pSemivariogramType = 1;
+ kriging.pA = 123537.0;
+ kriging.pNug = 0.0;
+ kriging.pS = 1.678383;
+
+ // Run the kriging for the time steps
+ while( reader.doProcess ) {
+ // Run this time step
+ reader.nextRecord();
+ HashMap id2ValueMap = reader.outData;
+ kriging.inData = id2ValueMap;
+ kriging.process();
+
+ // Get expected result
+ double[] values = id2ValueMap.get(1331);
+ Filter filter = CQL.toFilter(stationIdField + " = 1331");
+ SimpleFeatureCollection subCollection = stationsFC.subCollection(filter);
+ assertTrue(subCollection.size() == 1);
+
+ // Get result station coordinates
+ SimpleFeature station = subCollection.features().next();
+ Geometry geometry = (Geometry) station.getDefaultGeometry();
+ Coordinate stationCoordinate = geometry.getCoordinate();
+
+ // Extract and check calculated result
+ GridCoverage2D krigingRaster = kriging.outGrid;
+ double[] expected = krigingRaster.evaluate(new Point2D.Double(stationCoordinate.x, stationCoordinate.y),
+ (double[]) null);
+ assertEquals(expected[0], values[0], 0.01);
+
+ // Save grid for debug purposes
+ if (ENABLE_TEST_GEOTIFF_WRITING) {
+ byte[] geoTIFF = makeGeoTIFF(krigingRaster);
+ FileOutputStream stream = new FileOutputStream(new File("kriging-test-jgrasstools.tif"));
+ IOUtils.write(geoTIFF, stream);
+ }
+ }
+
+ // Close the time series reader
+ reader.close();
+ }
+
+ /**
+ * Creates a GeoTIFF file representing the input coverage and returns
+ * it as a byte array. It can then be persisted to disk (usually with
+ * .tif extension).
+ *
+ * The data in the image is guaranteed to be in 32 bits.
+ * @throws StriderException
+ *
+ * @author Rafael Almeida (@rafaelalmeida)
+ */
+ private byte[] makeGeoTIFF(GridCoverage2D coverage) throws IOException {
+ // Ensure the output will be a 32-bit image
+ coverage = transcodeCoverageTo32Bit(coverage);
+
+ // Convert the coverage to the specified format
+ GridCoverageFactory gridFactory = CoverageFactoryFinder.getGridCoverageFactory(null);
+ gridFactory.create(coverage.getName(), new float[1][1], coverage.getEnvelope());
+
+ // Convert the result into GeoTIFF format
+ ByteArrayOutputStream outputStream = new ByteArrayOutputStream();
+ GeoTiffWriter geoTIFFWriter = new GeoTiffWriter(outputStream);
+ geoTIFFWriter.getFormat();
+ geoTIFFWriter.write(coverage, null);
+
+ // Return the GeoTIFF as a byte array
+ return outputStream.toByteArray();
+ }
+
+ /**
+ * Creates a GridCoverage2D equivalent to the input, but with raster data
+ * in 32-bit IEEE 754 floats.
+ *
+ * @author Rafael Almeida (@rafaelalmeida)
+ */
+ private GridCoverage2D transcodeCoverageTo32Bit(GridCoverage2D coverage) {
+ // Aliases
+ int width = coverage.getRenderedImage().getWidth();
+ int height = coverage.getRenderedImage().getHeight();
+
+ // Prepare the raster
+ SampleModel model = new ComponentSampleModel(DataBuffer.TYPE_FLOAT,
+ width, height, 1, width, new int[] { 0 });
+ WritableRaster raster = RasterFactory.createWritableRaster(model, null);
+
+ // Populate the raster by transcoding the original coverage
+ for (int y = 0; y < height; y++) {
+ for (int x = 0; x < width; x++) {
+ // JLS, section 4.2.3, guarantees float primitives will always be
+ // 32-bit IEEE 754 floats.
+ raster.setPixel(x, y, coverage.evaluate(new GridCoordinates2D(
+ x, y), new float[1]));
+ }
+ }
+
+ // Create and return equivalent coverage with transcoded raster
+ return CoverageFactoryFinder.getGridCoverageFactory(null).create(
+ coverage.getName(), raster, coverage.getEnvelope());
+ }
+
+}
diff --git a/jgrassgears/pom.xml b/jgrassgears/pom.xml
index e9936060d..7621ea657 100644
--- a/jgrassgears/pom.xml
+++ b/jgrassgears/pom.xml
@@ -207,5 +207,10 @@
3.2
+
+
+ 1.7
+ 1.7
+