这篇文章主要介绍了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