gdal java 无法获取坐标_Java版gdal 初探

最新更新：2019.11.08：四.Spark并行化规整裁切

刚看了没几天的geotrellis就被叫去看gdal了，还好以前有一丢丢基础，借着这个机会复习一下。

一. gdal部署

部署方法可以参见另外两篇文章：windows下java版gdal部署和linux下java版gdal部署。

其中，在windows的部署文章中提到，要复制gdalconstjni.dll、gdaljni.dll、ogrjni.dll、osrjni.dll这四个dll，但本文使用的gdal版本是2.3.2，这四个dll已经被整合成了一个gdalalljni.dll。

二 .影像读取及基本信息查看

首先扔一个gdal的java API文档，该文档适用于gdal 1.7.0或以上版本。

首先是所有gdal程序都需要的注册语句：

gdal.AllRegister();

之前的版本还需要加这句话来支持中文路径：

gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8","YES");

但在2.3.2中，我没有加这句话，也能正常地在中文路径下读写影像。读取影像基本信息的代码如下：

//读取影像

Dataset rds = gdal.Open("影像路径", gdalconst.GA_ReadOnly);

//宽、高、波段数

int x = rds.getRasterXSize();

int y = rds.getRasterYSize();

int b = rds.getRasterCount();

//六参数信息

double[] geoTransform = rds.GetGeoTransform();

//影像左上角投影坐标

double[] ulCoord = new double[2];

ulCoord[0] = geoTransform[0];

ulCoord[1] = geoTransform[3];

//影像右下角投影坐标

double[] brCoord = new double[2];

brCoord[0] = geoTransform[0] + x * geoTransform[1] + y * geoTransform[2];

brCoord[1] = geoTransform[3] + x * geoTransform[4] + y * geoTransform[5];

//影像投影信息

String proj = rds.GetProjection();

这里有个数组geoTransform，容量为6，代表的是仿射变换六参数，其含义如下：

geoTransform[0]：左上角x坐标

geoTransform[1]：东西方向空间分辨率

geoTransform[2]：x方向旋转角

geoTransform[3]：左上角y坐标

geoTransform[4]：y方向旋转角

geoTransform[5]：南北方向空间分辨率

影像中任意一个像素的坐标可以用下式计算：

任意像素的坐标计算公式

其中六参数分别为{dx, a11, a12, dy, a21, a22}，正常的影像，正北方向朝上，两个旋转角的值都是0。

三. 影像裁切

在gdal中，影像裁切的基本思想，是把原始影像中想要保留的部分进行另存为。因此这部分要做的事情，就是在定好要裁切的部分后，读取相应的部分并保存即可，其中涉及到一些简单的位置计算。代码如下：

gdal.AllRegister();

Dataset rds = gdal.Open("影像路径", gdalconst.GA_ReadOnly);

//宽、高、波段数

int b = rds.getRasterCount();

//从波段中获取影像的数据类型，gdal中波段索引从1开始

int dataType = rds.GetRasterBand(1).GetRasterDataType();

//六参数信息

double[] geoTransform = rds.GetGeoTransform();

//设置要裁剪的起始像元位置，以及各方向的像元数

//这里表示从像元(5000, 5000)开始，x方向和y方向各裁剪1000个像元

//最终结果就是一幅1000*1000的影像

int startX = 5000;

int startY = 5000;

int clipX = 1000;

int clipY = 1000;

//计算裁剪后的左上角坐标

geoTransform[0] = geoTransform[0] + startX * geoTransform[1] + startY * geoTransform[2];

geoTransform[3] = geoTransform[3] + startX * geoTransform[4] + startY * geoTransform[5];

//创建结果图像

Driver driver = gdal.GetDriverByName("GTIFF");

Dataset outputDs = driver.Create("D:\\Javaworkspace\\gdal\\output\\clip1.tif", clipX, clipY, b, dataType);

outputDs.SetGeoTransform(geoTransform);

outputDs.SetProjection(rds.GetProjection());

//按band读取

for(int i = 0; i < clipY; i++){

//按行读取

for(int j = 1; j <= b; j++){

Band orgBand = rds.GetRasterBand(j);

int[] cache = new int[clipX];

//从位置x开始，只读一行

orgBand.ReadRaster(startX, startY + i, clipX, 1, cache);

Band desBand = outputDs.GetRasterBand(j);

//从左上角开始，只写一行

desBand.WriteRaster(0, i, clipX, 1, cache);

desBand.FlushCache();

}

}

//释放资源

rds.delete();

outputDs.delete();

除了Band对象中的ReadRaster和WriteRaster接口外，Dataset对象也有ReadRaster和WriteRaster的接口，但是用法有所不同。要使用Dataset中的接口，需要先构建波段数组，表示在读写过程中对哪些波段进行操作。

int[] bands = new int[]{1};

然后在实例化cache数组的时候，数组容量就不是简单的要读写的影像窗口的x*y，而是x*y*b，b表示波段数。而数组的类型也需要根据影像的数据类型进行调整。gdal中的数据类型以常量的形式存放在org.gdal.gdalconst这个类中。从接口文档中能查到有这么些个：

gdal数据类型

正常来讲，可以使用位数更高的数据类型的数组来读取位数较低的数据。比如在我自己的影像中，datatype的值为2，对应的数据类型为GDT_UInt16，即无符号16位整型，但我在使用int类型数组来读数据的时候，报了数据类型不匹配的错(仅限Dataset的接口)，但是输出的影像好像没有问题。报错信息如下：

ERROR 1: Java array type is not compatible with GDAL data type

当我改用short类型的数组来读写数据后，这个报错就消失了，不是很懂gdal的这一波操作。最后利用Dataset对象读写数据的代码如下，也是按行来执行：

int[] bands = new int[]{1};

for (int i = 0; i < clipY; i++)

{

short[] cache = new short[clipX * b];

rds.ReadRaster(startX, startY +