目录
任务描述
任务分解
GDAL读写tif文件函数介绍及注意点
具体步骤
①影像来源
②确定影像范围
③写入输入影像数据到输出影像中
参考文献
任务描述
将多张有重叠部分的影像经过影像拼接，最终得到一张完整的影像，如下：
➕
═
任务分解
①得到两张有重叠部分的影像；
②确定输出影像的范围（长x宽）；
③以此读取输入影像到输出影像中，保证输入影像在输出影像中的正确x、y偏移值（重叠部分默认用后一张影像值替代；无数据部分默认用黑色填充）。
GDAL读写tif文件函数介绍及注意点
读影像：Dataset.ReadRaster(int xOff,//参数1：读取数据的左上角x坐标；
int yOff, //参数2：读取数据的左上角y坐标；
int xSize, //参数3：输入影像的长（栅格列数）；
int ySize, //参数4：输入影像的宽（栅格行数）；
byte[] buffer, //参数5：存储读取到的输入影像数据的缓存空间，计算得到，具体为缓存长x缓存 宽x波段数；
int buf_xSize, //参数6：缓存空间长（栅格列数）；
int buf_ySize, //参数7：缓存空间高（栅格行数）；
int bandCount, //参数8：要获取的波段数，例如RGB为3；
int[] bandMap, //参数9：波段读取次数数组，{1,2,3}表示按照RGB的顺序读取；
int pixelSpace, //参数10：相邻像素间隔，单位为字节。例如RGB RGB，第一个像素R与第二个 像素R相隔GB，加上自身的一个R，所以间隔为3xbyte(type)；
int lineSpace, //参数11：相邻两行的间隔，单位为字节。如果是从左到右，从上到下的顺序，那 么上一行与下一行的间隔为3xbyte(type)xbuf_xSize；
int bandSpace)//参数12：相邻两波段的间隔，单位为字节。例如R与G两个波段的举例就是对应 buffer类型的字节数，即byte(type)
写影像：Dataset.WriteRaster(int xOff, int yOff, int xSize, int ySize, byte[] buffer, int buf_xSize, int buf_ySize, int bandCount, int[] bandMap, int pixelSpace, int lineSpace, int bandSpace)。
WriteRaster的参数和ReadRaster的参数完全一样，就不再重复介绍了。
本文的影像拼接功能需要先后使用上述两个函数，ReadRaster读取的数据会返回给buffer，因此，WriteRaster中的buffer一定和ReadRaster中的buffer是同一个数组。
需要注意的是，本文中由于输出影像是多张输入影像的拼接图，因此buffer的大小应该以输出影像数据的大小来计算，而不是以输入影像的大小来计算。
另外要注意分清楚参数3、4和参数6、7，前面两个是影像尺寸，后面两个是buffer尺寸。本文中的ReadRaster中这两者是不一样的，但在WriterRaster中是一样的，要特别注意，否则会生成错误的图像。
具体步骤
①影像来源
我在一个软件上下载的。一定要有坐标信息，不然没法判断重叠。
②确定影像范围
大致思路：比较所有输入的影像最大最小x坐标和y坐标，得到所有值中的最大x、最小x、最大y、最小y；
输出影像长（栅格列数）=（最大x-最小x）/ 像素横向分辨率；
输出影像宽（栅格行数）=（最大y-最小y）/ 像素纵向分辨率。
1）首先创建读取影像和输出影像的信息存储结构体
public struct OutTiffInfo//输出影像
{
public DataType dataType;//数据类型
public int xSize;//列
public int ySize;//行
public int bandCount;//波段数
public int depth;//图像深度
public string projection;//投影信息
public double maxX;
public double minX;
public double maxY;
public double minY;
public int[] bandMap;//波段读取顺序
public double[] geoTransform;//影像坐标信息
};
public struct InTiffInfo//输入影像
{
public int xSize;//列
public int ySize;//行
public int bdCount;//波段数
public int depth;//图像深度
public double[] geoTransform;//影像坐标信息
public int[] bandMap;//波段读取顺序
}
public struct OffXY//输入影像在输出影像中的坐标偏移
{
public int offX;
public int offY;
}2）上述的信息需要从打开的tiff影像中读取，主要用到GDAL中的Dataset类。这里需要特别注意的是获取的GeoTransform信息，也就是影像坐标转换信息，它返回一个长度为6的double型数组。6个值分别代表的含义如下：
[0][3]的值表示：影像左上角x坐标和左上角y坐标（tiff影像坐标原点在左下角，因此左上角y为最大y值，左上角x为最小x值）；
[1][5]的值表示：影像像素横向分辨率和纵向分辨率（tiff影像像素纵向分辨率为负值）；
[2][4]表示：影像旋转系数（一般影像为0）。
下面是获取影像信息的代码：
/// <summary>
/// 函数：读取影像
/// </summary>
/// <param name="fullPath">影像打开完整路径</param>
/// <param name="inTiffInfo">存储影像信息</param>
/// <param name="ds1">读取到的影像数据集</param>
/// <param name="i">影像读取顺序标识</param>
private void GdalReadTiff(string fullPath, out InTiffInfo inTiffInfo, out Dataset ds1, int i)
{
inTiffInfo = new InTiffInfo();
//读入第一幅图
string path = @fullPath;
ds1 = Gdal.Open(path, Access.GA_ReadOnly);
Band bd1 = ds1.GetRasterBand(1);//第一个波段
DataType dt1 = bd1.DataType;//数据类型
int cols1 = ds1.RasterXSize;//数据宽度
int rows1 = ds1.RasterYSize;//数据高度
int bdCount1 = ds1.RasterCount;//数据波段数
int depth = Gdal.GetDataTypeSize(dt1) / 8;//图像深度
double[] padTransform = new double[6];//获取影像坐标信息
ds1.GetGeoTransform(padTransform);
inTiffInfo.geoTransform = padTransform;
inTiffInfo.xSize = cols1;
inTiffInfo.ySize = rows1;
inTiffInfo.bdCount = bdCount1;
inTiffInfo.depth = depth;
if(i == 0)//第一幅图像
{
outTiffInfo.depth = depth;
outTiffInfo.dataType = dt1;
outTiffInfo.geoTransform = padTransform;
outTiffInfo.projection = ds1.GetProjection();
}
//获取图像角点坐标
double minX1 = padTransform[0];//左上角x
double maxY1 = padTransform[3];//左上角y
double pixelWidth1 = padTransform[1];//像素宽
double pixelHeight1 = padTransform[5];//像素高(负值)
double maxX1 = minX1 + (cols1 * pixelWidth1);//右下角x
double minY1 = maxY1 + (rows1 * pixelHeight1);//右下加y
StoreExtremeXYBandCount(maxX1, minX1, maxY1, minY1, bdCount1);//存储影像两角坐标与波段数
int[] bandMap = new int[bdCount1];//读取波段的次序；
for (int j = 0; j < bandMap.Length; j++ )
{
bandMap[j] = j + 1;
}
inTiffInfo.bandMap = bandMap;
}
//函数：存储影像最大与最小x、y和波段数
private void StoreExtremeXYBandCount(double maxX, double minX, double maxY, double minY, int bdCount)
{
xMax.Add(maxX);
xMin.Add(minX);
yMax.Add(maxY);
yMin.Add(minY);
bdCountMax.Add(bdCount);
}3）根据读取到的所有输入影像信息，确定输出影像范围和信息，代码如下：
/// <summary>
/// 函数：获取输出图像
/// </summary>
/// <param name="outTiffInfo">输出影像信息结构体</param>
/// <returns></returns>
private Dataset OutImg(ref OutTiffInfo outTiffInfo)
{
xMax.Sort();//所有输入影像最大x坐标升序排序
xMin.Sort();//所有输入影像最小x坐标升序排序
yMax.Sort();//所有输入影像最大y坐标升序排序
yMin.Sort();//所有输入影像最小y坐标升序排序
bdCountMax.Sort();//所有输入影像波段数升序排序
outTiffInfo.minX = xMin[0];
outTiffInfo.maxX = xMax[xMax.Count - 1];
outTiffInfo.minY = yMin[0];
outTiffInfo.maxY = yMax[yMax.Count - 1];
outTiffInfo.bandCount = bdCountMax[bdCountMax.Count - 1];
outTiffInfo.xSize = (int)Math.Ceiling((outTiffInfo.maxX - outTiffInfo.minX) / outTiffInfo.geoTransform[1]);
outTiffInfo.ySize = (int)Math.Ceiling((outTiffInfo.maxY - outTiffInfo.minY) / Math.Abs(outTiffInfo.geoTransform[5]));
outTiffInfo.bandMap = new int[outTiffInfo.bandCount];
for (int i = 0; i < outTiffInfo.bandCount; i++ )
{
outTiffInfo.bandMap[i] = i + 1;
}
//创建图像驱动
Driver dvOut = Gdal.GetDriverByName("GTiff");
string path = @pFileName;
//创建影像
Dataset dtOut = dvOut.Create(path, outTiffInfo.xSize, outTiffInfo.ySize, outTiffInfo.bandCount, outTiffInfo.dataType, null);
//设置影像坐标转换信息和投影信息
double[] geoTransform = new double[] { outTiffInfo.minX, outTiffInfo.geoTransform[1], outTiffInfo.geoTransform[2], outTiffInfo.maxY, outTiffInfo.geoTransform[4], outTiffInfo.geoTransform[5] };
dtOut.SetGeoTransform(geoTransform);
dtOut.SetProjection(outTiffInfo.projection);
MessageBox.Show("影像创建成功");
return dtOut;
}③写入输入影像数据到输出影像中
首先计算输入影像坐标偏移，然后通过tiff读写函数进行数据读写。要注意buffer的申请一定要是输出影像的大小哦。代码如下：
//函数：计算影像在输出图像中的偏移值
private void computeOff(InTiffInfo inTiffInfo, ref OutTiffInfo outTiffInfo, out OffXY offXY)
{
double minX1 = inTiffInfo.geoTransform[0];
double maxY1 = inTiffInfo.geoTransform[3];
double minX = outTiffInfo.minX;
double maxY = outTiffInfo.maxY;
double pixelWidth1 = outTiffInfo.geoTransform[1];
double pixelHeight1 = outTiffInfo.geoTransform[5];
//# 计算图1左上角的偏移值（在输出图像中）
int xOffset1 = (int)((minX1 - minX) / pixelWidth1);
int yOffset1 = (int)((maxY1 - maxY) / pixelHeight1);
offXY.offX = xOffset1;
offXY.offY = yOffset1;
}
/// <summary>
/// 函数：写入输出影像
/// </summary>
/// <param name="dtOut">输出影像数据集</param>
/// <param name="ds1">输入影像数据集</param>
/// <param name="inTiffInfo">输入影像信息结构体</param>
/// <param name="offXY">输入影像在输出影像中的左上角偏移值</param>
/// <param name="outTiffInfo"输出影像信息结构体></param>
private void GdalWriteTiff(ref Dataset dtOut, Dataset ds1, InTiffInfo inTiffInfo, OffXY offXY, OutTiffInfo outTiffInfo)
{
byte[] buf = new byte[outTiffInfo.xSize * outTiffInfo.ySize * outTiffInfo.bandCount];
ds1.ReadRaster(0, 0, inTiffInfo.xSize, inTiffInfo.ySize, buf, outTiffInfo.xSize, outTiffInfo.ySize, inTiffInfo.bdCount, inTiffInfo.bandMap, outTiffInfo.bandCount * outTiffInfo.depth, outTiffInfo.xSize * outTiffInfo.bandCount * outTiffInfo.depth, outTiffInfo.depth);//读取影像数据
dtOut.WriteRaster(offXY.offX, offXY.offY, inTiffInfo.xSize, inTiffInfo.ySize, buf, outTiffInfo.xSize, outTiffInfo.ySize, outTiffInfo.bandCount, outTiffInfo.bandMap, outTiffInfo.bandCount * outTiffInfo.depth, outTiffInfo.xSize * outTiffInfo.bandCount * outTiffInfo.depth, outTiffInfo.depth);
dtOut.FlushCache();//刷新磁盘
}大功告成，在使用GDAL的时候，在调用函数里记得首先要注册GDAL，写上下面一行代码就行了：
Gdal.AllRegister();
