From 34c340b6ab71720ca91e078534bff002dfd83512 Mon Sep 17 00:00:00 2001 From: Rafael Almeida Date: Wed, 16 Sep 2015 10:18:37 -0300 Subject: [PATCH 1/3] Fixing compilation problem on utils class --- .../jgrasstools/gears/utils/math/NumericsUtilities.java | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/jgrassgears/src/main/java/org/jgrasstools/gears/utils/math/NumericsUtilities.java b/jgrassgears/src/main/java/org/jgrasstools/gears/utils/math/NumericsUtilities.java index 2776d8356..d87b08fce 100644 --- a/jgrassgears/src/main/java/org/jgrasstools/gears/utils/math/NumericsUtilities.java +++ b/jgrassgears/src/main/java/org/jgrasstools/gears/utils/math/NumericsUtilities.java @@ -121,7 +121,7 @@ public static boolean dEq( double a, double b ) { return true; } double diffAbs = abs(a - b); - return a == b ? true : diffAbs < D_TOLERANCE ? true : diffAbs / max(abs(a), abs(b)) < D_TOLERANCE; + return a == b ? true : diffAbs < D_TOLERANCE ? true : diffAbs / Double.max(abs(a), abs(b)) < D_TOLERANCE; } /** @@ -138,7 +138,7 @@ public static boolean dEq( double a, double b, double epsilon ) { return true; } double diffAbs = abs(a - b); - return a == b ? true : diffAbs < epsilon ? true : diffAbs / max(abs(a), abs(b)) < epsilon; + return a == b ? true : diffAbs < epsilon ? true : diffAbs / Double.max(abs(a), abs(b)) < epsilon; } /** @@ -155,7 +155,7 @@ public static boolean fEq( float a, float b ) { return true; } float diffAbs = abs(a - b); - return a == b ? true : diffAbs < F_TOLERANCE ? true : diffAbs / max(abs(a), abs(b)) < F_TOLERANCE; + return a == b ? true : diffAbs < F_TOLERANCE ? true : diffAbs / Float.max(abs(a), abs(b)) < F_TOLERANCE; } /** @@ -172,7 +172,7 @@ public static boolean fEq( float a, float b, float epsilon ) { return true; } float diffAbs = abs(a - b); - return a == b ? true : diffAbs < epsilon ? true : diffAbs / max(abs(a), abs(b)) < epsilon; + return a == b ? true : diffAbs < epsilon ? true : diffAbs / Float.max(abs(a), abs(b)) < epsilon; } /** From 2eab225f51a4229a0592f8f8d4432cbb9b37da3d Mon Sep 17 00:00:00 2001 From: Rafael Almeida Date: Mon, 21 Sep 2015 17:48:20 -0300 Subject: [PATCH 2/3] Fixing blank column problem in kriging (issue 12) - Changed strategy to pre-calculate grid positions, avoiding rounding errors; - Revived one of the unit tests for the kriging module; - Some documentation fixes and improvements; - Minor refactorings for cleaner code. --- hortonmachine/pom.xml | 3 +- .../hortonmachine/i18n/HortonMessages.java | 4 +- .../statistics/kriging/OmsGridPoint.java | 57 ++ .../statistics/kriging/OmsKriging.java | 171 ++-- .../hortonmachine/models/hm/TestKriging.java | 882 ++++-------------- 5 files changed, 359 insertions(+), 758 deletions(-) create mode 100644 hortonmachine/src/main/java/org/jgrasstools/hortonmachine/modules/statistics/kriging/OmsGridPoint.java diff --git a/hortonmachine/pom.xml b/hortonmachine/pom.xml index f4c9b01cb..6111cc63f 100644 --- a/hortonmachine/pom.xml +++ b/hortonmachine/pom.xml @@ -7,9 +7,8 @@ 0.7.9-SNAPSHOT - org.jgrasstools jgt-hortonmachine - 0.7.9-SNAPSHOT + 0.7.9-SNAPSHOT jar The Horton Machine diff --git a/hortonmachine/src/main/java/org/jgrasstools/hortonmachine/i18n/HortonMessages.java b/hortonmachine/src/main/java/org/jgrasstools/hortonmachine/i18n/HortonMessages.java index e4916497b..291d442bb 100644 --- a/hortonmachine/src/main/java/org/jgrasstools/hortonmachine/i18n/HortonMessages.java +++ b/hortonmachine/src/main/java/org/jgrasstools/hortonmachine/i18n/HortonMessages.java @@ -126,8 +126,8 @@ public class HortonMessages { public static final String OMSKRIGING_pA_DESCRIPTION = "The range if the models runs with the gaussian variogram."; public static final String OMSKRIGING_pS_DESCRIPTION = "The sill if the models runs with the gaussian variogram."; public static final String OMSKRIGING_pNug_DESCRIPTION = "Is the nugget if the models runs with the gaussian variogram."; - public static final String OMSKRIGING_outGrid_DESCRIPTION = "The interpolated gridded data (for mode 2 and 3."; - public static final String OMSKRIGING_outData_DESCRIPTION = "The interpolated data (for mode 0 and 1)."; + public static final String OMSKRIGING_outGrid_DESCRIPTION = "The interpolated gridded data (for pMode == 1)."; + public static final String OMSKRIGING_outData_DESCRIPTION = "The interpolated data (for pMode == 0)."; public static final String OMSNETNUMBERING_DESCRIPTION = "Assigns the numbers to the network's links."; public static final String OMSNETNUMBERING_DOCUMENTATION = "OmsNetNumbering.html"; diff --git a/hortonmachine/src/main/java/org/jgrasstools/hortonmachine/modules/statistics/kriging/OmsGridPoint.java b/hortonmachine/src/main/java/org/jgrasstools/hortonmachine/modules/statistics/kriging/OmsGridPoint.java new file mode 100644 index 000000000..32bb92706 --- /dev/null +++ b/hortonmachine/src/main/java/org/jgrasstools/hortonmachine/modules/statistics/kriging/OmsGridPoint.java @@ -0,0 +1,57 @@ +package org.jgrasstools.hortonmachine.modules.statistics.kriging; + +import org.geotools.coverage.grid.GridCoordinates2D; +import org.geotools.coverage.grid.GridGeometry2D; + +import com.vividsolutions.jts.geom.Coordinate; + +/** + * A world coordinate associated with its respective grid point in the context + * of a certain grid geometry. + * + *

+ * 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 == 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. - */ + * 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()); + } + +} From 7f87a855cb1d93e0df2f146f97f481ba7cc4c065 Mon Sep 17 00:00:00 2001 From: Rafael Almeida Date: Tue, 22 Sep 2015 12:18:41 -0300 Subject: [PATCH 3/3] Making compilation work both on Travis and on Eclipse --- hortonmachine/pom.xml | 5 +++++ jgrassgears/pom.xml | 5 +++++ .../jgrasstools/gears/utils/math/NumericsUtilities.java | 8 ++++---- 3 files changed, 14 insertions(+), 4 deletions(-) diff --git a/hortonmachine/pom.xml b/hortonmachine/pom.xml index 6111cc63f..2cb9fc26c 100644 --- a/hortonmachine/pom.xml +++ b/hortonmachine/pom.xml @@ -87,4 +87,9 @@ --> + + + 1.7 + 1.7 + 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 + diff --git a/jgrassgears/src/main/java/org/jgrasstools/gears/utils/math/NumericsUtilities.java b/jgrassgears/src/main/java/org/jgrasstools/gears/utils/math/NumericsUtilities.java index d87b08fce..05f6f2132 100644 --- a/jgrassgears/src/main/java/org/jgrasstools/gears/utils/math/NumericsUtilities.java +++ b/jgrassgears/src/main/java/org/jgrasstools/gears/utils/math/NumericsUtilities.java @@ -121,7 +121,7 @@ public static boolean dEq( double a, double b ) { return true; } double diffAbs = abs(a - b); - return a == b ? true : diffAbs < D_TOLERANCE ? true : diffAbs / Double.max(abs(a), abs(b)) < D_TOLERANCE; + return a == b ? true : diffAbs < D_TOLERANCE ? true : diffAbs / Math.max(abs(a), abs(b)) < D_TOLERANCE; } /** @@ -138,7 +138,7 @@ public static boolean dEq( double a, double b, double epsilon ) { return true; } double diffAbs = abs(a - b); - return a == b ? true : diffAbs < epsilon ? true : diffAbs / Double.max(abs(a), abs(b)) < epsilon; + return a == b ? true : diffAbs < epsilon ? true : diffAbs / Math.max(abs(a), abs(b)) < epsilon; } /** @@ -155,7 +155,7 @@ public static boolean fEq( float a, float b ) { return true; } float diffAbs = abs(a - b); - return a == b ? true : diffAbs < F_TOLERANCE ? true : diffAbs / Float.max(abs(a), abs(b)) < F_TOLERANCE; + return a == b ? true : diffAbs < F_TOLERANCE ? true : diffAbs / Math.max(abs(a), abs(b)) < F_TOLERANCE; } /** @@ -172,7 +172,7 @@ public static boolean fEq( float a, float b, float epsilon ) { return true; } float diffAbs = abs(a - b); - return a == b ? true : diffAbs < epsilon ? true : diffAbs / Float.max(abs(a), abs(b)) < epsilon; + return a == b ? true : diffAbs < epsilon ? true : diffAbs / Math.max(abs(a), abs(b)) < epsilon; } /**