基於Landsat-8 OLI影像的鄱陽湖資訊提取（python實現）
一、背景
之前寫過基於雙峰閾值分割的冰湖提取演算法，近期需要做一個湖泊提取的簡單程式，就以鄱陽湖為例吧。本文從零開始介紹如何提取鄱陽湖資訊，並製作shp檔案。
二、資料獲取及預處理
為了獲取鄱陽湖的Landsat-8 OLI影像，首先需要知道鄱陽湖的位置，利用百度直接搜尋，可以查詢到鄱陽湖的經緯度資訊：
鄱陽湖位於北緯28°22′至29°45′，東經115°47′至116°45′。根據鄱陽湖的經緯度資訊，查詢其相應的Landsat-8 OLI的行列號資訊，查詢地址為：https://landsat.usgs.gov/wrs-2-pathrow-latitudelongitude-converter
查詢結果如下：
我用北緯29°和東經116°進行查詢，可知，鄱陽湖所在地區的行列號約為121/40左右，為了下載資料方便，這裡我使用了地理空間資料雲進行資料下載，當然也可以用USGS，只不過慢一點，開啟地理空間資料雲進行資料的查詢：
從縮圖種可以看出，行列號為121/40的影像幾乎包含了鄱陽湖的所有範圍，因此我們只下載這一景影像即可，下載時候需注意，雲量要儘可能的少，因此我找了2017年11月01日的Landsat-8 OLI影像，運量為0.04，直接下載即可。
下載好並解壓之後，我們可以在ENVI軟體中開啟看一下：
影象的質量還不錯，水體非常清晰，因此接下來的工作是關於影象預處理了。
（1）像元DN值轉TOA（top of atmosphere）
因為像元DN值沒有任何物理意義，只表示了影象的明暗程度。我們常常使用的水體指數，其計算的是波段反射率，因此我們需要將DN值轉換為TOA，一般來說，地表引數真實反演的情況下需要計算地表真實反射率，只是提取水體的話，只需要計算感測器入瞳處的反射率，即天頂反射率TOA即可。計算的公式為：
其中，gain表示增益，offset表示偏置，這兩個引數都能夠從Landsat-8的標頭檔案中查詢到：
可知，所有可見光波段的gain=0.00002，offset=-0.1，在ENVI中對其進行轉換即可。轉換完成後檢視其各個波段的反射率統計，反射率基本介於[0 1]說明轉換成功。
（2）影象裁剪
為了去除影象中的冗餘資訊，需要對影象進行簡單裁剪。這裡我只是簡單的用了一下矩形裁剪。
儲存方式選擇tiff即可：
儲存好之後，即可開始我們的提取試驗了。
三、載入影象並呼叫方法進行提取main.py
首先定義一個main.py檔案作為主要檔案，這個檔案中需要定義兩個函式，一個是載入影像函式，一個是根據提取湖泊二值圖製作shp檔案的函式。
首先需要讀取預處理後的tiff檔案，相應的程式碼為：
from osgeo import ogr, osr # 匯入處理shp檔案的庫
from osgeo import gdal, gdal_array # 匯入讀取遙感影像的庫
from NDWI import * # 匯入湖泊提取方法，這裡以NDWI為例
def main(img_dir):
# 載入影像，使用gdal將其載入到numpy中
img = gdal_array.LoadFile(img_dir)
# 呼叫湖泊提取方法，返回一個二值影像
extracted_img = NDWI(img)
# 儲存二值影像
gdal_array.SaveArray(extracted_img.astype(gdal_array.numpy.uint8),
'extracted_img.tif', format="GTIFF", prototype='')
# 對提取的結果進行去噪處理，並將其轉化為shp檔案
raster2shp()
這裡需要注意的是NDWI是接下來定義的湖泊提取方法，因為先寫的main.py檔案，所以即使這裡先報錯也不要急，後面會補上NDWI檔案。
四、編寫湖泊提取演算法NDWI.py
一般的NDWI是利用的green和NIR波段，也就是Landsat-8 的第三和第五波段，根據這個原理，我們直接載入第三和第五波段來計算NDWI，這裡需要注意的就是計算完的NDWI中有很多噪音，因此這裡考慮利用形態學開運算去除噪音。下面直接給出程式碼：
import numpy as np
import cv2
def NDWI(img, threshold=0.4):
"""
該函式採用最簡單的NDWI進行水體的提取。
需要輸入的引數為載入好的影像img和閾值threshold
返回為提取好的水體掩模
"""
# NDWI用到了Landsat 8 OLI的第3和第5波段，先找到這兩個波段
green = img[2]
nir = img[4]
# 計算NDWI並建立掩模
ndwi = (green - nir) / (green + nir)
# 根據閾值來確定掩模的值
for row in range(ndwi.shape[0]):
for col in range(ndwi.shape[1]):
if ndwi[row, col] >= threshold:
# ndwi_mask[row, col] = 1
ndwi[row, col] = 1
else:
ndwi[row, col] = 0
# 最後對影象進行開運算進行去噪，即先腐蝕後膨脹
kernel = np.ones((5, 5), np.uint8)
opening = cv2.morphologyEx(ndwi, cv2.MORPH_OPEN, kernel)
return opening
五、提取湖泊並製作shp檔案main.py
最後一步是我們根據NDWI提取好的掩模結果來製作shp檔案，之前我是matlab做這個部分的，如果是matlab可以考慮用contour製作，python的製作方法我是參考我之前的文章：遙感影象處理中常用的python操作中的第五節來做的，基本沒有做任何改動，下面直接給出main.py後面部分的程式碼：
def raster2shp(src="extracted_img.tif"):
"""
函式輸入的是一個二值影像，利用這個二值影像，建立shp檔案
"""
# src = "extracted_img.tif"
# 輸出的shapefile檔名稱
tgt = "extract.shp"
# 圖層名稱
tgtLayer = "extract"
# 開啟輸入的柵格檔案
srcDS = gdal.Open(src)
# 獲取第一個波段
band = srcDS.GetRasterBand(1)
# 讓gdal庫使用該波段作為遮罩層
mask = band
# 建立輸出的shapefile檔案
driver = ogr.GetDriverByName("ESRI Shapefile")
shp = driver.CreateDataSource(tgt)
# 拷貝空間索引
srs = osr.SpatialReference()
srs.ImportFromWkt(srcDS.GetProjectionRef())
layer = shp.CreateLayer(tgtLayer, srs=srs)
# 建立dbf檔案
fd = ogr.FieldDefn("DN", ogr.OFTInteger)
layer.CreateField(fd)
dst_field = 0
# 從圖片中自動提取特徵
extract = gdal.Polygonize(band, mask, layer, dst_field, [], None)
if __name__ == '__main__':
main('poyangLake.tif')六、試驗結果
製作好main.py和NDWI.py，就可以直接執行main.py檔案，執行的結果會生成一個二值掩模影象和shp檔案，在ENVI中開啟最終的效果：
效果算不上好，還有很大的提升空間。後續有空我再進行優化。
七、分析
1.雖然寫的方法比較簡單，不過提供了一種思路。這裡我認為可以改進的地方有：（1）是湖泊提取方法的改進，NDWI的閾值分割法是最最初級的湖泊提取方法，對於眾多的水體指數而言思路的差別不大。（2）提取結果的後處理，可以看到提取結果非常破碎，很多細小的河流被提取了出來，這部分可以考慮用形狀指數進行限制，還有就是湖泊內部的很多斑點，需要進一步填充處理。
2.製作shp檔案的地理參照過程，我個人感覺還有一點小問題，這似乎與你開啟影像的左上角座標有關。個人感覺可以利用contour來繪製湖泊邊界。另外在預處理過程中，對原始影像進行裁剪的時候，裁剪資料的右上角其實是空值，在後續進行處理的時候需要對其置空處理。關於使用contour方法提取二值影象邊界的程式碼：
# 可以參考https://www.cnblogs.com/denny402/p/5160955.html
from skimage import measure
# 檢測所有圖形的輪廓
contours = measure.find_contours(img, 0.5)3.如果改進了演算法，也可以直接套用這個思路，需要改的地方只有main.py檔案中的第3和第11行
4.所有檔案的結構為：
-- poyangLake.tif # 預處理好的影像資料
-- main.py
{
from xxx import...
def main(img_dir):...
def raster2shp(src="extracted_img.tif"):...
if __name__ == '__main__':...
}
-- NDWI.py
{
import...
def NDWI(img, threshold=0.8):...
}