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

該公眾號將為您奉上Python地學分析、爬蟲、數據分析、Web開發、機器學習、深度學習等熱門源代碼。

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

https://github.com/xbr2017

CSDN也在同步更新：

https://blog.csdn.net/XBR_2014

---

“?遙感影像的特點之一就是同一個圖像文件可以存儲若干個波段的圖像，本節內容主要介紹GDAL將多個尺寸相同的圖像疊加到一起，以GeoTiff格式輸出，這有利于不同波段之間進行數值運算?！?/p>

本節以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數據類型常量

常量	數據類型
GDT_Unknown	未知
GDT_Byte	無符號8位整數（字節）
GDT_UInt16	無符號16位整數
GDT_Int16	有符號的16位整數
GDT_UInt32	無符號32位整數
GDT_Int32	有符號的32位整數
GDT_Float32	32位浮點
GDT_Float64	64位浮點
GDT_CInt16	16位復數整數
GDT_CInt32	32位復數整數
GDT_CFloat32	32位復雜浮點
GDT_CFloat64	64位復雜浮點
GDT_TypeCount	可用數據類型的數量

你從輸入數據集中獲取投影（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高清彩色合成圖：

?

?

總結

以上是生活随笔為你收集整理的Python地学分析 — GDAL将多个遥感图像叠加保存为tif文件的全部內容，希望文章能夠幫你解決所遇到的問題。

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