Java 调用gdal API(二)——栅格裁剪_Java教程

Java 调用gdal API(二)——栅格裁剪

gdal可以说是GIS数据处理比较好的工具之一，虽然也提供了Java API，但是官方文档确实太过简单，用起来确实太难受，每次都需要去参考对应的C++api，然后在对应使用。

因此小编决定从这篇文章开始，将自己以前以及今后用到的API，都记录下，以帮助更多从事java GIS服务开发者。

今天主要来记录栅格数据的裁剪功能，gdal中提供了多种裁剪方法，今天主要介绍以下两种：

第一中是通过gdal自身的Warp方法，其实现有以下几个接口；

public static int Warp(Dataset dstDS, Dataset[] object_list_count, WarpOptions warpAppOptions, ProgressCallback callback)； public static int Warp(Dataset dstDS, Dataset[] object_list_count, WarpOptions warpAppOptions)； public static Dataset Warp(String dest, Dataset[] object_list_count, WarpOptions warpAppOptions, ProgressCallback callback)； public static Dataset Warp(String dest, Dataset[] object_list_count, WarpOptions warpAppOptions)；

其实这个几个接口基本差不多

1、dstDs/dest 其实就是目标数据，两种方式，要么传入目标路径，要么是目标Dataset

2、object_list_count 这个是源数据的Dataset 数据

3、warpAppOptions 这个参数很重要，是我们裁剪需要传递的重要参数，后面代码会专门解释

4、callback 处理进度

下面我们结合最后一个方法，来处理一个栅格裁剪

private static void cutTif() throws FactoryException, IOException, ParseException {
        Dataset dataSet = gdal.Open("E:\\数据\\影像数据\\tiff\\xian18.tif");
        //通过最大最小经纬度生成一个面矢量
        createCutShapeRegion(108.583740,33.928528,108.883740,34.260242);
        Vector vector = new Vector();
        vector.add("-cutline");
        vector.add("D:\\data\\cut2.shp");
        vector.add("-crop_to_cutline");
        Dataset ds = gdal.Warp("D:\\data\\xian18_1.tif",new Dataset[]{dataSet},new WarpOptions(vector));
        if (ds != null){
            System.out.println("成功！");
        }
    }

createCutShapeRegion是用通过bbox来生成shape数据的，我这边是用矩形来裁剪的，当然你也可以用多边形，通过传入多边形的顶点来生层一个面矢量，记住一定是面矢量，

D:\\data\\cut2.shp是createCutShapeRegion生成的矢量文件，最主要的两个参数一个是-cutline 一个是-crop_to_cutline，其中-cutline后面跟裁剪的矢量文件，-crop_to_cutline不需要具体值，这些具体参数可参考gdal的warp的参数介绍，其实和安装版gdalwarp的参数是基本一致的。

但是这里说明下，并不是gdalwarp的所有参数都在gdal的java api中进行了适配，就比如gdalwarp这个命令中除了上面使用shp作为裁剪之外，还有一个-te参数，这个参数后面跟一个经纬度的范围
主要作用就是将原始数据裁剪为制定的范围，就比如我们使用安装的gdal中可以这样调用裁剪
gdalwarp.exe -te 108.583740 33.928528 108.883740 34.260242  E:\\数据\\影像数据\\tiff\\xian18.tif D:\\data\\xian18_2.tif
然而这个-te参数如果使用代码来调用，就会报没有-te这个option，就比如下面的这个例子其实是跑不通的

public static void main(String[] args) throws NoSuchFieldException, IllegalAccessException, FactoryException, IOException, ParseException {
        gdal.AllRegister();
        Dataset dataSet = gdal.Open("E:\\数据\\影像数据\\tiff\\xian18.tif");
        Vector vector = new Vector();
        vector.add("-te");
        vector.add("108.583740 33.928528 108.883740 34.260242"); //裁剪范围
        vector.add("-crop_to_cutline");
        Dataset ds = gdal.Warp("D:\\data\\xian18_1.tif",new Dataset[]{dataSet},new WarpOptions(vector));
    }

上述代码是跑不通的，虽然安装板的gdalwarp可以使用-te来裁剪，而且效率比上面使用shp要更快，但是java代码的api却没有集成这个-te参数，当然我不知道是我调用的问题，还是确实没有适配

如果有了解的朋友，麻烦告知一声。

除了上述使用gdalwarp的接口外，我们还可以通过使用driver来快速裁剪出一个tif,具体代码如下

private static void clipRaster(Dataset dataset, Envelope boundingBox){
        //获取当前数据范围
        ReferencedEnvelope originalEnvelope = (ReferencedEnvelope) GdalDataAnalysis.readBboxFromFile(dataset);
        //判断是否相交，不想交则直接返回
        if (!originalEnvelope.intersects(boundingBox)){
            return;
        }
        //求请求与原始数据的范围
        ReferencedEnvelope interEnvelope = originalEnvelope.intersection(boundingBox);
        //原始数据上偏移位置
        int offset_x = (int) ((interEnvelope.getMinX() - originalEnvelope.getMinX()) / GdalDataAnalysis.getResoulution(dataset,true));
        int offset_y = (int) ((interEnvelope.getMaxY() - originalEnvelope.getMaxY()) / GdalDataAnalysis.getResoulution(dataset,false));
        //需要多少个原始像素
        int block_xsize = (int) ((interEnvelope.getMaxX() - interEnvelope.getMinX()) / GdalDataAnalysis.getResoulution(dataset,true));
        int block_ysize = (int)Math.abs((interEnvelope.getMaxY() - interEnvelope.getMinY()) / GdalDataAnalysis.getResoulution(dataset,false));
        //实际瓦片的长度，裁剪数据的分辨率和原始图像一致
        int image_Xbuf = block_xsize;
        int image_Ybuf =  block_ysize;
        //写入瓦片的偏移位置
        int imageOffSetX = 0;
        int imageOffSetY = 0;
        //循环读写
        Driver driver = gdal.GetDriverByName("GTiff");
        int dataType = dataset.GetRasterBand(1).getDataType();
        int bands = dataset.getRasterCount();
        byte[] pixelArray = new byte[image_Xbuf*image_Ybuf];
        Dataset dsDataset = driver.Create("D:\\data\\xian18_1.tif",image_Xbuf,image_Ybuf,bands, dataType);
        for (int i = 1;i<=bands;i++){
            Band band = dataset.GetRasterBand(i);
            band.ReadRaster(offset_x,offset_y,block_xsize,block_ysize,image_Xbuf,image_Ybuf,dataType,pixelArray);
            //写数据
            dsDataset.GetRasterBand(i).WriteRaster(imageOffSetX,imageOffSetY,block_xsize,block_ysize,pixelArray);
            //清理删除缓存
            dsDataset.GetRasterBand(i).FlushCache();
            dsDataset.GetRasterBand(i).delete();
        }
        dsDataset.SetProjection(dataset.GetProjection());
        double[] rsGeoTransform = dataset.GetGeoTransform();
        rsGeoTransform[0] = interEnvelope.getMinX();
        rsGeoTransform[3] = interEnvelope.getMaxY();
        dsDataset.SetGeoTransform(rsGeoTransform);
        dsDataset.FlushCache();
        dsDataset.delete();
    }