java调用GDAL实现栅格数据的重采样的一种方法

1.关于重采样

1.1概念

1.2用途

1.3常见算法

2.关于GDAL

2.1GDAL中的重采样算法

3.实现重采样

3.1思路

3.2完整代码

3.3使用QGIS验证效果

1.关于重采样

1.1概念

重采样是以原始图像的像元值或者导出的值填充到新的图像的每个像元的的过程。

1.2用途

在地理信息栅格数据的处理中有着比较重要的用途，如地图制图中不同比例尺下需要有不同的分辨率显示；在卫星遥感影像的数据的处理中，重采样也是必不可少的一个环节；在空间分析时，保持不同来源的数据的分辨率一致也是必要的操作。

1.3常见算法

重采样常用的算法有：邻近点插值、双线性插值、三次卷积插值等。他们各自有其应用场景。这里暂时不对这些算法展开深入讨论。

2.关于GDAL

GDAL的Warp程序不仅能实现栅格数据的重投影和镶嵌，也可以实现栅格数据的重采样。本文简单介绍一个使用java调用GDAL的API实现栅格数据的重投影的方法。

2.1GDAL中的重采样算法

根据GDAL的文档介绍，其支持的重采样算法有以下类型：

near为邻近点插值，它是将最近的原始像元值赋值给新像元。

bilinear为双线性插值，它是通过4个相邻像元的距离加权平均值估算新像元值的方法。

cubic为三次卷积插值，它是通过16个相邻像元的距离加权平均值估算新像元值的方法。

其他方法暂时不介绍。

3.实现重采样

3.1思路

在GDAL的java绑定库中，gdal.Warp允许传入Dataset、Dataset[]、WarpOptions对像作为入参来进行调用，结合GDAL中对C语言的GDALWarp函数的描述文档，再经过验证，可以确认Dataset为要输出的栅格数据集，Dataset[]是输入的栅格数据集，包含原栅格和要输出的栅格。

  public static int Warp(Dataset dstDS, Dataset[] object_list_count, WarpOptions warpAppOptions) {
    return gdalJNI.Warp__SWIG_2(Dataset.getCPtr(dstDS), dstDS, object_list_count, WarpOptions.getCPtr(warpAppOptions), warpAppOptions);
  }

根据上述参数，我们可以根据输入文件路径读取源文件为Dataset，根据输出文件路径创建Dataset，创建一个Vector对象来构建WarpOptions，Vector的内容按照gdalwarp命令行程序的参数要求来输入，最后调用gdal.Warp来给新创建的数据集填充数据，填充数据后保存即可。

以下代码是对上述文字描述的代码的一个示例（省去了创建数据集的代码）。

        Vector<String> vector = new Vector<>();
        vector.add("-r");
        vector.add("bilinear");//执行双线性插值法重采样
        WarpOptions warpOptions = new WarpOptions(vector);
        gdal.Warp(datasetWarp, new Dataset[]{dataset}, warpOptions);

3.2完整代码

以下代码是一个完整的方法，输入要重采样的栅格数据的路径文本、要输出的栅格数据的路径文本、输入东西、南北方向的分辨率，其中南北方向的分辨率需要输入负数。

 /**
     * 对栅格重采样
     *
     * @param inputPath  输入的栅格路径
     * @param outputPath 输出的栅格路径
     * @param r1         东西方向像元大小（地图单位）
     * @param r2         南北方向像元大小（地图单位）
     */
    public static void reSample(String inputPath, String outputPath, double r1, double r2) {
        // 打开输入栅格文件
        Dataset dataset = gdal.Open(inputPath, gdalconst.GA_ReadOnly);
        double[] gt1 = dataset.GetGeoTransform();

        int xSize = new Double(dataset.getRasterXSize() * gt1[1] / r1).intValue();
        int ySize = new Double(dataset.getRasterYSize() * gt1[5] / r2).intValue();
        ySize = Math.abs(ySize);

        double[] gt2 = {gt1[0], r1, 0.0, gt1[3], 0.0, r2};
        Driver driver = gdal.GetDriverByName("GTiff");
        Dataset datasetWarp = driver.Create(outputPath, xSize, ySize, gdalconst.GA_Update, gdalconst.GDT_Float32);
        datasetWarp.SetGeoTransform(gt2);
        datasetWarp.SetSpatialRef(dataset.GetSpatialRef());
        Vector<String> vector = new Vector<>();
        vector.add("-r");
        vector.add("bilinear");//执行双线性插值法重采样
        WarpOptions warpOptions = new WarpOptions(vector);
        gdal.Warp(datasetWarp, new Dataset[]{dataset}, warpOptions);

        datasetWarp.FlushCache();
        datasetWarp.delete();
        dataset.delete();
    }