利用归一化植被指数ndvi实现植被分类制图(gdal实现)

遥感与GIS分析 专栏收录该内容
6 篇文章 9 订阅

一、背景

本文主要参考《python地理空间分析指南》。目的在于利用归一化植被指数NDVI来实现植被的分类制图。分类的方法比较简单,就是阈值切片。我对书中的内容进行了改进,并在分析中会指出书中存在的问题,欢迎进一步探讨。

二、数据

进行实验的话数据很重要,本实验的数据是tiff格式的影像和一个shp矢量文件。原始数据的下载地址为:

http://git.io/v3fS9

下载好的数据放到项目文件夹下,直接解压即可。下面我们来看一下下载到的数据是什么样的:

下载到的数据虽然看着很多,但其实我们能够用到的只有两个,即原始影像数据文件farm.tif和农田的矢量文件field.shp,用ENVI软件打开这两个文件:

可以看到:原始数据中的farm.tif是3波段数据(打开波段查看器可以看的更清楚),而且田中的田地和水体都是真彩色影像,所以可以推测出farm.tif的三个波段对应的是RGB三个波段;而field.shp对应的则是图中间的一小块区域;另外查看像元值可以知道,图像为8bit图像,各个波段的灰度范围都是0-255。同时该图像还带了空间投影及地理信息,是一幅比较典型的遥感影像。

三、实验

这里需要重点强调的一点是,实验采用的植被指数NDVI的计算公式为:

                                                                         NDVI=\frac{\rho _{red}-\rho _{NIR}}{\rho _{red}+\rho _{NIR}}

即利用红光和近红外波段的反射率进行运算。

但是,我个人认为原实验中有两处小错误:①NDVI的计算公式利用的是反射率而非图像DN值,用DN值算下来的结果其实没有任何物理意义,实际上我们在进行计算的时候,根据需求,往往采用入瞳处的反射率或者地表真实反射率,计算这两种反射率都需要用到头文件;②原始图像应该是真彩色图像,所以是根本没有近红外波段的。

不过这只是一个简单的小实验,不影响我们的处理思路。

总体思路只有两步,计算NDVI和影像分类:

1.计算NDVI

这一步主要的目的是利用gdal_array将农田影像载入numpy数组(gdal_array只能载入图像,而gdal可以获取图像的地理参照),然后分别获取近红外和红光波段;之后是将shp文件栅格化,填充并转化为numpy数组;最后是裁剪相应的波段,计算NDVI并保存。

先直接给出主要代码(ndvi.py):

from osgeo import gdal_array
from PIL import Image, ImageDraw


def imageToArray(i):
    """
    function: 将PIL图像库的数组转换为一个gdal_array图片
    """
    a = gdal_array.numpy.fromstring(i.tobytes(), 'b')
    a.shape = i.im.size[1], i.im.size[0]
    return a


def world2Pixel(geoMatrix, x, y):
    """
    function: 该函数使用gdal的geoMatrix方法,将目标图片的地理坐标转换成像素坐标
    """
    ulX = geoMatrix[0]
    ulY = geoMatrix[3]
    xDist = geoMatrix[1]
    yDist = geoMatrix[5]
    pixel = int((x - ulX) / xDist)
    line = int((ulY - y) / abs(yDist))
    return (pixel, line)


import gdalnumeric

# 获取数据部分-----------------------------------------------------
# 需要读取的影像
source = "./NDVI-update/NDVI-update/farm.tif"

# 最终生成的影像
target = "ndvi.tif"

# 将载入的数据转换成gdal_array数组
srcArray = gdal_array.LoadFile(source)

# 同时将源数据作为gdal图片载入,获取地理信息参照
srcImage = gdal.Open(source)
geoTrans = srcImage.GetGeoTransform()

# 获取红光和近红外波段
r = srcArray[1]
ir = srcArray[2]        # 这个应该是蓝光波段

# 裁剪影像部分-----------------------------------------------------
# 根据shapefile文件创建一个OGR图层
field = ogr.Open("./NDVI-update/NDVI-update/field.shp")
# 定义一个图层确保OGR格式完好
lyr = field.GetLayer("field")
# shapefile文件中只有一个多边形
poly = lyr.GetNextFeature()

# 将图层坐标转化为图片的像素坐标
minX, maxX, minY, maxY = lyr.GetExtent()
ulX, ulY = world2Pixel(geoTrans, minX, maxY)    # 左上角坐标
lrX, lrY = world2Pixel(geoTrans, maxX, minY)    # 右下角坐标

# 计算新图片的像素尺寸
pxWidth = int(lrX - ulX)
pxHeight = int(lrY - ulY)

# 为遮罩层创建一个合适尺寸的空白图片
clipped = gdal_array.numpy.zeros((3, pxHeight, pxWidth), gdal_array.numpy.uint8)

# 裁剪波段
rclip = r[ulY:lrY, ulX:lrX]
irclip = ir[ulY:lrY, ulX:lrX]

# 为图片创建一个新的geomatrix对象
geoTrans = list(geoTrans)
geoTrans[0] = minX
geoTrans[3] = maxY

# 在8位黑白遮罩图片上把点映射为像素用于绘制区域边界
points = []
pixels = []
# 获取多边形几何图形
geom = poly.GetGeometryRef()
pts = geom.GetGeometryRef(0)

# 遍历几何图像,将点转换为易于管理的list对象
for p in range(pts.GetPointCount()):
    points.append((pts.GetX(p), pts.GetY(p)))

# 遍历点集并将之映射为像素,然后把这些像素添加到list中
for p in points:
    pixels.append(world2Pixel(geoTrans, p[0], p[1]))

# 用“L”模式创建一个栅格多边形的黑白图片,白色值=1
rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)

# 创建一个PIL绘制对象
rasterize = ImageDraw.Draw(rasterPoly)

# 将像素点存储为图片上的多边形, 黑色=0
rasterize.polygon(pixels, 0)

# 将图片转换为gdal_array以便将其当作一个遮罩数组来用
mask = imageToArray(rasterPoly)

# 裁剪波段---------------------------------------------------------------
# 使用遮罩分别裁剪红光和近红外波段
rclip = gdal_array.numpy.choose(mask, (rclip, 0)).astype(gdal_array.numpy.uint8)
irclip = gdal_array.numpy.choose(mask, (irclip, 0)).astype(gdal_array.numpy.uint8)

# 计算NDVI---------------------------------------------------------------
# 忽略裁剪过程中numpy遇到的任何NaN值
gdal_array.numpy.seterr(all='ignore')

# 按照公式计算
ndvi = 1.0 * ((irclip - rclip) / (irclip + rclip + 1.0))

# 将结果中的所有NaN值转换为0
ndvi = gdal_array.numpy.nan_to_num(ndvi)

# 将NDVI图像进行保存,并调整地理参照信息
gdalnumeric.SaveArray(ndvi, target, format="GTiff", prototype=srcImage)

# 并调整地理参照信息
update = gdal.Open(target, 1)
update.SetGeoTransform(list(geoTrans))
update.GetRasterBand(1).SetNoDataValue(0.0)
update = None

直接执行上述程序,会生成计算好的NDVI影像,下面直接展示结果:

对于结果进行了分析,NDVI大部分的值还是介于0-1.2之间的,但是仍有像元的值为150。。。。。。考虑了一下,可能是程序中的计算公式,为了防止分母为0,在分母中又多加1,这样就有可能出现异常值。

2.影像分类

这一步并没有用分类方法对其进行分类,而是使用棕色到深绿色7个类别对其进行分类。首先读取处理好的ndvi影像,然后进行直方图均衡化处理;之后是创建7个类别,对影像数据进行分类。

下面直接给出分类的代码(classification.py):

from osgeo import gdal_array as gd
import operator
from functools import reduce

def histogram(a, bins=list(range(256))):
    """
    多维数组的直方图函数
    """
    # 直方图的数组均衡化、排序和分割
    fa = a.flat
    n = gd.numpy.searchsorted(gd.numpy.sort(fa), bins)
    n = gd.numpy.concatenate([n, [len(fa)]])
    hist = n[1:] - n[:-1]
    return hist

def stretch(a):
    """
    在gdal_array数组图像上执行直方图平衡化操作
    """
    hist = histogram(a)
    lut = []
    for b in range(0, len(hist), 256):
        # 创建等间隔区间
        step = reduce(operator.add, hist[b:b+256]) / 255
        # 创建均衡化查找表
        n = 0
        for i in range(256):
            lut.append(n / step)
            n += hist[i+b]
    gd.numpy.take(lut, a, out=a)
    return a

# 载入数据-----------------------------------------------------------
# 原图片为NDVI脚本输出
source = "./ndvi.tif"

# 分类图片的目标文件名
target = "ndvi_color.tif"

# 将图片载入数组
ndvi = gd.LoadFile(source).astype(gd.numpy.uint8)

# 预解析NDVI
# 执行直方图均衡化操作以便能够使用所有类别
ndvi = stretch(ndvi)

# 创建一个尺寸大小和NDVI相同的3波段图像
rgb = gd.numpy.zeros((3, len(ndvi), len(ndvi[0])), gd.numpy.uint8)

# 创建类-------------------------------------------------------------
# 类别列表中NDVI的区间上限值,需要注意上限和下限值在列表两端
classes = [58, 73, 110, 147, 184, 220, 255]

# 颜色查找表。lut必须从深棕色到深绿色代表的RGB元组匹配
lut = [[120, 69, 25], [255, 178, 74], [255, 237, 166], [173, 232, 94],
       [135, 181, 64], [3, 156, 0], [1, 100, 0]]

# 第一种类别的起始值
start = 1

# 获取每个类别的区间值,获取区间值,然后使用遮罩过滤
for i in range(len(classes)):
    mask = gd.numpy.logical_and(start <= ndvi, ndvi <= classes[i])
    for j in range(len(lut[i])):
        rgb[j] = gd.numpy.choose(mask, (rgb[j], lut[i][j]))
    start = classes[i] + 1

# 保存彩色的NDVI数据为geotiff图片
output = gd.SaveArray(rgb, target, format="GTiff", prototype=source)
output = None

直接执行上述代码,最终的实验结果为:

个人感觉实验的结果并不好,一眼看过去,主要都是翠绿色,和一点墨绿色,其他颜色没有。。。虽然分了7类,不过这个分类的实验瑕疵挺大的。

四、分析

1.这个实验的主要问题表现在两个方面:①NDVI的计算公式利用的是反射率而非图像DN值,而原始实验中采用的是像元DN值②原始图像应该是真彩色图像,所以是根本没有近红外波段的。

2.关于分类,一般比较简单的是采用阈值分割,不过这个阈值切片的分类效果是真的很差。

3.实验的目的主要在于提供图像处理的思路,不过效果是真的不好。

  • 7
    点赞
  • 6
    评论
  • 66
    收藏
  • 一键三连
    一键三连
  • 扫一扫,分享海报

©️2021 CSDN 皮肤主题: 技术黑板 设计师:CSDN官方博客 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值