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();
}
这个函数中dataset是原始数据，boundingBox是裁剪矩形范围，当然这个方法，无法裁剪不规则多边形，也就是它只能裁剪矩形规则范围。
相比使用shape文件裁剪效率，这个直接拷贝像素的效率更高，当然这些裁剪方法中，还是直接使用gdalwarp.exe的-te参数效率最高，但是
奈何gdal的java api中好像没有适配这个参数，不知道为什么？
好了今天的内容就讲解到这里，如有问题，欢迎评论区讨论。