python基础教程：python+gdal+遥感图像拼接(mosaic)的实

这篇文章主要介绍了python+gdal+遥感图像拼接(mosaic)的实例，具有很好的参考价值，希望对大家有所帮助。一起跟随小编过来看看吧
作为摄影测量与遥感的从业者，笔者最近开始深入研究gdal，为工作打基础！个人觉得gdal也是没有什么技术含量，调用别人的api。但是想想这也是算法应用的一个技能，多学无害！

关于遥感图像的镶嵌，主要分为6大步骤：

step1：

1）对于每一幅图像，计算其行与列；

2）获取左上角X,Y

3）获取像素宽和像素高

4）计算max X 和 min Y，切记像素高是负值

maxX1 = minX1 + (cols1 * pixelWidth)
minY1 = maxY1 + (rows1 * pixelHeight）

step2 ：计算输出图像的min X ,max X,min Y,max Y

minX = min(minX1, minX2, …)
maxX = max(maxX1, maxX2, …)

y坐标同理

step3：计算输出图像的行与列

cols = int((maxX – minX) / pixelWidth)
rows = int((maxY – minY) / abs(pixelHeight)

step 4：创建一个输出图像

driver.create()

step 5：

1）计算每幅图像左上角坐标在新图像的偏移值

2）依次读入每幅图像的数据并利用1)计算的偏移值将其写入新图像中

step6 ：对于输出图像

1）刷新磁盘并计算统计值

2）设置输出图像的几何和投影信息

3）建立金字塔

下面附上笔者的代码：

#mosica 两张图像
import os, sys, gdal
from gdalconst import *
os.chdir('c:/temp/****')#改变文件夹路径
# 注册gdal(required)
gdal.AllRegister()
 
# 读入第一幅图像
ds1 = gdal.Open('**.img')
band1 = ds1.GetRasterBand(1)
rows1 = ds1.RasterYSize
cols1 = ds1.RasterXSize
 
# 获取图像角点坐标
transform1 = ds1.GetGeoTransform()
minX1 = transform1[0]
maxY1 = transform1[3]
pixelWidth1 = transform1[1]
pixelHeight1 = transform1[5]#是负值（important）
maxX1 = minX1 + (cols1 * pixelWidth1)
minY1 = maxY1 + (rows1 * pixelHeight1)
 
# 读入第二幅图像
ds2 = gdal.Open('**.img')
band2 = ds2.GetRasterBand(1