歡迎關注博主的微信公眾號：“智慧遙感”。

該公眾號將為您奉上Python地學分析、爬蟲、資料分析、Web開發、機器學習、深度學習等熱門原始碼。

本人的GitHub程式碼資料主頁（持續更新中，多給Star，多Fork）：

https://github.com/xbr2017

CSDN也在同步更新：

https://blog.csdn.net/XBR_2014

“ 遙感影像的特點之一就是同一個影象檔案可以儲存若干個波段的影象，本節內容主要介紹GDAL將多個尺寸相同的影象疊加到一起，以GeoTiff格式輸出，這有利於不同波段之間進行數值運算。”

本節以LandSat影像作為案列，來實現多波段疊加功能。美國NASA的陸地衛星(Landsat)計劃從1972年7月23日以來，已發射8顆。目前Landsat1—4均相繼失效，Landsat 5於2013年6月退役。Landsat 7於1999年4月15日發射升空。Landsat8於2013年2月11日發射升空，經過100天測試執行後開始獲取影像。

GDAL實現多波段疊加

現在理論已經不在了，讓我們學習如何使用GDAL處理這些資料集。柵格資料存在許多不同的檔案格式，GDAL是一個非常流行且強大的庫，用於讀寫許多檔案格式。GDAL庫是開源的，但具有許可，因此即使許多商業軟體包也使用它。

GDAL庫以其讀寫較多不同格式的資料而聞名，但它也包含一些資料處理功能，如最鄰近法。在許多情況下，你仍然需要編寫自己的處理程式碼，但這是對於許多型別的分析來說相對容易。有一個名為NumPy的Python模組，用於處理大型資料陣列，你可以使用GDAL直接將資料讀入NumPy陣列。根據需要操作資料後，使用NumPy或其他適用於這些陣列的模組，你可以將陣列作為柵格資料集寫入磁碟，這是一個非常簡單的過程。

為了說明如何使用GDAL讀取和寫入柵格資料，讓我們從一個將三個Landsat波段組合成一個疊加影象的示例開始，如圖1所示。

圖1  紅色（A），綠色（B）和藍色（C）Landsat帶以黑色和白色顯示，你可以看到它們看起來有點不同。D圖顯示這三個波段疊加成RGB影象。

Landsat計劃是美國地質調查局（USGS）與美國國家航空航天局（NASA）之間的聯合倡議，自1972年以來一直在全球範圍內採集中等解析度的衛星影象。Landsat影象由USGS作為集合分發GeoTIFFs，每個資料集均有波段。除了波段6（熱）和8（全色）之外，每個波段都有30米的解析度，因為它們來自同一個Landsat場景，尺寸相同。這使得事情變得簡單，因為波段直接相互疊加而不需要重取樣之類的處理。下面的程式實現瞭如何建立具有相同尺寸的三波段資料集，然後將波段3，2和1複製到其中。這三個波段分別對應於可見光的紅色，綠色和藍色波長，因此將它們按此順序排列將產生RGB（紅色，綠色，藍色）影象，其外觀與你自己的眼睛非常相似。圖1顯示了黑色和白色的單獨紅色，藍色和綠色波段，以及生成的三波段自然彩色影象。如果你在ArcGIS中預覽影象，它們可能看起來與此類似，因為ArcGIS很可能會拉伸它們。下面來看看Python程式碼的實現：

# _*_ coding: utf-8 _*_
__author__ = 'xbr'
__date__ = '2018/12/3 11:57'

import os
from osgeo import gdal

# 當前所在路徑
os.chdir(r'D:\osgeopy-data\Landsat\Washington')
band1_fn = 'p047r027_7t20000730_z10_nn10.tif'
band2_fn = 'p047r027_7t20000730_z10_nn20.tif'
band3_fn = 'p047r027_7t20000730_z10_nn30.tif'

in_ds = gdal.Open(band1_fn)
in_band = in_ds.GetRasterBand(1)

gtiff_driver = gdal.GetDriverByName('GTiff')
out_ds = gtiff_driver.Create('nat_color.tif',
         in_band.XSize, in_band.YSize, 3, in_band.DataType)
out_ds.SetProjection(in_ds.GetProjection())
out_ds.SetGeoTransform(in_ds.GetGeoTransform())

# 讀取第1波段資料
in_data = in_band.ReadAsArray()
out_band = out_ds.GetRasterBand(3)
out_band.WriteArray(in_data)

# 讀取第2波段資料
in_ds = gdal.Open(band2_fn)
out_band = out_ds.GetRasterBand(2)
out_band.WriteArray(in_ds.ReadAsArray())

# 讀取第3波段資料
out_ds.GetRasterBand(1).WriteArray(
    gdal.Open(band3_fn).ReadAsArray())

out_ds.FlushCache()
for i in range(1, 4):
    out_ds.GetRasterBand(i).ComputeStatistics(False)

out_ds.BuildOverviews('average', [2, 4, 8, 16, 32])
del out_ds

程式碼中首先需要匯入gdal模組，然後設定當前目錄並指定哪個檔案對應哪個Landsat波段。然後通過將檔名傳遞給gdal.Open開啟包含第一個波段的GeoTIFF。你還可以獲取資料集中第一個也是唯一一個波段的控制代碼，儘管你尚未讀取任何資料。請注意，使用索引1而不是0來獲取第一個波段。當你使用GetRasterBand時，帶編號始終以1開頭，很多人會忘記這個知識點，容易出錯。無論如何，在建立輸出影象之前，你需要此波段物件，因為它具有輸出檔案所需要的基本資訊（尺寸大小、投影、座標等）。

友情提示  請記住，波段索引從1開始而不是0。

接下來，你將建立一個新資料集以將像元資料複製到其中。你必須使用驅動程式物件來建立新資料集，因此你可以找到GeoTIFF驅動程式，然後使用其Create函式。這是該功能的完整引數：

driver.Create(filename, xsize, ysize, [bands], [data_type], [options])

• filename是要建立的資料集的路徑。
• xsize是新資料集中的列數。
• ysize是新資料集中的行數。
• bands是新資料集中的波段數。預設值為1。
• data_type是將儲存在新資料集中的資料型別。預設值為GDT_Byte。
• options是建立選項字串的列表。可能的值取決於正在建立的資料集的型別。

因為你使用GeoTIFF驅動程式，所以無論你提供什麼副檔名，輸出檔案都將是GeoTIFF。但是，副檔名不會自動新增，因此你需要提供該擴充套件程式。在這種情況下，你將其命名為nat_color.tif並將其儲存在指定資料夾中，因為這是使用os.chdir設定的當前資料夾。你還需要提供列數和行數建立新資料集，因此分別使用XSize和YSize屬性從輸入波段獲取該資訊。Open的下一個引數是band的數量，你希望這個新的raster有三個。下一個可選引數是資料型別，它必須是表1中的值之一。你可以從輸入波段獲取此資訊，但在這種情況下你可能會忽略它，因為這些影象使用預設型別的DT_Byte。

            表1  GDAL資料型別常量

你從輸入資料集中獲取投影（SRS）並將其複製到新資料集，然後對地理轉換執行相同操作。地理轉換很重要，因為它提供原點座標和像元大小。如上所述，在將資料集放置在正確的空間位置時，原點和像元大小非常重要。雖然你不必在新增像元值之前新增投影和地理轉換資訊，但在建立新資料集後立即將其解除。

設定資料集後，就可以新增像元值了。因為你已經擁有Landsat band 1的GeoTIFF中的波段物件，所以你可以將其中的像元值讀入NumPy陣列。如果不向ReadAsArray提供任何引數，則所有像元值都以二維陣列的形式返回，其尺寸與柵格本身相同。此時，你的in_data變數包含一個二維像元值陣列：

in_data = in_band.ReadAsArray()

現在，因為Landsat影象的波段1是藍色波段，你需要將其放入輸出影象的第三個波段以獲得RGB順序的波段。接下來要做的是從out_ds獲取第三個波段，然後使用WriteArray將in_data陣列中的值複製到新資料集的第三個波段：

out_band = out_ds.GetRasterBand(3)
out_band.WriteArray(in_data)

你仍然需要將綠色和紅色Landsat波段新增到資料集中，然後開啟第二個波段的GeoTIFF。請注意，你沒有從資料集中獲取band物件，因為這次你將直接從資料集本身讀取像元資料。因為第二個Landsat波段是綠色波段，你可以獲得堆疊資料集中第二個（綠色）波段的控制代碼，並將Landsat檔案中的資料複製到堆疊資料集：

in_ds = gdal.Open(band2_fn)
out_band = out_ds.GetRasterBand(2)
out_band.WriteArray(in_ds.ReadAsArray())

當你在資料集上呼叫ReadAsArray時，如果你正在讀取的資料集具有多個波段，則會獲得三維陣列。由於Landsat檔案只有一個波段，因此資料集上的ReadAsArray返回與波段物件相同的二維陣列。而不是將資料儲存到中間變數，這次你立即將其傳送到輸出波段。然後，你對紅色波段像元值執行相同的操作，但將其壓縮為更少的程式碼。但效果是一樣的：

out_ds.GetRasterBand(1).WriteArray(gdal.Open(band3_fn).ReadAsArray())

在下一段程式碼中，你可以計算資料集中每個波段的統計資訊。這不是必須的，但它使一些軟體更容易顯示它。統計資料包括平均值，最小值，最大值和標準差。GIS可以使用此資訊來拉伸螢幕上的資料並使其看起來更好。你將在後面的章節中看到如何手動拉伸資料的示例。在計算統計資料之前，你必須確保資料已寫入磁碟而不是僅快取在記憶體中，因此這是對FlushCache的呼叫。然後迴圈遍歷波段並計算每個波段的統計資料。將False傳遞給此函式會告訴你需要實際統計資訊而不是估計值，它可能來自概覽圖層（尚不存在）或從像元子集中取樣。如果估計值可以接受，那麼你可以傳遞True；這也將使計算更快，因為不是每個像元都需要檢查：

out_ds.FlushCache()
for i in range(1, 4):
  out_ds.GetRasterBand(i).ComputeStatistics(False)

你要做的最後一件事是為資料集構建概覽圖層。由於這些像元值是連續資料，因此使用平均插值而不是預設的最近鄰法。你還可以指定要構建的五個級別的概檢視。碰巧有五個級別是你需要為此特定影象獲取大小為256的影象塊：

out_ds.BuildOverviews('average', [2, 4, 8, 16, 32])

最後，不要忘記刪除輸出的資料集。當變數超出範圍時，這將自動發生，但如果你使用互動式Python環境，則可能不會在指令碼完成執行時發生。在平時處理資料時，這經常會發生。它們不會重新整理快取或刪除變數，並且當指令碼完成時，它們的IDE不會釋放資料集物件，因此它們最終會顯示空影象並且不知道是什麼原因。

最後來張LandSat高清彩色合成圖：