/******************************************************************************
 * $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 <xmin ymin xmax ymax>|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"
                        + "    <src_datasource> <dst_filename>\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"
                        + "        <metric name>: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<Double> X, List<Double> Y,
            List<Double> 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<Double> X = new ArrayList<Double>();
        List<Double> Y = new ArrayList<Double>();
        List<Double> Z = new ArrayList<Double>();

        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);
    }
};