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：创建⼀个输出图像
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]python正则表达式判断
pixelWidth1 = transform1[1]
pixelHeight1 = transform1[5]#是负值（important）
maxX1 = minX1 + (cols1 * pixelWidth1)
minY1 = maxY1 + (rows1 * pixelHeight1)
# 读⼊第⼆幅图像
ds2 = gdal.Open('**.img')
band2 = ds2.GetRasterBand(1)
rows2 = ds2.RasterYSize
cols2 = ds2.RasterXSize
# 获取图像⾓点坐标
transform2 = ds2.GetGeoTransform()
minX2 = transform2[0]
maxY2 = transform2[3]
pixelWidth2 = transform2[1]
pixelHeight2 = transform2[5]
maxX2 = minX2 + (cols2 * pixelWidth2)
minY2 = maxY2 + (rows2 * pixelHeight2)
# 获取输出图像坐标
minX = min(minX1, minX2)
maxX = max(maxX1, maxX2)
minY = min(minY1, minY2)
maxY = max(maxY1, maxY2)
#获取输出图像的⾏与列
cols = int((maxX - minX) / pixelWidth1)
rows = int((maxY - minY) / abs(pixelHeight1))
# 计算图1左上⾓的偏移值（在输出图像中）
xOffset1 = int((minX1 - minX) / pixelWidth1)
yOffset1 = int((maxY1 - maxY) / pixelHeight1)
# 计算图2左上⾓的偏移值（在输出图像中）
xOffset2 = int((minX2 - minX) / pixelWidth1)
yOffset2 = int((maxY2 - maxY) / pixelHeight1)
# 创建⼀个输出图像
driver = ds1.GetDriver()
dsOut = driver.Create('mosiac.img', cols, rows, 1, band1.DataType)#1是bands，默认
bandOut = dsOut.GetRasterBand(1)
# 读图1的数据并将其写到输出图像中
data1 = band1.ReadAsArray(0, 0, cols1, rows1)
bandOut.WriteArray(data1, xOffset1, yOffset1)
#读图2的数据并将其写到输出图像中
data2 = band2.ReadAsArray(0, 0, cols2, rows2)
bandOut.WriteArray(data2, xOffset2, yOffset2)
''' 写图像步骤'''
# 统计数据
bandOut.FlushCache()#刷新磁盘
stats = bandOut.GetStatistics(0, 1)#第⼀个参数是1的话，是基于⾦字塔统计，第⼆个
#第⼆个参数是1的话：整幅图像重度，不需要统计
# 设置输出图像的⼏何信息和投影信息
geotransform = [minX, pixelWidth1, 0, maxY, 0, pixelHeight1]
dsOut.SetGeoTransform(geotransform)
dsOut.SetProjection(ds1.GetProjection())
# 建⽴输出图像的⾦字塔
gdal.SetConfigOption('HFA_USE_RRD', 'YES')
dsOut.BuildOverviews(overviewlist=[2,4,8,16])#4层
补充知识：运⽤Python的第三⽅库：GDAL进⾏遥感数据的读写
0 背景及配置环境
0.1 背景
GDAL(Geospatial Data Abstraction Library)是⼀个在X/MIT许可协议下的开源栅格空间数据转换库。它利⽤抽象数据模型来表达所⽀持的各种⽂件格式。它还有⼀系列命令⾏⼯具来进⾏数据转换和处理。
这个开源栅格空间数据转换库拥有许多和其他语⾔的接⼝，对于python，他有对应的第三⽅包GDAL，下载安装已在上篇⽂章中提到。
⽬的: 可以使⽤Python的第三⽅包：GDAL进⾏遥感数据的读写，⽅便批处理。
0.2 配置环境
电脑系统： win7x64
Python版本: 3.6.4
GDAL版本: 2.3.2
1 读
1.1 TIFF格式
标签图像⽂件格式（Tag Image File Format，简写为TIFF）是⼀种灵活的位图格式，主要⽤来存储包括照⽚和艺术图在内的图像。它最初由Aldus公司与微软公司⼀起为PostScript打印开发。TIFF与JPEG和PNG⼀起成为流⾏的⾼位彩⾊图像格式。
TIFF⽂件以.tif为扩展名。
def tif_read(tifpath, bandnum):
"""
Use GDAL to read data and transform them into arrays.
:param tifpath:tif⽂件的路径
:param bandnum：需要读取的波段
:return：该波段的数据，narray格式。len(narray)是⾏数，len(narray[0])列数
"""
image = gdal.Open(tifpath) # 打开该图像
if image == None:
print(tifpath + "该tif不能打开!")
return
lie = image.RasterXSize # 栅格矩阵的列数
hang = image.RasterYSize # 栅格矩阵的⾏数
im_bands = image.RasterCount # 波段数
im_proj = image.GetProjection() # 获取投影信息
im_geotrans = image.GetGeoTransform() # 仿射矩阵
print('该tif:{}个⾏，{}个列，{}层波段, 取出第{}层.'.format(hang, lie, im_bands, bandnum))
band = image.GetRasterBand(bandnum) # Get the information of band num.
band_array = band.ReadAsArray(0,0,lie,hang) # Getting data from zeroth rows and 0 columns
# band_df = pd.DataFrame(band_array)
del image # 减少冗余
return band_array, im_proj, im_geotrans
2 写
2.1 TIFF格式
TIFF格式的数据格式有：Byete、int16、uint16、int32、uint32、float32、float64等7余种。
⾸先，要判断数据的格式，才能按需求写出。
def tif_write(self, filename, im_data, im_proj, im_geotrans):
"""
gdal数据类型包括
gdal.GDT_Byte,
gdal.GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,
gdal.GDT_Float32, gdal.GDT_Float64
:
param filename: 存出⽂件名
:param im_data: 输⼊数据
:param im_proj: 投影信息
:param im_geotrans: 放射变换信息
:return: 0
"""
if 'int8' in im_data.dtype.name: # 判断栅格数据的数据类型
datatype = gdal.GDT_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
# 判读数组维数
if len(im_data.shape) == 3:
im_bands, im_height, im_width = im_data.shape
else:
im_bands, (im_height, im_width) = 1,im_data.shape # 多维或1.2维
#创建⽂件
driver = gdal.GetDriverByName("GTiff")  #数据类型必须有，因为要计算需要多⼤内存空间
dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
dataset.SetGeoTransform(im_geotrans)    #写⼊仿射变换参数
dataset.SetProjection(im_proj)    #写⼊投影
if im_bands == 1:
dataset.GetRasterBand(1).WriteArray(im_data) #写⼊数组数据
else:
for i in range(im_bands):
dataset.GetRasterBand(i+1).WriteArray(im_data[i])
del dataset
3 展⽰
3.1 TIFF格式
# 这个展⽰的效果并不是太好，当做⽰意图⽤
def tif_display(self,im_data):
"""
:
param im_data: 影像数据,narray
:return: 展出影像
"""
# plt.imshow(im_data,'gray') # 必须规定为显⽰的为什么图像
plt.imshow(im_data) # 必须规定为显⽰的为什么图像
plt.show() # 显⽰出来，不要也可以，但是⼀般都要了
以上这篇python+gdal+遥感图像拼接(mosaic)的实例就是⼩编分享给⼤家的全部内容了，希望能给⼤家⼀个参考，也希望⼤家多多⽀持。