遥感影像批量拼接

遥感图像处理

第一章：遥感图形拼接

---

---

前言

GDAL处理遥感影像最大的优势是可以地理信息数据的读取、编码操作。读取遥感影像主要是要获取影像中的
1）数据主要指数组；
2）地理信息中的坐标、投影信息；
本文根据笔者最近写的一个遥感影像批量处理代码中的图像拼接代码举例；文中仅采使用了两张影响拼接进行举例，该方法可以进一步进行扩展。
在实现该功能过程中要注意的点是：
1）数据读取后转换为栅格数据行列的变化：
根据GDAL读取后获取的影像信息可以得到影像的（通道，行：RasterYSize，列：RasterXSize）
但但要注意将数据转换为数组后数据的行，列数；
2）当数据裁剪后很多信息可能已经丢失了，因此可以在裁剪影像时在文件名中记录相关的地理编码信息如：

• 原始图像的大小 坐标系参数（主要是图像左上角x,y坐标信息，其他信息可以通过读取裁剪后的影像获得；
• 裁剪图像时计算的索引信息，要注意对边缘的单独进行处理（边缘裁剪的图像可能不是固定大小）

---

一、拼接步骤

1.数据读取

代码如下：

from osgeo import gdal
import numpy
image1 = gdal.Open(r"1.tif")
image2 = gdal.Open(r"2.tif")
# 获取第一张影像参数
geotransform1 = image1.GetGeoTransform()
projection1 = image1.GetProjection()
# 列
cols1 = image1.RasterXSize
# 行
rows1 = image1.RasterYSize
# 获取第二张影像参数
geotransform2 = image2.GetGeoTransform()
projection2 = image2.GetProjection()
cols2 = image2.RasterXSize
rows2 = image2.RasterYSiz

2.读取数据

# 数据
band1 = image1.GetRasterBand(1)
band2 = image2.GetRasterBand(1)

data1 = band1.ReadAsArray()
data2 = band2.ReadAsArray()
#cows1,cols1 = data.shape
# 创建数组用于存储数据
output_band = output_image.GetRasterBand(1)
output_data = numpy.zeros((new_rows, new_cols), dtype=numpy.float32)

output_data[0:rows1, 0:cols1] = data1
output_data[0:rows2, cols1:new_cols] = data2

3.保存数据

# 写入数据
output_band.WriteArray(output_data)

output_band.FlushCache()
output_image = None
image1 = None
image2 = None

---

二、总结

这是数据拼接的一个思路，肯定还有其他思路，大家遇到了可以一起交流一下

下面对影像的裁剪与拼接进行一下总结：

裁剪

1.先对整张图像划分裁剪的网格，根据裁剪图像的大小；
2.计算每个图像对应的位置，通过索引提取对应区域的图像
3.根据位置索引x,y和坐标转换信息计算裁剪区域的坐标转换信息
4.保存图像，这里我们要在文件名中保存原始图像的大小，坐标系参数（主要是图像左上角x,y坐标信息，其他信息可以通过读取裁剪后的影像获得；
5.用于裁剪图像计算的索引信息，这需要对边缘的单独进行处理（边缘裁剪的图像可能不是固定大小）；

• 倘若我们存储时进行了填充，统一了裁剪后图像的大小。主要解决思路：利用裁剪索引值，用过对原始图像和裁剪后图像索引对应的局部区域并赋值

拼接操作

1）读取图像，获取数据与坐标转换信息
2）创建数组、利用数组索引将图像放回原来的初始位置
3）写入数据，注意写入时数据行列
图像裁剪的代码后面会陆续更新