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 + i, clipX, 1, clipX, 1, dataType, cache, bands);
outputDs.WriteRaster(0, i, clipX, 1, clipX, 1, dataType, cache, bands);
outputDs.FlushCache();
}
补充：
后来发现，有一种比较通用的方式可以读取各种数据类型的影像，即使用byte类型的cache数组，但其长度还需要再乘以目标数据类型与byte类型之间的倍数。举个例子，假设要存储一个像素的short类型的数据，正常来讲只需要实例化一个长度为1的short数组；如果用byte型，一个byte是8位，而一个short是16位，因此应该实例化一个长度1*2=2的byte数组。gdal中也提供了接口可以用来计算这个倍数，代码如下：
// 首先获取影像数据类型 int dataType = dataset.GetRasterBand(1).GetRasterDataType(); // 计算该类型的字节数 int typeSize = gdal.GetDataTypeSize(dataType); // 一个byte8位，计算倍数只需除以8即可 int byteToType = typeSize / 8; // 假设要读取宽x，高y，波段数b的数据 byte[] datas = new byte[x * y * b * byteToType]; // 读取时只需将cache数组替换，读取的宽度、高度、数据类型等参数都不需要更改，写入同理 dataset.Read(0, 0, x, y, x, y, dataType, datas, bandList);
这种方法在不对数据进行具体计算，且又想在不确定数据的类型时使程序变得通用的情况下相当有效。但对于某些压缩过的数据，用short存储，需要通过元数据信息进行计算复原的话就不适用了。
四.Spark并行化规整裁切
这次会使用gdal，主要就是为了这个需求。利用spark，将影像裁切成一个个规整的方块。比如我设定x方向和y方向的切片数量分别为5的话，结果应该是这样：
影像规整裁切结果
并且，这些影像都应该带有位置信息，在ENVI或ArcGIS等平台中加载，能够拼成完整的一幅影像。通过上面裁切的demo，想要实现这个功能并不难，只需要多套两层循环即可。现在需要考虑的，是如何利用Spark，进行分布式的影像裁切。
放了一波寒假，一段时间没碰这个了。。。。。。现在考虑的思路是，建立一个格网索引，将影像分发到它们所在的格网中，并裁出格网与影像相交的部分，最后在每一个格网内将影像进行重叠、拼接，重叠部分则计算平均值。
更新，很长一段时间没有在看gdal了，这次把这部分内容的坑填上吧。Spark拥有众多算子，想要实现一个目标可以有很多种写方法（光是wordcount我所知道的就有2~3种）。因此，对于影像的规整裁切，我所采用的思路、方法可能不是最优的，但这里也先贴出来供参考。
- 输入及元数据信息提取
Spark不能直接读取GeoTiff数据，这里还是要靠GDAL。在我的思路中输入影像应该是可以有多幅的，最终的结果应该呈现出所有输入影像拼接后再裁切的样子（当然实际上并没有拼接的过程）。所以首先统计了输入影像的各类元数据信息，来为进一步切分做准备。
- 首先定义一个需要获得的信息的类
public class DatasetInfo implements Serializable {
public String projection; // 投影WKT
public int dataType; // 数据类型
public int bandCount; // 波段数
public double xPixel; // x方向分辨率
public double yPixel; // y方向分辨率
public double xRotate; // x方向旋转因子
public double yRotate; // y 方向旋转因子
public double noData; // NODATA值
public Extent extent; // 外包框
public DatasetInfo(String projection, int dataType, int bandCount, double xPixel, double yPixel,
double xRotate, double yRotate, double noData, Extent extent) {
this.projection = projection;
this.dataType = dataType;
this.bandCount = bandCount;
this.xPixel = xPixel;
this.yPixel = yPixel;
this.xRotate = xRotate;
this.yRotate = yRotate;
this.noData = noData;
this.extent = extent;
}
}
// 外包框类
public class Extent implements Serializable {
public double minX;
public double minY;
public double maxX;
public double maxY;
public Extent(double minX, double minY, double maxX, double maxY) {
this.minX = minX;
this.minY = minY;
this.maxX = maxX;
this.maxY = maxY;
}
// 外包框重叠
public Extent overlap(Extent extent){
double newMinX = Math.max(minX, extent.minX);
double newMaxX = Math.min(maxX, extent.maxX);
double newMinY = Math.max(minY, extent.minY);
double newMaxY = Math.min(maxY, extent.maxY);
if (newMinX >= newMaxX || newMinY >= newMaxY) {
throw new RuntimeException("Not overlap extents!");
}
return new Extent(newMinX, newMinY, newMaxX, newMaxY);
}
// 外包框合并
public Extent union(Extent extent) {
double newMinX = Math.min(minX, extent.minX);
double newMinY = Math.min(minY, extent.minY);
double newMaxX = Math.max(maxX, extent.maxX);
double newMaxY = Math.max(maxY, extent.maxY);
return new Extent(newMinX, newMinY, newMaxX, newMaxY);
}
}
- 根据输入路径参数（输入文件路径的List）整合元数据
思路是先利用map算子读取每一幅输入影像的元数据信息，再在reduce中统合。统合原则是，外包框进行合并，其他信息则要确保一致，否则抛出异常。
// 持久化路径RDD，以后还要使用
JavaRDD<String> pathsRDD = sc.parallelize(inputs).cache();
DatasetInfo datasetInfo = pathsRDD.map(file -> {
// 首先读取每一幅影像的元数据
// 在每一个Executor中都需要先进行注册
gdal.AllRegister();
Dataset ds = gdal.Open(file, gdalconst.GA_ReadOnly);
double[] geoTransform = ds.GetGeoTransform();
// 构造geoTransform数组，其含义见上文
int width = ds.getRasterXSize();
int height = ds.getRasterYSize();
double leftUpX = geoTransform[0];
double leftUpY = geoTransform[3];
double rightBottomX = leftUpX + width * geoTransform[1] + height * geoTransform[2];
double rightBottomY = leftUpY + width * geoTransform[4] + height * geoTransform[5];
Extent extent = new Extent(leftUpX, rightBottomY, rightBottomX, leftUpY);
String proj = ds.GetProjection();
int dataType = ds.GetRasterBand(1).GetRasterDataType();
int bandCount = ds.getRasterCount();
Double[] val = new Double[1];
// 获取NoData值，若没有，默认置0
ds.GetRasterBand(1).GetNoDataValue(val);
double noData = DEFAULT_NO_DATA_VALUE;
if (val[0] != null) {
noData = val[0];
}
ds.delete();
return new DatasetInfo(proj, dataType, bandCount, geoTransform[1], geoTransform[5],
geoTransform[2], geoTransform[4], noData, extent);
}).reduce((info1, info2) -> {
// 利用reduce算子整合所有元数据
// 确保外包框以外的信息都一样
if (!info1.projection.equals(info2.projection)
|| info1.dataType != info2.dataType
|| info1.bandCount != info2.bandCount
|| info1.xPixel != info2.xPixel
|| info1.yPixel != info2.yPixel
|| info1.xRotate != info2.xRotate
|| info1.yRotate != info2.yRotate
|| info1.noData != info2.noData) {
throw new RuntimeException
("Input images should have same projection, dataType, bands, resolution, rotation and noData");
}
// 合并外包框
Extent extent = info1.extent.union(info2.extent);
return new DatasetInfo(info1.projection, info1.dataType, info1.bandCount,
info1.xPixel, info1.yPixel, info2.xRotate, info2.yRotate, info1.noData, extent);
});
}
- 利用统合元数据信息构建格网索引。假设格网是个正方形，其分片数（slices）为格网数量的平方。比如你最后想输出5x5的切片，那么格网数量就是25，分片数就是5。
public class GridCell implements Serializable {
// 格网范围
private Extent extent;
// 分片数
private int slices;
// 每个大格子的平均边长
private int cellXSize;
private int cellYSize;
// 格网表示成影像时其长宽像元数
private int allXSize;
private int allYSize;
// 影像的分辨率
private double xPixel;
private double yPixel;
public GridCell(Extent extent, int slices, double xPixel, double yPixel) {
this.extent = extent;
this.slices = slices;
this.xPixel = xPixel;
this.yPixel = yPixel;
this.allXSize = (int) ((extent.maxX - extent.minX) / xPixel);
this.allYSize = (int) ((extent.minY - extent.maxY) / yPixel);
this.cellXSize = (int) Math.ceil(allXSize * 1.0 / slices);
this.cellYSize = (int) Math.ceil(allYSize * 1.0 / slices);
}
public int getXIndex(double x) {
int widthSize = (int) ((x - extent.minX) / xPixel);
int xIndex = widthSize / cellXSize;
return Math.min(xIndex, slices - 1);
}
public int getYIndex(double y) {
int heightSize = (int) ((y - extent.maxY) / yPixel);
int yIndex = heightSize / cellYSize;
return Math.min(yIndex, slices - 1);
}
public Extent getCellExtent(int xIndex, int yIndex) {
if (xIndex >= slices || yIndex >= slices) {
return null;
}
int minXIndex = xIndex * cellXSize;
int maxXIndex = Math.min(minXIndex + cellXSize, allXSize);
int minYIndex = yIndex * cellYSize;
int maxYIndex = Math.min(minYIndex + cellYSize, allYSize);
return new Extent(extent.minX + minXIndex * xPixel,
extent.maxY + maxYIndex * yPixel,
extent.minX + maxXIndex * xPixel,
extent.maxY + minYIndex * yPixel);
}
}
// 利用元数据信息构造格网
GridCell gridCell = new GridCell(datasetInfo.extent, slices, datasetInfo.xPixel, datasetInfo.yPixel);
看代码可能有点抽象，具体来讲，假设我有一个5x5的影像，想切成3x3的格网
GridCell各属性含义
- 进行裁切
准备工作都做完了，接下来开始进行裁切。首先，为了减少对象的复制次数，把必要的东西先进行广播。
// 构建读写tiff时需要用到的波段数组
int[] bands = new int[datasetInfo.bandCount];
for (int i = 0; i < bands.length; i++) {
bands[i] = i + 1;
}
Broadcast<Integer> byteSizeBroadcast = sc.broadcast(BYTESIZE);
Broadcast<String> outputBroadcast = sc.broadcast(outputDir);
Broadcast<DatasetInfo> datasetInfoBroadcast = sc.broadcast(datasetInfo);
Broadcast<GridCell> gridCellBroadcast = sc.broadcast(gridCell);
Broadcast<int[]> bandsBroadcast = sc.broadcast(bands);
Broadcast<byte[]> noDataBroadcast = sc.broadcast(nodata);
接下来，是裁切的主体逻辑部分。简单来讲就是先对输入数据进行索引，生成key为格网号，value为原始数据的PairRDD，其中一个输入可以和多个格网相交，因此可以对应多个格网，所以这里采用的算子是flatMapToPair。然后，将生成的PairRDD按照格网号聚集（groupByKey），并在foreach算子内对每一个格网内的数据进行处理。对于每一个格网，先生成一个统一的data空数组，遍历每一个与之相交的数据，裁剪出相交的部分，并将该部分数据写进data中。然后利用gdal将data数组写到本地文件中即可。如下所示：
裁切过程
然后是代码，比较长。
pathsRDD.flatMapToPair(path -> {
gdal.AllRegister(); // 也要先注册驱动
Dataset dataset = gdal.Open(path, gdalconst.GA_ReadOnly);
double[] geoTransform = dataset.GetGeoTransform();
// 计算与该影像相交的格网
double minX = geoTransform[0], maxX, minY, maxY = geoTransform[3];
maxX = minX + dataset.getRasterXSize() * geoTransform[1] + dataset.getRasterYSize() * geoTransform[2];
minY = maxY + dataset.getRasterXSize() * geoTransform[4] + dataset.getRasterYSize() * geoTransform[5];
int minXIndex = gridCellBroadcast.value().getXIndex(minX);
int maxXIndex = gridCellBroadcast.value().getXIndex(maxX);
int minYIndex = gridCellBroadcast.value().getYIndex(maxY);
int maxYIndex = gridCellBroadcast.value().getYIndex(minY);
List<Tuple2<Tuple2<Integer, Integer>, String>> grids = new ArrayList<>();
for (int x = minXIndex; x <= maxXIndex; x++) {
for (int y = minYIndex; y <= maxYIndex; y++) {
grids.add(new Tuple2<>(new Tuple2<>(x, y), path));
}
}
dataset.delete();
return grids.iterator();
}).groupByKey().foreach(grids -> {
gdal.AllRegister();
DatasetInfo curInfo = datasetInfoBroadcast.value();
Tuple2<Integer, Integer> grid = grids._1();
Iterator<String> datasets = grids._2().iterator();
// 创建临时输出文件夹（也可指定一个路径并广播）
Path path = Files.createTempDirectory("gdalTemp");
String datasetPath = Paths.get(path.toString(), grid._2() + "_" + grid._1() + ".tif").toString();
Driver driver = gdal.GetDriverByName("GTiff");
// 获取当前格网的范围
Extent gridExtent = gridCellBroadcast.value().getCellExtent(grid._1(), grid._2());
int gridWidth = (int) ((gridExtent.maxX - gridExtent.minX) / curInfo.xPixel);
int gridHeight = (int) ((gridExtent.minY - gridExtent.maxY) / curInfo.yPixel);
// 创建结果Tiff文件，设置空间信息，此时文件内无数据。
Dataset gridDs = driver.Create(datasetPath, gridWidth, gridHeight, curInfo.bandCount, curInfo.dataType);
gridDs.SetProjection(curInfo.projection);
gridDs.SetGeoTransform(new double[]{gridExtent.minX, curInfo.xPixel, curInfo.xRotate, gridExtent.maxY,
curInfo.yRotate, curInfo.yPixel});
// 创建buffer数组，即上文中的data，并用NoData值填充
int byteToType = gdal.GetDataTypeSize(curInfo.dataType) / byteSizeBroadcast.value();
byte[] dst = new byte[gridWidth * gridHeight * curInfo.bandCount * byteToType];
fillWithNodata(dst, noDataBroadcast.value());
//遍历相交的原始Tiff
while (datasets.hasNext()) {
Dataset dataset = gdal.Open(datasets.next(), gdalconst.GA_ReadOnly);
double[] geoTransform = dataset.GetGeoTransform();
double minX = geoTransform[0];
double maxX = minX + dataset.getRasterXSize() * geoTransform[1] + dataset.getRasterYSize() * geoTransform[2];
double maxY = geoTransform[3];
double minY = maxY + dataset.getRasterXSize() * geoTransform[4] + dataset.getRasterYSize() * geoTransform[5];
// 计算相交范围
Extent overlap = gridExtent.overlap(new Extent(minX, minY, maxX, maxY));
int overlapWidth = (int) ((overlap.maxX - overlap.minX) / curInfo.xPixel);
int overlapHeight = (int) ((overlap.minY - overlap.maxY) / curInfo.yPixel);
// r读取相交部分的数据
int srcXOff = (int) ((overlap.minX - minX) / curInfo.xPixel);
int srcYOff = (int) ((overlap.maxY - maxY) / curInfo.yPixel);
if (srcXOff + overlapWidth > dataset.getRasterXSize()) {
overlapWidth = dataset.getRasterXSize() - srcXOff;
}
if (srcYOff + overlapHeight > dataset.getRasterYSize()) {
overlapHeight = dataset.getRasterYSize() - srcYOff;
}
byte[] data = new byte[overlapWidth * overlapHeight * curInfo.bandCount * byteToType];
dataset.ReadRaster(srcXOff, srcYOff, overlapWidth, overlapHeight, overlapWidth, overlapHeight,
curInfo.dataType, data, bandsBroadcast.value());
// 相交部分的数据写入buffer中
int dstXOff = (int) ((overlap.minX - gridExtent.minX) / curInfo.xPixel);
int dstYOff = (int) ((overlap.maxY - gridExtent.maxY) / curInfo.yPixel);
dataMerge(data, overlapWidth, overlapHeight, dst, gridWidth, gridHeight, dstXOff, dstYOff, byteToType,
noDataBroadcast.value());
}
// 写入本地文件
gridDs.WriteRaster(0, 0, gridWidth, gridHeight, gridWidth, gridHeight,
curInfo.dataType, dst, bandsBroadcast.value());
gridDs.FlushCache();
gridDs.delete();
});
涉及的辅助函数
// 数组赋值NoData（按波段）
private void fillWithNodata(byte[] arr, byte[] noData) {
int count = arr.length / noData.length;
for (int i = 0; i < count; i++) {
System.arraycopy(noData, 0, arr, noData.length * i, noData.length);
}
}
// 原Tiff的buffer数据写入新Tiff的Buffer中，涉及坐标变换
private void dataMerge(byte[] src, int srcXSize, int srcYSize, byte[] dst, int dstXSize, int dstYSize,
int xoff, int yoff, int byteToType, byte[] noData) {
int srcRowPixels = srcXSize * byteToType;
int srcBandPixels = srcRowPixels * srcYSize;
int dstRowPixels = dstXSize * byteToType;
int dstBandPixels = dstRowPixels * dstYSize;
int bandCount = src.length / srcBandPixels;
for (int i = 0; i < bandCount; i++) {
for (int r = 0; r < srcYSize; r++) {
int dstR = r + yoff;
if (dstR >= dstYSize) {
continue;
}
for (int c = 0; c < srcXSize; c++) {
int dstC = c + xoff;
if (dstC >= dstXSize) {
continue;
}
int srcFirstIndex = srcBandPixels * i + srcRowPixels * r + byteToType * c;
int dstFirstIndex = dstBandPixels * i + dstRowPixels * dstR + byteToType * dstC;
// 同一像元的新旧值不同，若其中之一为NoData，设置为另一个的值；否则，取平均
if (ArrayUtils.
isEquals(noData, ArrayUtils.subarray(src, srcFirstIndex, srcFirstIndex + byteToType))) {
continue;
} else if (ArrayUtils.
isEquals(noData, ArrayUtils.subarray(dst, dstFirstIndex, dstFirstIndex + byteToType))) {
for (int b = 0; b < byteToType; b++) {
dst[dstFirstIndex + b] = src[srcFirstIndex + b];
}
} else {
for (int b = 0; b < byteToType; b++) {
dst[dstFirstIndex + b] = (byte) ((dst[dstFirstIndex + b] + src[srcFirstIndex + b]) / 2);
}
}
}
}
}
}
- 一些注意事项
1）GDAL中的对象都没有实现Serializable接口，不能在driver和Executor、各Executor之间传递，所以传递时使用了原始数据路径而不是Dataset对象。此外，也尽量在同一个task中完成Dataset的创建和写入。
2）文中的merge方法在取平均时逻辑并不严密，假如同一像元先后要写入三个值：a、b、c，得到的结果是(((a + b) / 2 ) + c ) / 2。要得到(a + b + c) / 3的结果，需要另外声明一个与buff同样长度的count数组来记录个数。
3）如前文所说，这个并行逻辑可能不是最佳的，应该会有更高效的写法。欢迎交流。
五.影像重投影
gdal提供了重投影的接口，且有多个重载：
gdal重投影接口
我自己一般使用
ReprojectImage(Dataset src_ds, Dataset dst_ds, java.lang.String
src_wkt, java.lang.String dst_wkt, int resampleAlg)
在使用该接口之前，需要自行构建dst_ds即目标数据集，其波段数、数据类型可以与原数据src_ds保持一致，投影与目标投影dst_wkt相同，但其宽、高、空间范围等都需要重新计算。
首先计算空间范围，方法是先单独对原始影像的四个角点进行重投影，得到新影像的角点值。
// 初始化投影工具类CoordinateTransformation
SpatialReference src = new SpatialReference();
src.ImportFromWkt(srcPrj);
SpatialReference dst = new SpatialReference();
dst.ImportFromWkt(dstPrj);
CoordinateTransformation coordinateTransformation = osr.CreateCoordinateTransformation(src, dst);
// 计算角点坐标
double[][] coords = new double{/*利用geoTransform计算四个角的坐标值*/}
// 批量点投影
coordinateTransformation.TransformPoints(coords);
// 得到新的坐标范围
double minX = Double.MAX_VALUE, minY = Double.MAX_VALUE, maxX = Double.MIN_VALUE, maxY = Double.MIN_VALUE;
for (int i = 0; i < coords.length; i++) {
if (coords[i][0] < minX) minX = coords[i][0];
if (coords[i][0] > maxX) maxX = coords[i][0];
if (coords[i][1] < minY) minY = coords[i][1];
if (coords[i][1] > maxY) maxY = coords[i][1];
}
遗憾的是，目标影像的宽高需要自己根据实际情况指定，也可以指定一个分辨率，并根据坐标范围来计算。这样一来，就可以通过输出路径、影像宽高、波段数、数据类型来创建一个新的Dataset，利用上面计算出来的左上角坐标和自己指定的分辨率来设置该Dataset的geoTransform（y方向的分辨率是x方向的负数，rotate一般为0），并设置好参考系即可。
最后是投影过程中采用的重采样方法resampleAlg，重采样的含义及常用方法可自行百度，在gdalconst中提供了如下方法：
gdal中的重采样方法
其中， GRA_NearestNeighbour（最邻近法）是重投影时的默认方法。所有参数都创建好后，直接调用即可：
gdal.ReprojectImage(srcDs, dstDs, srcPrj, dstPrj, reasampleAlg); dstDs.FlushCache();
咕咕咕