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）
 

• 1
• 2

step2 ：計算輸出影象的min X ,max X,min Y,max Y

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

• 1
• 2

y座標同理 
step3：計算輸出影象的行與列

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

• 1
• 2

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)
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層