如何用Python| 製作遙感影像拼接

阿新 • • 發佈：2020-11-07

本文的文字及圖片來源於網路,僅供學習、交流使用,不具有任何商業用途,版權歸原作者所有,如有問題請及時聯絡我們以作處理

以下文章來源於騰訊雲，作者：bugsuse。

0.前言

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

在上一篇推文中，我展示瞭如何使用Python結合Landsat製作遙感影像圖（Python乾貨 | 製作遙感影像圖）。

對於Landsat資料來說，對某個區域的重訪週期為16天，每個位置使用全球參考系（WRS）進行索引，即每一個位置都會對應一個Path和Row，相鄰的影像之間會有部分割槽域是重疊的。

Fig.1 World Reference System

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

單張影像是這樣。

本文合併後是這樣。

1.準備工作

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

import rasterio as rio
import 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, prefix
def 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，如果是其他遙感資料，需對應進行修改。

```python
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, top
        out_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 * 255
            data[np.where(data == 255)] = 0
            # 影像重合部分處理，重合部分取最大值
            xSize = fnBand.XSize
            ySize = fnBand.YSize
            outData = out_band.ReadAsArray(x, y, xSize, ySize)
            data = np.maximum(data, outData)
            out_band.WriteArray(data, x, y)
        del out_band, out_ds
        file2read = ibands + 'mosaic.tif'
        arr[i, :, :] = tiff.imread(file2read)
        os.remove(file2read)
    plot_rgb(arr, rgb=(0, 1, 2))