基于Landsat-8 OLI影像的鄱阳湖信息提取(python实现)

基于Landsat-8 OLI影像的鄱阳湖信息提取(python实现)一、背景之前写过基于双峰阈值分割的冰湖提取算法,近期需要做一个湖泊提取的简单程序,就以鄱阳湖为例吧。本文从零开始介绍如何提取鄱阳湖信息,并制作shp文件。二、数据获取及预处理为了获取鄱阳湖的Landsat-8OLI影像,首先需要知道鄱阳湖的位置,利用百度直接搜索,可以查询到鄱阳湖的经纬度信息:鄱阳湖位于北纬28°22′至29°45′,东经115°47′至116°45′。根据鄱…

一、背景

之前写过基于双峰阈值分割的冰湖提取算法,近期需要做一个湖泊提取的简单程序,就以鄱阳湖为例吧。本文从零开始介绍如何提取鄱阳湖信息,并制作shp文件。

二、数据获取及预处理

为了获取鄱阳湖的Landsat-8 OLI影像,首先需要知道鄱阳湖的位置,利用百度直接搜索,可以查询到鄱阳湖的经纬度信息:

基于Landsat-8 OLI影像的鄱阳湖信息提取(python实现)

鄱阳湖位于北纬28°22′至29°45′,东经115°47′至116°45′。根据鄱阳湖的经纬度信息,查询其相应的Landsat-8 OLI的行列号信息,查询地址为:https://landsat.usgs.gov/wrs-2-pathrow-latitudelongitude-converter

查询结果如下:

基于Landsat-8 OLI影像的鄱阳湖信息提取(python实现)


2019-12-30更新:

上述查询地址似乎失效了,可以用这个查询地址:https://landsat.usgs.gov/landsat_acq#convertPathRow

输入经纬度点击convert就能转换了:

基于Landsat-8 OLI影像的鄱阳湖信息提取(python实现)


我用北纬29°和东经116°进行查询,可知,鄱阳湖所在地区的行列号约为121/40左右,为了下载数据方便,这里我使用了地理空间数据云进行数据下载,当然也可以用USGS,只不过慢一点,打开地理空间数据云进行数据的查询:

基于Landsat-8 OLI影像的鄱阳湖信息提取(python实现)

从缩略图种可以看出,行列号为121/40的影像几乎包含了鄱阳湖的所有范围,因此我们只下载这一景影像即可,下载时候需注意,云量要尽可能的少,因此我找了2017年11月01日的Landsat-8 OLI影像,运量为0.04,直接下载即可。

下载好并解压之后,我们可以在ENVI软件中打开看一下:

基于Landsat-8 OLI影像的鄱阳湖信息提取(python实现)

图像的质量还不错,水体非常清晰,因此接下来的工作是关于图像预处理了。

(1)像元DN值转TOA(top of atmosphere)

因为像元DN值没有任何物理意义,只表示了图像的明暗程度。我们常常使用的水体指数,其计算的是波段反射率,因此我们需要将DN值转换为TOA,一般来说,地表参数真实反演的情况下需要计算地表真实反射率,只是提取水体的话,只需要计算传感器入瞳处的反射率,即天顶反射率TOA即可。计算的公式为:

                                                                      TOA = gain\times DN + offset

其中,gain表示增益,offset表示偏置,这两个参数都能够从Landsat-8的头文件中查找到:

基于Landsat-8 OLI影像的鄱阳湖信息提取(python实现)

可知,所有可见光波段的gain=0.00002,offset=-0.1,在ENVI中对其进行转换即可。转换完成后查看其各个波段的反射率统计,反射率基本介于[0 1]说明转换成功。

(2)图像裁剪

为了去除图像中的冗余信息,需要对图像进行简单裁剪。这里我只是简单的用了一下矩形裁剪。

基于Landsat-8 OLI影像的鄱阳湖信息提取(python实现)

保存方式选择tiff即可:

基于Landsat-8 OLI影像的鄱阳湖信息提取(python实现)

保存好之后,即可开始我们的提取试验了。

三、加载图像并调用方法进行提取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中打开最终的效果:

基于Landsat-8 OLI影像的鄱阳湖信息提取(python实现)

效果算不上好,还有很大的提升空间。后续有空我再进行优化。

七、分析

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):...
    }

 

今天的文章基于Landsat-8 OLI影像的鄱阳湖信息提取(python实现)分享到此就结束了,感谢您的阅读,如果确实帮到您,您可以动动手指转发给其他人。

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。
如需转载请保留出处:https://bianchenghao.cn/25505.html

(0)
编程小号编程小号

相关推荐

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注