一、背景
之前写过基于双峰阈值分割的冰湖提取算法,近期需要做一个湖泊提取的简单程序,就以鄱阳湖为例吧。本文从零开始介绍如何提取鄱阳湖信息,并制作shp文件。
二、数据获取及预处理
为了获取鄱阳湖的Landsat-8 OLI影像,首先需要知道鄱阳湖的位置,利用百度直接搜索,可以查询到鄱阳湖的经纬度信息:
鄱阳湖位于北纬28°22′至29°45′,东经115°47′至116°45′。根据鄱阳湖的经纬度信息,查询其相应的Landsat-8 OLI的行列号信息,查询地址为:https://landsat.usgs.gov/wrs-2-pathrow-latitudelongitude-converter
查询结果如下:
2019-12-30更新:
上述查询地址似乎失效了,可以用这个查询地址:https://landsat.usgs.gov/landsat_acq#convertPathRow
输入经纬度点击convert就能转换了:
我用北纬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):...
}