0.前言

因為沒有喝上“秋天的第一份奶茶”，準備來更新一篇推送。

在上一篇推文中，我展示了如何使用Python結合Landsat制作遙感影像圖(Python干貨 | 制作遙感影像圖)。

對于Landsat數據來說，對某個區域的重訪周期為16天，每個位置使用全球參考系(WRS)進行索引，即每一個位置都會對應一個Path和Row，相鄰的影像之間會有部分區域是重疊的。

Fig.1 World Reference System

在某些遙感影像的應用場景中，如果我們關注的區域正好處于兩景影像的交界處，如下圖中的象山港，那我們就需要將影像拼接起來才可以使用。

單張影像是這樣。

本文合并后是這樣。

1.準備工作

相較于上一篇推送，我們這次為了實現遙感影像的鑲嵌拼接，我們使用到了兩個庫，rasterio和gdal。

import rasterio as rioimport gdal

先介紹一下我們實現兩組遙感影像拼接的思路，首先選取兩景相鄰的影像，分別得到他們的空間范圍，再得到兩景組合到一起之后的空間范圍，使用gdal新建一個tif文件(數據中轉用)，分別得到原來兩景影像在新建的tif文件中的起始位置，將對應的數據寫入新的tif文件中，即實現鑲嵌拼接。

上面說的是兩景影像的拼接，如果是更多影像拼接同樣適用，但是現階段的方法如果拼接多的影像的話，需要的內存空間很大，容易導致內存溢出，感興趣的朋友可以思考一下如何高效實現多景影像的拼接。

其中還有兩處關鍵處理，一是如何去除重疊區域的無效信息，二是重疊區域的數據如何選擇。希望各位看官能從代碼里面找到答案。

2.動起手來

得到輸入影像的四個角點。

def?tiffileList2filename(tiffileList):filename = []prefix = []for ifile in tiffileList:file0 = ifile.split("\\")[-1]prefix.append(os.path.join(ifile, file0))filename.append(os.path.join(ifile, file0) + "_B1.TIF")return?filename,?prefixdef get_extent(tiffileList):filename, prefix = tiffileList2filename(tiffileList)rioData = rio.open(filename[0])left = rioData.bounds[0]bottom = rioData.bounds[1]right = rioData.bounds[2]top?=?rioData.bounds[3]for ifile in filename[1:]:rioData = rio.open(ifile)left = min(left, rioData.bounds[0])bottom = min(bottom, rioData.bounds[1])right = max(right, rioData.bounds[2])top?=?max(top,?rioData.bounds[3])return?left,?bottom,?right,?top,?filename,?prefix

得到新建tif文件的size，這里已知Landsat空間分辨率為30m，如果是其他遙感數據，需對應進行修改。

def getRowCol(left, bottom, right, top):cols = int((right - left) / 30.0)rows = int((top - bottom) / 30.0)return cols, rows

主程序，其中plot_rgb為上一篇推送中用到的函數。

if?__name__?==?'__main__':tiffileList = [r'PathofLandsat8\LC08_L1TP_118039_20160126_20170330_01_T1',r'PathofLandsat8\LC08_L1TP_118040_20160126_20170330_01_T1']left, bottom, right, top, filename, prefix = get_extent(tiffileList)cols, rows= getRowCol(left, bottom, right, top)bands = ['B7', 'B5', 'B3']n_bands = len(bands)arr = np.zeros((n_bands, rows, cols), dtype=np.float)# 打開一個tif文件in_ds = gdal.Open(filename[0])for i in range(len(bands)):ibands = bands[i]# 新建一個tif文件driver = gdal.GetDriverByName('gtiff')out_ds = driver.Create(ibands + 'mosaic.tif', cols, rows)# 設置tif文件的投影out_ds.SetProjection(in_ds.GetProjection())out_band = out_ds.GetRasterBand(1)# 設置新tif文件的地理變換gt = list(in_ds.GetGeoTransform())gt[0], gt[3] = left, topout_ds.SetGeoTransform(gt)# 對要拼接的影像進行循環讀取for ifile in prefix:in_ds = gdal.Open(ifile + '_' + ibands + '.tif')# 計算新建的tif文件及本次打開的tif文件之間的坐標漂移trans = gdal.Transformer(in_ds, out_ds, [])# 得到偏移起始點success, xyz = trans.TransformPoint(False, 0, 0)x, y, z = map(int, xyz)# 讀取波段信息fnBand = in_ds.GetRasterBand(1)data = fnBand.ReadAsArray()# 寫入tif文件之前，最大值設置為255，這一步很關鍵data = data / 65535 * 255data[np.where(data == 255)] = 0# 影像重合部分處理，重合部分取最大值xSize = fnBand.XSizeySize = fnBand.YSizeoutData = out_band.ReadAsArray(x, y, xSize, ySize)data = np.maximum(data, outData)out_band.WriteArray(data, x, y)del out_band, out_dsfile2read = ibands + 'mosaic.tif'arr[i, :, :] = tiff.imread(file2read)os.remove(file2read)plot_rgb(arr,?rgb=(0,?1,?2))

3.小結

大功告成！

因遙感數據較大，可后臺回復 Landsat拼接獲取數據及源碼。

數據處理、分析、可視化

行業資訊、學習資料

長按關注不迷路

總結

以上是生活随笔為你收集整理的python 拼接 遥感影像_Python干货 | 遥感影像拼接的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。