GDAL切割影像数据
正在学习GDAL,发现网上的资料,java的资料很少,大部分是python与C++进行的数据处理。下面是我在学习过程中,使用GDAL编写的第一个处理影像的测试程序,使用Band与Dataset两种接口方式切割影像数据:
package GDALExample.TIFRead;
import org.gdal.gdal.Band;
import org.gdal.gdal.Dataset;
import org.gdal.gdal.Driver;
import org.gdal.gdal.gdal;
import org.gdal.gdalconst.gdalconstConstants;
/**
* 影像数据切割
*
*/
public class TifClip {
static int bandNum = -1;
static String fileName_tif = "D:\\DATA\\p1dom28.tif";
static int rCount = -1;
static int dataType = -1;
public static void main(String[] args) {
if (gdal.GetDriverCount() < 1)
gdal.AllRegister();
// int iCount = gdal.GetDriverCount();
String outPath = "D:\\GDAL\\clip";
Dataset dsTIF = gdal.Open(fileName_tif, gdalconstConstants.GA_ReadOnly);
// 切割为 1000*1000的影像块
// clipTIF(dsTIF, outPath, 1000, 1000); // 使用Band接口
clipTIFOnDSReadWrite(dsTIF, outPath, 1000, 1000); // 使用Dataset接口
}
/**
* 切割影像 使用Band接口读写影像
*
* @param dsTIF
* @param outPath
* @param clipWidth
* @param clipHeight
*/
public static void clipTIF(Dataset dsTIF, String outPath, int clipWidth, int clipHeight) {
Driver driver = gdal.GetDriverByName("GTIFF");
int pWidth = clipWidth;
int pHeight = clipHeight;
// 波段数
rCount = dsTIF.getRasterCount();
// 从波段中获取影像的数据类型,gdal中波段索引从1开始
dataType = dsTIF.GetRasterBand(1).GetRasterDataType();
// 六参数信息
double[] geoTransform = dsTIF.GetGeoTransform();
// 计算切割为多少列、多少行
int iBlockCols = dsTIF.getRasterXSize() / clipWidth;
int iBlockRows = dsTIF.getRasterYSize() / clipHeight;
if (dsTIF.getRasterXSize() % clipWidth > 0)
iBlockCols++;
if (dsTIF.getRasterYSize() % clipHeight > 0)
iBlockRows++;
int iBlock = 1;
double dCurrYcoor = geoTransform[3];
for (int iCol = 0; iCol < iBlockCols; iCol++) {
// 设置要裁剪的起始像元位置,以及各方向的像元数
geoTransform[3] = dCurrYcoor;// 重置为左上角的坐标(Y坐标),按列输出
if (iCol > 0)
geoTransform[0] = geoTransform[0] + (clipWidth * geoTransform[1] + clipWidth * geoTransform[2]);
for (int iRow = 0; iRow < iBlockRows; iRow++) {
// 计算裁剪后的左上角坐标,设定裁剪后影像块的坐标位置
if (iRow > 0)
geoTransform[3] = geoTransform[3] + (clipHeight * geoTransform[4] + clipHeight * geoTransform[5]);
if (dsTIF.getRasterXSize() % clipWidth > 0 && iCol == iBlockCols - 1) {
clipWidth = dsTIF.getRasterXSize() % pWidth;
System.out.println(clipWidth);
} else {
clipWidth = pWidth;
System.out.println(clipWidth);
}
if (dsTIF.getRasterYSize() % clipHeight > 0 && iRow == iBlockRows - 1) {
clipHeight = dsTIF.getRasterYSize() % clipHeight;
} else {
clipHeight = pHeight;
}
// 创建结果图像
Dataset outputDs = driver.Create(outPath + iBlock++ + ".tif", clipWidth, clipHeight, rCount, dataType);
outputDs.SetGeoTransform(geoTransform);
outputDs.SetProjection(dsTIF.GetProjection());
// 按band读取
for (int i = 0; i < clipHeight; i++) {
// 按行读取
for (int j = 1; j <= rCount; j++) {
Band orgBand = dsTIF.GetRasterBand(j);
int[] cache = new int[clipWidth];
// 从位置x开始,只读一行
orgBand.ReadRaster(iCol * pWidth, iRow * pHeight + i, clipWidth, 1, cache);
Band desBand = outputDs.GetRasterBand(j);
// 从左上角开始,只写一行
desBand.WriteRaster(0, i, clipWidth, 1, cache);
desBand.FlushCache();
}
}
outputDs.delete();
// if (iBlock > 5) break;
}
}
// 释放资源
dsTIF.delete();
}
/**
* 切割影像 使用Dataset接口读写影像
*
* @param dsTIF
* @param outPath
* @param clipWidth
* @param clipHeight
*/
public static void clipTIFOnDSReadWrite(Dataset dsTIF, String outPath, int clipWidth, int clipHeight) {
Driver driver = gdal.GetDriverByName("GTIFF");
int pWidth = clipWidth;
int pHeight = clipHeight;
// 宽、高、波段数
rCount = dsTIF.getRasterCount();
// 从波段中获取影像的数据类型,gdal中波段索引从1开始
dataType = dsTIF.GetRasterBand(1).GetRasterDataType();
// 六参数信息
double[] geoTransform = dsTIF.GetGeoTransform();
int iBlockCols = dsTIF.getRasterXSize() / clipWidth;
int iBlockRows = dsTIF.getRasterYSize() / clipHeight;
if (dsTIF.getRasterXSize() % clipWidth > 0)
iBlockCols++;
if (dsTIF.getRasterYSize() % clipHeight > 0)
iBlockRows++;
int iBlock = 1;
double dCurrYcoor = geoTransform[3];
for (int iCol = 0; iCol < iBlockCols; iCol++) {
geoTransform[3] = dCurrYcoor;// 重置为左上角的坐标(Y坐标),按列输出
if (iCol > 0)
geoTransform[0] = geoTransform[0] + (clipWidth * geoTransform[1] + clipWidth * geoTransform[2]);
for (int iRow = 0; iRow < iBlockRows; iRow++) {
// 计算裁剪后的左上角坐标,设定裁剪后影像块的坐标位置
if (iRow > 0)
geoTransform[3] = geoTransform[3] + (clipHeight * geoTransform[4] + clipHeight * geoTransform[5]);
if (dsTIF.getRasterXSize() % clipWidth > 0 && iCol == iBlockCols - 1) {
clipWidth = dsTIF.getRasterXSize() % pWidth;
// System.out.println(clipWidth);
} else {
clipWidth = pWidth;
// System.out.println(clipWidth);
}
if (dsTIF.getRasterYSize() % clipHeight > 0 && iRow == iBlockRows - 1) {
clipHeight = dsTIF.getRasterYSize() % clipHeight;
} else {
clipHeight = pHeight;
}
// 创建结果图像
Dataset outputDs = driver.Create(outPath + iBlock++ + ".tif", clipWidth, clipHeight, rCount, dataType);
outputDs.SetGeoTransform(geoTransform);
outputDs.SetProjection(dsTIF.GetProjection());
int[] bands = new int[] { 1, 2, 3 };
// 按Dataset读取
for (int i = 0; i < clipHeight; i++) {
byte[] cache = new byte[clipHeight * rCount];
dsTIF.ReadRaster(iCol * pWidth, iRow * pHeight + i, clipWidth, 1, clipHeight, 1, dataType, cache,
bands);
outputDs.WriteRaster(0, i, clipWidth, 1, clipHeight, 1, dataType, cache, bands);
outputDs.FlushCache();
}
outputDs.delete();
}
}
// 释放资源
dsTIF.delete();
}
}