/****************************************************************************** * $Id: GDALGrid.java $ * * Project: GDAL Java applications * Purpose: GDAL scattered data gridding (interpolation) tool * Author: Ivan Lucena, ivan.lucena@pmldnet.com, * translated from gdal_grid.cpp * originally written by Andrey Kiselev, dron@ak4719.spb.edu ****************************************************************************** * Copyright (c) 2010, Ivan Lucena * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal * in the Software without restriction, including without limitation the rights * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell * copies of the Software, and to permit persons to whom the Software is * furnished to do so, subject to the following conditions: * * The above copyright notice and this permission notice shall be included in * all copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE * SOFTWARE. ****************************************************************************/ import java.util.ArrayList; import java.util.List; import java.nio.ByteBuffer; import org.gdal.gdal.Band; import org.gdal.gdal.Dataset; import org.gdal.gdal.Driver; import org.gdal.gdal.ProgressCallback; import org.gdal.gdal.TermProgressCallback; import org.gdal.gdal.gdal; import org.gdal.gdalconst.gdalconstConstants; import org.gdal.ogr.DataSource; import org.gdal.ogr.Feature; import org.gdal.ogr.Geometry; import org.gdal.ogr.Layer; import org.gdal.ogr.ogr; import org.gdal.ogr.ogrConstants; import org.gdal.osr.SpatialReference; public class GDALGrid { public static void Usage() { System.out .println("Usage: gridcreate [--help-general] [--formats]\n" + " [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/\n" + " CInt16/CInt32/CFloat32/CFloat64}]\n" + " [-of format] [-co \"NAME=VALUE\"]\n" + " [-zfield field_name]\n" + " [-a_srs srs_def] [-spat xmin ymin xmax ymax]\n" + " [-clipsrc |WKT|datasource|spat_extent]\n" + " [-clipsrcsql sql_statement] [-clipsrclayer layer]\n" + " [-clipsrcwhere expression]\n" + " [-l layername]* [-where expression] [-sql select_statement]\n" + " [-txe xmin xmax] [-tye ymin ymax] [-outsize xsize ysize]\n" + " [-a algorithm[:parameter1=value1]*]\n" + " [-q]\n" + " \n" + "\n" + "Available algorithms and parameters with their's defaults:\n" + " Inverse distance to a power (default)\n" + " invdist:power=2.0:smoothing=0.0:radius1=0.0:radius2=0.0:" + "angle=0.0:max_points=0:min_points=0:nodata=0.0\n" + " Moving average\n" + " average:radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0\n" + " Nearest neighbor\n" + " nearest:radius1=0.0:radius2=0.0:angle=0.0:nodata=0.0\n" + " Various data metrics\n" + " :radius1=0.0:radius2=0.0:angle=0.0:min_points=0:nodata=0.0\n" + " possible metrics are:\n" + " minimum\n" + " maximum\n" + " range\n" + " count\n" + " average_distance\n" + " average_distance_pts\n"); System.exit(-1); } /* * Print algorithm and options. */ public static void PrintAlgorithmAndOptions(String algorithmAndOptions) { if (algorithmAndOptions == null) { System.out .println("Algorithm name: not selected, using default Inverse Distance"); } else { int firstColon = algorithmAndOptions.indexOf(':'); System.out.println("Algorithm name: " + algorithmAndOptions.substring(0, firstColon)); System.out.println("Options are \"" + algorithmAndOptions.substring(firstColon + 1) + "\""); } } /* * ProcessGeometry * * Extract point coordinates from the geometry reference and set the Z value * as requested. Test whether we are in the clipped region before * processing. */ public static void ProcessGeometry(Geometry point, Geometry clipSrc, int burnField, double burnValue, List X, List Y, List Z) { if (clipSrc != null && point.Within(clipSrc) == false) return; X.add(point.GetX()); Y.add(point.GetY()); if (burnField < 0) { Z.add(point.GetZ()); } else { Z.add(burnValue); } } /* * Process all the features in a layer selection, collecting geometries and * burn values. */ public static void ProcessLayer(Layer srcLayer, Dataset dstDS, Geometry clipSrc, int sizeX, int sizeY, int bandIndex, boolean[] isXExtentSet, boolean[] isYExtentSet, double[] minX, double[] maxX, double[] minY, double[] maxY, String burnAttribute, int type, String algorithmAndOptions, boolean quiet, ProgressCallback progressCallback) { /* * Get field index, and check. */ int burnFieldIndex = -1; if (burnAttribute != null) { burnFieldIndex = srcLayer.GetLayerDefn().GetFieldIndex( burnAttribute); if (burnFieldIndex == -1) { System.out.println("Failed to find field " + burnAttribute + " on layer " + srcLayer.GetLayerDefn().GetName()); return; } } /* * Collect the geometries from this layer, and build list of values to * be interpolated. */ Feature feature; List X = new ArrayList(); List Y = new ArrayList(); List Z = new ArrayList(); srcLayer.ResetReading(); while ((feature = srcLayer.GetNextFeature()) != null) { Geometry geometry = feature.GetGeometryRef(); if (geometry != null) { int geomtype = geometry.GetGeometryType() & (~ogrConstants.wkb25DBit); double burnValue = 0.0; if (burnFieldIndex >= 0) { burnValue = feature.GetFieldAsDouble(burnFieldIndex); } if (geomtype == ogr.wkbMultiPoint) { int geomIndex = 0; int geomCount = geometry.GetGeometryCount(); for (geomIndex = 0; geomIndex < geomCount; geomIndex++) { ProcessGeometry(geometry.GetGeometryRef(geomIndex), clipSrc, burnFieldIndex, burnValue, X, Y, Z); } } else { ProcessGeometry(geometry, clipSrc, burnFieldIndex, burnValue, X, Y, Z); } } feature.delete(); } if (X.size() == 0) { System.out.println("No point geometry found on layer " + srcLayer.GetLayerDefn().GetName() + ", skipping."); return; } /* * Compute grid geometry. */ if (isXExtentSet[0] == false) { minX[0] = X.get(0); maxX[0] = X.get(0); for (int i = 1; i < X.size(); i++) { if (minX[0] > X.get(i)) minX[0] = X.get(i); if (maxX[0] < X.get(i)) maxX[0] = X.get(i); } isXExtentSet[0] = true; } if (isYExtentSet[0] == false) { minY[0] = Y.get(0); maxY[0] = Y.get(0); for (int i = 1; i < Y.size(); i++) { if (minY[0] > Y.get(i)) minY[0] = Y.get(i); if (maxY[0] < Y.get(i)) maxY[0] = Y.get(i); } isYExtentSet[0] = true; } /* * Perform gridding. */ Double deltaX = (maxX[0] - minX[0]) / sizeX; Double deltaY = (maxY[0] - minY[0]) / sizeY; if (!quiet) { System.out.println("Grid data type is " + gdal.GetDataTypeName(type)); System.out.println("Grid size = (" + sizeX + " " + sizeY + ")."); System.out.println("Corner coordinates = (" + (minX[0] - deltaX / 2) + " " + (maxY[0] + deltaY / 2) + ") - (" + (maxX[0] + deltaX / 2) + " " + (minY[0] - deltaY / 2) + ")"); System.out.println("Grid cell size = (" + deltaX + " " + deltaY + ")."); System.out.println("Source point count = " + X.size() + " ."); PrintAlgorithmAndOptions(algorithmAndOptions); } Band band = dstDS.GetRasterBand(bandIndex); if (X.size() == 0) { Double val[] = new Double[1]; band.GetNoDataValue(val); if (val[0] != null) { band.Fill(val[0]); } else { band.Fill(0.0); } return; } int offsetX = 0; int offsetY = 0; int[] blockXSize = new int[1]; int[] blockYSize = new int[1]; band.GetBlockSize(blockXSize, blockYSize); int bufferSize = blockXSize[0] * blockYSize[0] * gdal.GetDataTypeSize(type) / 8; ByteBuffer data = ByteBuffer.allocateDirect(bufferSize); int blockIndex = 0; int blockCount = ((sizeX + blockXSize[0] - 1) / blockXSize[0]) * ((sizeY + blockYSize[0] - 1) / blockYSize[0]); GDALGridScaledProgress griddingProgress = null; for (offsetY = 0; offsetY < sizeY; offsetY += blockYSize[0]) { for (offsetX = 0; offsetX < sizeX; offsetX += blockXSize[0]) { int requestX = blockXSize[0]; if (offsetX + requestX > sizeX) { requestX = sizeX - offsetX; } int requestY = blockYSize[0]; if (offsetY + requestY > sizeY) { requestY = sizeY - offsetY; } /* * Reformat arguments */ double[][] points = new double[X.size()][3]; for (int i = 0; i < X.size(); i++) points[i][0] = X.get(i); for (int i = 0; i < Y.size(); i++) points[i][1] = Y.get(i); for (int i = 0; i < Z.size(); i++) points[i][2] = Z.get(i); /* * Create Scaled progress report */ if (quiet == false) { griddingProgress = new GDALGridScaledProgress( blockIndex * 1.0 / blockCount, (blockIndex + 1) * 1.0 / blockCount, progressCallback); } /* * Create Grid */ gdal.GridCreate(algorithmAndOptions, points, minX[0] + deltaX * offsetX, minX[0] + deltaX * (offsetX + requestX), minY[0] + deltaY * offsetY, minY[0] + deltaY * (offsetY + requestY), requestX, requestY, type, data, griddingProgress); /* * Write grid to raster output */ band.WriteRaster_Direct(offsetX, offsetY, requestX, requestY, requestX, requestY, type, data); if (quiet == false) { griddingProgress.delete(); } blockIndex++; } } } /* * LoadGeometry * * Read geometries from the given dataset using specified filters and * returns a collection of read geometries. */ static Geometry LoadGeometry(String srcDS, String srcSQL, String srcLyr, String srcWhere) { DataSource DS; Layer lyr; Feature feat; Geometry geom = null; DS = ogr.Open(srcDS, false); if (DS == null) { return null; } if (srcSQL != null) { lyr = DS.ExecuteSQL(srcSQL, null, null); } else if (srcLyr != null) { lyr = DS.GetLayerByName(srcLyr); } else { lyr = DS.GetLayer(0); } if (lyr == null) { System.err .println("Failed to identify source layer from datasource."); DS.delete(); return null; } if (srcWhere != null) { lyr.SetAttributeFilter(srcWhere); } while ((feat = lyr.GetNextFeature()) != null) { Geometry srcGeom = feat.GetGeometryRef(); if (srcGeom != null) { int srcType = srcGeom.GetGeometryType() & (~ogrConstants.wkb25DBit); if (geom == null) { geom = new Geometry(ogr.wkbMultiPolygon); } if (srcType == ogr.wkbPolygon) { geom.AddGeometry(srcGeom); } else if (srcType == ogr.wkbMultiPolygon) { int geomIndex = 0; int geomCount = srcGeom.GetGeometryCount(); for (geomIndex = 0; geomIndex < geomCount; geomIndex++) { geom.AddGeometry(srcGeom.GetGeometryRef(geomIndex)); } } else { System.err .println("FAILURE: Geometry not of polygon type."); if (srcSQL != null) { DS.ReleaseResultSet(lyr); } DS.delete(); return null; } } } if (srcSQL != null) { DS.ReleaseResultSet(lyr); } DS.delete(); return geom; } public static void main(String[] args) { String sourceFilename = null; String outputFilename = null; String outputFormat = "GTiff"; String layers = null; String burnAttribute = null; String where = null; int outputType = gdalconstConstants.GDT_Float64; String createOptions = null; boolean quiet = false; double[] minX = new double[1]; double[] maxX = new double[1]; double[] minY = new double[1]; double[] maxY = new double[1]; int sizeX = 0; int sizeY = 0; String SQL = null; boolean hasClipSrc = false; String clipSQL = null; String clipLayer = null; String clipWhere = null; String outputSRS = null; String algorithmAndOptions = null; Geometry clipSrc = null; Geometry spatialFilter = null; ProgressCallback progressCallback = null; String clipSrcDS = null; boolean[] isXExtentSet = new boolean[1]; boolean[] isYExtentSet = new boolean[1]; minX[0] = 0.0; maxX[0] = 0.0; minY[0] = 0.0; maxY[0] = 0.0; isXExtentSet[0] = false; isYExtentSet[0] = false; /* * Register GDAL and OGR format(s) */ gdal.AllRegister(); ogr.RegisterAll(); /* * Parse arguments */ args = ogr.GeneralCmdLineProcessor(args); if (args.length < 2) { Usage(); } for (int i = 0; i < args.length; i++) { if (args[i].equals("---utility_version")) { System.out.println("Running against GDAL " + gdal.VersionInfo()); return; } else if (args[i].equals("-of") && args.length > i) { outputFormat = args[++i]; } else if (args[i].equals("-q") || args[i].equals("-quiet")) { quiet = true; } else if (args[i].equals("-ot") && args.length > i) { outputType = gdal.GetDataTypeByName(args[++i]); if (outputType == gdalconstConstants.GDT_Unknown) { System.err.println("FAILURE: Unknown output pixel type: " + args[i]); Usage(); } } else if (args[i].equals("-txe") && args.length > i + 1) { minX[0] = Double.parseDouble(args[++i]); maxX[0] = Double.parseDouble(args[++i]); isXExtentSet[0] = true; } else if (args[i].equals("-tye") && args.length > i + 1) { minY[0] = Double.parseDouble(args[++i]); maxY[0] = Double.parseDouble(args[++i]); isYExtentSet[0] = true; } else if (args[i].equals("-outsize") && args.length > i + 1) { sizeX = Integer.parseInt(args[++i]); sizeY = Integer.parseInt(args[++i]); } else if (args[i].equals("-co") && args.length > i) { if (createOptions == null) { createOptions = args[++i]; } else { createOptions += ':' + args[++i]; } } else if (args[i].equals("-zfield") && args.length > i) { burnAttribute = args[++i]; } else if (args[i].equals("-where") && args.length > i) { where = args[++i]; } else if (args[i].equals("-l") && args.length > i) { if (layers == null) { layers = args[++i]; } else { layers += ':' + args[++i]; } } else if (args[i].equals("-sql") && args.length > i) { SQL = args[++i]; } else if (args[i].equals("-spat") && args.length > i + 3) { double clipMinX = Double.parseDouble(args[i + 1]); double clipMinY = Double.parseDouble(args[i + 2]); double clipMaxX = Double.parseDouble(args[i + 3]); double clipMaxY = Double.parseDouble(args[i + 4]); i += 4; clipSrc = new Geometry(ogr.wkbPolygon); clipSrc.AddPoint(clipMinX, clipMinY); clipSrc.AddPoint(clipMinX, clipMaxY); clipSrc.AddPoint(clipMaxX, clipMaxY); clipSrc.AddPoint(clipMaxX, clipMinY); clipSrc.AddPoint(clipMinX, clipMinY); } else if (args[i].equals("-clipsrc") && args.length > i) { hasClipSrc = true; try { double clipMinX = Double.parseDouble(args[i + 1]); double clipMinY = Double.parseDouble(args[i + 2]); double clipMaxX = Double.parseDouble(args[i + 3]); double clipMaxY = Double.parseDouble(args[i + 4]); i += 4; clipSrc = new Geometry(ogr.wkbPolygon); clipSrc.AddPoint(clipMinX, clipMinY); clipSrc.AddPoint(clipMinX, clipMaxY); clipSrc.AddPoint(clipMaxX, clipMaxY); clipSrc.AddPoint(clipMaxX, clipMinY); clipSrc.AddPoint(clipMinX, clipMinY); } catch (NumberFormatException e) { if (args[i].substring(0, 6).equals("POLYGON") || args[i].substring(0, 11).equals("MULTIPOLYGON")) { clipSrc = ogr.CreateGeometryFromWkt(args[++i]); } else if (args[i].equals("spat_extent")) { ++i; } else { clipSrcDS = args[++i]; } } } else if (args[i].equals("-clipsrcsql") && args.length > i) { clipSQL = args[++i]; } else if (args[i].equals("-clipsrclayer") && args.length > i) { clipLayer = args[++i]; } else if (args[i].equals("-clipsrcwhere") && args.length > i) { clipLayer = args[++i]; } else if (args[i].equals("-a_srs") && args.length > i) { SpatialReference enteredSRS = new SpatialReference(); if (enteredSRS.SetFromUserInput(args[i + 1]) != 0) { System.err.println("Failed to process SRS definition: " + args[i + 1]); Usage(); } i++; outputSRS = enteredSRS.ExportToWkt(); } else if (args[i].equals("-a") && args.length > i) { if (algorithmAndOptions == null) { algorithmAndOptions = args[++i]; } else { algorithmAndOptions += ':' + args[++i]; } } else if (args[i].substring(0, 1).equals("-")) { System.err.println("FAILURE: Option " + args[i] + "incomplete, or not recognised"); Usage(); } else if (sourceFilename == null) { sourceFilename = args[i]; } else if (outputFilename == null) { outputFilename = args[i]; } else { Usage(); } } if (sourceFilename == null || outputFilename == null || (SQL == null && layers == null)) { Usage(); } /* * Open Input */ if (hasClipSrc && clipSrcDS != null) { clipSrc = LoadGeometry(clipSrcDS, clipSQL, clipLayer, clipWhere); if (clipSrc == null) { System.out.println("FAILURE: cannot load source clip geometry"); Usage(); } } else if (hasClipSrc && clipSrcDS == null) { clipSrc = spatialFilter.Clone(); if (clipSrc == null) { System.out .println("FAILURE: " + "-clipsrc must be used with -spat option or \n" + "a bounding box, WKT string or datasource must be specified"); Usage(); } } /* * Find the output driver. */ Driver driver = gdal.GetDriverByName(outputFormat); if (driver == null) { System.out.println("FAILURE: Output driver '" + outputFormat + "' not recognized."); Usage(); } /* * Open input datasource. */ DataSource srcDS = ogr.Open(sourceFilename); if (srcDS == null) { System.out.println("Unable to open input datasource '" + sourceFilename); System.out.println(gdal.GetLastErrorMsg()); System.exit(3); } /* * Create target raster file. */ Dataset dstDS = null; String[] layerList = layers.split(":"); int layerCount = layerList.length; int bandCount = layerCount; if (SQL != null) { bandCount++; } if (sizeX == 0) { sizeX = 256; } if (sizeY == 0) { sizeY = 256; } String[] optionList = createOptions == null ? null : createOptions .split(":"); dstDS = driver.Create(outputFilename, sizeX, sizeY, 1, outputType, optionList); if (dstDS == null) { System.out.println("Unable to create dataset '" + outputFilename); System.out.println(gdal.GetLastErrorMsg()); System.exit(3); } /* * Use terminal progress report */ if (quiet == false) { progressCallback = new TermProgressCallback(); } /* * Process SQL request. */ if (SQL != null) { Layer srcLayer = srcDS.ExecuteSQL(SQL); if (srcLayer != null) { ProcessLayer(srcLayer, dstDS, clipSrc, sizeX, sizeY, 1, isXExtentSet, isYExtentSet, minX, maxX, minY, maxY, burnAttribute, outputType, algorithmAndOptions, quiet, progressCallback); } } /* * Process each layer. */ for (int i = 0; i < layerList.length; i++) { Layer srcLayer = srcDS.GetLayerByName(layerList[i]); if (srcLayer == null) { System.out.println("Unable to find layer '" + layerList[i]); continue; } if (where != null) { if (srcLayer.SetAttributeFilter(where) != gdalconstConstants.CE_None) { break; } } if (spatialFilter != null) { srcLayer.SetSpatialFilter(spatialFilter); } if (outputSRS == null) { SpatialReference srs = srcLayer.GetSpatialRef(); if (srs != null) { outputSRS = srs.ExportToWkt(); } } ProcessLayer(srcLayer, dstDS, clipSrc, sizeX, sizeY, i + 1 + bandCount - layerCount, isXExtentSet, isYExtentSet, minX, maxX, minY, maxY, burnAttribute, outputType, algorithmAndOptions, quiet, progressCallback); } /* * Apply geotransformation matrix. */ double[] geoTransform = new double[6]; geoTransform[0] = minX[0]; geoTransform[1] = (maxX[0] - minX[0]) / sizeX; geoTransform[2] = 0.0; geoTransform[3] = minY[0]; geoTransform[4] = 0.0; geoTransform[5] = (maxY[0] - minY[0]) / sizeY; dstDS.SetGeoTransform(geoTransform); /* * Apply SRS definition if set. */ if (outputSRS != null) { dstDS.SetProjection(outputSRS); } /* * Cleanup. */ srcDS.delete(); dstDS.delete(); gdal.GDALDestroyDriverManager(); } } class GDALGridScaledProgress extends ProgressCallback { private double pctMin; private double pctMax; private ProgressCallback mainCbk; public GDALGridScaledProgress(double pctMin, double pctMax, ProgressCallback mainCbk) { this.pctMin = pctMin; this.pctMax = pctMax; this.mainCbk = mainCbk; } public int run(double dfComplete, String message) { return mainCbk.run(pctMin + dfComplete * (pctMax - pctMin), message); } };