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

最新推荐文章于 2024-09-20 16:06:30 发布

全部梭哈迟早暴富

最新推荐文章于 2024-09-20 16:06:30 发布

阅读量2.5w

点赞数 11

分类专栏：遥感与GIS分析

本文链接：https://blog.csdn.net/z704630835/article/details/83305746

版权

遥感与GIS分析专栏收录该内容

6 篇文章 42 订阅

订阅专栏

一、背景

本文主要参考《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

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