GeoTiff及GDAL切图（java）

一、GeoTiff

GeoTiff是很常见的一种遥感影像数据格式，它是一种栅格数据，这种数据是以像素点矩阵的形式存储的。

TIF文件中的空间信息包括投影信息和仿射变换信息。其中，仿射变换（经纬度等）的信息分别为：左上角的坐标（经度和维度）、东西和南北方向的像素分辨率（也就是每个像素点的经纬度偏移量）以及图像的旋转系数（正北方向时为零）。这样，每个像素点的坐标就可以通过它与左上角的偏移量和每个像素的分辨率值计算出来。

切图产生的新图的投影信息是不变的，但仿射信息会发生变化，因为左上角的原点可能跟原图不一样了。因此，需要重新计算仿射信息并设置（投影信息也需要设置，但不需要重新计算）。

二、波段

波段又称为波谱段或波谱带，在遥感技术中，通常把电磁波谱划分为大大小小的段落，大的成为波段区，如可见区、红外区等；中等的如近红外、远红外等；小的称为波段。

影像数据根据波段的多少可以分为单波段影像和多波段影像两种，单波段影像一般用黑白色的灰度图来描述，多波段常用RGB 合成象素值的彩色图来描述，就是将三个波段的数据分别通过红、绿、蓝三个通道加载，从而渲染出。

将多波段影像数据添加到地图中之后，可使用多波段栅格数据集中的任意三个可用波段的组合来创建 RGB 合成图。与仅处理一个波段相比，通过将多个波段共同显示为RGB 合成图通常可从数据集收集到更多信息。

三、GDAL切图

Tif文件的读取
根据Tif文件中的波段数量，读取相应的波段(Band)。
根据裁剪的起始点和矩形框的大小，并读取相应的数据
获取创建文件并创建新的目标文件
计算新图的仿射信息，并设置投影和仿射信息
将提取出来的数据写入新文件并刷新缓存到磁盘
释放文件句柄资源

源码如下：

import org.gdal.gdal.Dataset;
import org.gdal.gdal.Driver;
import org.gdal.gdal.gdal;
import org.gdal.gdalconst.gdalconstConstants;public class GDALClip {public static void clip(String srcFile, String dstFile, int offsetX, int offsetY, int sizeX, int sizeY) {// 读取要切的原图Dataset srcDs = gdal.Open(srcFile, gdalconstConstants.GA_ReadOnly);if (srcDs == null) {System.err.println("GDALOpen failed - " + gdal.GetLastErrorNo());System.err.println(gdal.GetLastErrorMsg());return;}int rasterCount = srcDs.GetRasterCount();System.out.println("波段数 = " + rasterCount);if (rasterCount < 1) {System.err.println("异常：没有波段。");return;}// 获取Tif的驱动，为创建切出来的图文件做准备Driver GTIFFDiver = gdal.GetDriverByName("GTIFF");// 创建切出来的要存的文件Dataset dstDs = GTIFFDiver.Create(dstFile, sizeX, sizeY, rasterCount, srcDs.GetRasterBand(1).getDataType());if (dstDs == null) {System.err.println("Create failed - " + gdal.GetLastErrorNo());System.err.println(gdal.GetLastErrorMsg());return;}// 复制每个波段的数据for (int index = 1; index <= rasterCount; index++) {copyData(srcDs, dstDs, index, offsetX, offsetY, sizeX, sizeY);}// 仿射变换系数，即6个参数，左上角地理坐标，纵横方向上的分辨率，以及旋转系数。// GeoTransform[0] /* 左上角x坐标  */// GeoTransform[1] /* 东西方向分辨率resolution */// GeoTransform[2] /* 旋转角度, 0表示图像 "北方朝上" */// GeoTransform[3] /* 左上角y坐标 */// GeoTransform[4] /* 旋转角度, 0表示图像 "北方朝上" */// GeoTransform[5] /* 南北方向分辨率resolution */ // 获取原图的原点坐标信息double[] srcTransform = srcDs.GetGeoTransform();System.out.println("原点坐标   = " + srcTransform[0] + "," + srcTransform[3]);System.out.println("像素坐标差 = " + srcTransform[1] + "," + srcTransform[5]);// 根据反射变换参数计算新图的原点坐标double top_left_x = srcTransform[0] + offsetX * srcTransform[1];double top_left_y = srcTransform[3] + offsetY * srcTransform[5];// 将计算后的值组装为一个元组double[] dstTransform = {top_left_x, srcTransform[1], srcTransform[2], top_left_y, srcTransform[4], srcTransform[5]};// 设置裁剪出来图的原点坐标dstDs.SetGeoTransform(dstTransform);// 设置SRS属性（投影信息）// System.out.println("Projection = " + srcDs.GetProjection());dstDs.SetProjection(srcDs.GetProjection());// 将缓存写入磁盘dstDs.FlushCache();dstDs.delete();srcDs.delete();}private static void copyData(Dataset srcDs, Dataset dstDs, int bandNum, int offsetX, int offsetY, int sizeX, int sizeY) {int bufR[] = new int[sizeX * sizeY];//从波段中切需要的矩形框内的数据(注意读取的矩形框不能超过原图大小)srcDs.GetRasterBand(bandNum).ReadRaster(offsetX, offsetY, sizeX, sizeY, bufR);// 写入目标文件dstDs.GetRasterBand(bandNum).WriteRaster(0, 0, sizeX, sizeY, bufR);}public static void main(String[] args) {gdal.AllRegister();GDALClip.clip("D:\\test\\cut_img.tif", "D:\\test\\clip.tif", 0, 0, 100, 100);gdal.GDALDestroyDriverManager();}
}