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

发布时间 2023-07-23 17:32:45作者: 开放GIS

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

这个函数中dataset是原始数据,boundingBox是裁剪矩形范围,当然这个方法,无法裁剪不规则多边形,也就是它只能裁剪矩形规则范围。

相比使用shape文件裁剪效率,这个直接拷贝像素的效率更高,当然这些裁剪方法中,还是直接使用gdalwarp.exe的-te参数效率最高,但是

奈何gdal的java api中好像没有适配这个参数,不知道为什么?

好了今天的内容就讲解到这里,如有问题,欢迎评论区讨论。