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

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

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

https://github.com/xbr2017

CSDN也在同步更新：

https://blog.csdn.net/XBR_2014

---

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

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

GDAL實現(xiàn)多波段疊加

現(xiàn)在理論已經不在了，讓我們學習如何使用GDAL處理這些數(shù)據(jù)集。柵格數(shù)據(jù)存在許多不同的文件格式，GDAL是一個非常流行且強大的庫，用于讀寫許多文件格式。GDAL庫是開源的，但具有許可，因此即使許多商業(yè)軟件包也使用它。

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

為了說明如何使用GDAL讀取和寫入柵格數(shù)據(jù)，讓我們從一個將三個Landsat波段組合成一個疊加圖像的示例開始，如圖1所示。

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

Landsat計劃是美國地質調查局（USGS）與美國國家航空航天局（NASA）之間的聯(lián)合倡議，自1972年以來一直在全球范圍內采集中等分辨率的衛(wèi)星圖像。Landsat圖像由USGS作為集合分發(fā)GeoTIFFs，每個數(shù)據(jù)集均有波段。除了波段6（熱）和8（全色）之外，每個波段都有30米的分辨率，因為它們來自同一個Landsat場景，尺寸相同。這使得事情變得簡單，因為波段直接相互疊加而不需要重采樣之類的處理。下面的程序實現(xiàn)了如何創(chuàng)建具有相同尺寸的三波段數(shù)據(jù)集，然后將波段3，2和1復制到其中。這三個波段分別對應于可見光的紅色，綠色和藍色波長，因此將它們按此順序排列將產生RGB（紅色，綠色，藍色）圖像，其外觀與你自己的眼睛非常相似。圖1顯示了黑色和白色的單獨紅色，藍色和綠色波段，以及生成的三波段自然彩色圖像。如果你在ArcGIS中預覽圖像，它們可能看起來與此類似，因為ArcGIS很可能會拉伸它們。下面來看看Python代碼的實現(xiàn)：

#?_*_?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波段數(shù)據(jù) in_data?=?in_band.ReadAsArray() out_band?=?out_ds.GetRasterBand(3) out_band.WriteArray(in_data)#?讀取第2波段數(shù)據(jù) in_ds?=?gdal.Open(band2_fn) out_band?=?out_ds.GetRasterBand(2) out_band.WriteArray(in_ds.ReadAsArray())#?讀取第3波段數(shù)據(jù) 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。你還可以獲取數(shù)據(jù)集中第一個也是唯一一個波段的句柄，盡管你尚未讀取任何數(shù)據(jù)。請注意，使用索引1而不是0來獲取第一個波段。當你使用GetRasterBand時，帶編號始終以1開頭，很多人會忘記這個知識點，容易出錯。無論如何，在創(chuàng)建輸出圖像之前，你需要此波段對象，因為它具有輸出文件所需要的基本信息（尺寸大小、投影、坐標等）。

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

接下來，你將創(chuàng)建一個新數(shù)據(jù)集以將像元數(shù)據(jù)復制到其中。你必須使用驅動程序對象來創(chuàng)建新數(shù)據(jù)集，因此你可以找到GeoTIFF驅動程序，然后使用其Create函數(shù)。這是該功能的完整參數(shù)：

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

• filename是要創(chuàng)建的數(shù)據(jù)集的路徑。
• xsize是新數(shù)據(jù)集中的列數(shù)。
• ysize是新數(shù)據(jù)集中的行數(shù)。
• bands是新數(shù)據(jù)集中的波段數(shù)。默認值為1。
• data_type是將存儲在新數(shù)據(jù)集中的數(shù)據(jù)類型。默認值為GDT_Byte。
• options是創(chuàng)建選項字符串的列表。可能的值取決于正在創(chuàng)建的數(shù)據(jù)集的類型。

因為你使用GeoTIFF驅動程序，所以無論你提供什么文件擴展名，輸出文件都將是GeoTIFF。但是，擴展名不會自動添加，因此你需要提供該擴展程序。在這種情況下，你將其命名為nat_color.tif并將其保存在指定文件夾中，因為這是使用os.chdir設置的當前文件夾。你還需要提供列數(shù)和行數(shù)創(chuàng)建新數(shù)據(jù)集，因此分別使用XSize和YSize屬性從輸入波段獲取該信息。Open的下一個參數(shù)是band的數(shù)量，你希望這個新的raster有三個。下一個可選參數(shù)是數(shù)據(jù)類型，它必須是表1中的值之一。你可以從輸入波段獲取此信息，但在這種情況下你可能會忽略它，因為這些圖像使用默認類型的DT_Byte。

? ? ? ? ? ? 表1 ?GDAL數(shù)據(jù)類型常量

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

你從輸入數(shù)據(jù)集中獲取投影（SRS）并將其復制到新數(shù)據(jù)集，然后對地理轉換執(zhí)行相同操作。地理轉換很重要，因為它提供原點坐標和像元大小。如上所述，在將數(shù)據(jù)集放置在正確的空間位置時，原點和像元大小非常重要。雖然你不必在添加像元值之前添加投影和地理轉換信息，但在創(chuàng)建新數(shù)據(jù)集后立即將其解除。

設置數(shù)據(jù)集后，就可以添加像元值了。因為你已經擁有Landsat band 1的GeoTIFF中的波段對象，所以你可以將其中的像元值讀入NumPy數(shù)組。如果不向ReadAsArray提供任何參數(shù)，則所有像元值都以二維數(shù)組的形式返回，其尺寸與柵格本身相同。此時，你的in_data變量包含一個二維像元值數(shù)組：

in_data?=?in_band.ReadAsArray()

現(xiàn)在，因為Landsat圖像的波段1是藍色波段，你需要將其放入輸出圖像的第三個波段以獲得RGB順序的波段。接下來要做的是從out_ds獲取第三個波段，然后使用WriteArray將in_data數(shù)組中的值復制到新數(shù)據(jù)集的第三個波段：

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

你仍然需要將綠色和紅色Landsat波段添加到數(shù)據(jù)集中，然后打開第二個波段的GeoTIFF。請注意，你沒有從數(shù)據(jù)集中獲取band對象，因為這次你將直接從數(shù)據(jù)集本身讀取像元數(shù)據(jù)。因為第二個Landsat波段是綠色波段，你可以獲得堆疊數(shù)據(jù)集中第二個（綠色）波段的句柄，并將Landsat文件中的數(shù)據(jù)復制到堆疊數(shù)據(jù)集：

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

當你在數(shù)據(jù)集上調用ReadAsArray時，如果你正在讀取的數(shù)據(jù)集具有多個波段，則會獲得三維數(shù)組。由于Landsat文件只有一個波段，因此數(shù)據(jù)集上的ReadAsArray返回與波段對象相同的二維數(shù)組。而不是將數(shù)據(jù)保存到中間變量，這次你立即將其發(fā)送到輸出波段。然后，你對紅色波段像元值執(zhí)行相同的操作，但將其壓縮為更少的代碼。但效果是一樣的：

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

在下一段代碼中，你可以計算數(shù)據(jù)集中每個波段的統(tǒng)計信息。這不是必須的，但它使一些軟件更容易顯示它。統(tǒng)計數(shù)據(jù)包括平均值，最小值，最大值和標準差。GIS可以使用此信息來拉伸屏幕上的數(shù)據(jù)并使其看起來更好。你將在后面的章節(jié)中看到如何手動拉伸數(shù)據(jù)的示例。在計算統(tǒng)計數(shù)據(jù)之前，你必須確保數(shù)據(jù)已寫入磁盤而不是僅緩存在內存中，因此這是對FlushCache的調用。然后循環(huán)遍歷波段并計算每個波段的統(tǒng)計數(shù)據(jù)。將False傳遞給此函數(shù)會告訴你需要實際統(tǒng)計信息而不是估計值，它可能來自概覽圖層（尚不存在）或從像元子集中采樣。如果估計值可以接受，那么你可以傳遞True；這也將使計算更快，因為不是每個像元都需要檢查：

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

你要做的最后一件事是為數(shù)據(jù)集構建概覽圖層。由于這些像元值是連續(xù)數(shù)據(jù)，因此使用平均插值而不是默認的最近鄰法。你還可以指定要構建的五個級別的概視圖。碰巧有五個級別是你需要為此特定圖像獲取大小為256的圖像塊：

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

最后，不要忘記刪除輸出的數(shù)據(jù)集。當變量超出范圍時，這將自動發(fā)生，但如果你使用交互式Python環(huán)境，則可能不會在腳本完成運行時發(fā)生。在平時處理數(shù)據(jù)時，這經常會發(fā)生。它們不會刷新緩存或刪除變量，并且當腳本完成時，它們的IDE不會釋放數(shù)據(jù)集對象，因此它們最終會顯示空圖像并且不知道是什么原因。

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

?

?

總結

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

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