gdal+python计算NDVI

---

前言

本文采用gdal+python 对遥感影像进行处理，并计算ndvi

---

1.引入库

代码如下：

import pathlib
import numpy as np
from osgeo import gdal

2.读入数据

利用GDAL读取数据：

def read_img(filename):
    dataset = gdal.Open(filename)  # 打开文件
    if dataset == None:
        raise Exception(f"cant find/open {filename}")
    im_width = dataset.RasterXSize  # 栅格矩阵的列数
    im_height = dataset.RasterYSize  # 栅格矩阵的行数
    im_Band = dataset.RasterCount  # 栅格矩阵的行数

    im_geotrans = dataset.GetGeoTransform()
    im_proj = dataset.GetProjection()
    im_data = dataset.ReadAsArray(0, 0, im_width, im_height)  # 将数据写成数组，对应栅格矩阵

    del dataset  # 关闭对象，文件dataset
    return im_proj, im_geotrans, im_data, im_width, im_height, im_Band

3.计算NDVI

利用GDAL读取数据
$$NDVI=\frac{NIR-RED}{NIR+RED}.$$

def read_img(filename):
    dataset = gdal.Open(filename)  # 打开文件
    if dataset == None:
        raise Exception(f"cant find/open {filename}")
    im_width = dataset.RasterXSize  # 栅格矩阵的列数
    im_height = dataset.RasterYSize  # 栅格矩阵的行数
    im_Band = dataset.RasterCount  # 栅格矩阵的行数

    im_geotrans = dataset.GetGeoTransform()
    im_proj = dataset.GetProjection()
    im_data = dataset.ReadAsArray(0, 0, im_width, im_height)  # 将数据写成数组，对应栅格矩阵

    del dataset  # 关闭对象，文件dataset
    return im_proj, im_geotrans, im_data, im_width, im_height, im_Band

---

3.写入数据

利用GDAL写入数据：

DType2GDAL = {"uint8": gdal.GDT_Byte,
              "uint16": gdal.GDT_UInt16,
              "int16": gdal.GDT_Int16,
              "uint32": gdal.GDT_UInt32,
              "int32": gdal.GDT_Int32,
              "float32": gdal.GDT_Float32,
              "float64": gdal.GDT_Float64,
              "cint16": gdal.GDT_CInt16,
              "cint32": gdal.GDT_CInt32,
              "cfloat32": gdal.GDT_CFloat32,
              "cfloat64": gdal.GDT_CFloat64}


def write_img(filename, im_proj, im_geotrans, im_data):
    # 判断栅格数据的数据类型
    if im_data.dtype.name in DType2GDAL:
        datatype = DType2GDAL[im_data.dtype.name]
    else:
        datatype = gdal.GDT_Float32
    # 判读数组维数
    if len(im_data.shape) == 3:
        im_bands, im_height, im_width = im_data.shape
    else:
        im_bands, (im_height, im_width) = 1, im_data.shape
    # 创建文件
    if not pathlib.Path(filename).parent.exists():
        pathlib.Path(filename).parent.mkdir(parents=True, exist_ok=True)
    driver = gdal.GetDriverByName("GTiff")  # 数据类型必须有，因为要计算需要多大内存空间
    dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

    dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
    dataset.SetProjection(im_proj)  # 写入投影
    if im_bands == 1:
        dataset.GetRasterBand(1).WriteArray(im_data)  # 写入数组数据
    else:
        for i in range(im_bands):
            dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
    del dataset

---

4.完整代码

# -*- coding: utf-8 -*-
"""
 @Time    : 2022/10/23 16:49
 @Author  : QYD
 @About   :
"""

import pathlib

import numpy as np
from osgeo import gdal

DType2GDAL = {"uint8": gdal.GDT_Byte,
              "uint16": gdal.GDT_UInt16,
              "int16": gdal.GDT_Int16,
              "uint32": gdal.GDT_UInt32,
              "int32": gdal.GDT_Int32,
              "float32": gdal.GDT_Float32,
              "float64": gdal.GDT_Float64,
              "cint16": gdal.GDT_CInt16,
              "cint32": gdal.GDT_CInt32,
              "cfloat32": gdal.GDT_CFloat32,
              "cfloat64": gdal.GDT_CFloat64}


def write_img(filename, im_proj, im_geotrans, im_data):
    # 判断栅格数据的数据类型
    if im_data.dtype.name in DType2GDAL:
        datatype = DType2GDAL[im_data.dtype.name]
    else:
        datatype = gdal.GDT_Float32
    # 判读数组维数
    if len(im_data.shape) == 3:
        im_bands, im_height, im_width = im_data.shape
    else:
        im_bands, (im_height, im_width) = 1, im_data.shape
    # 创建文件
    if not pathlib.Path(filename).parent.exists():
        pathlib.Path(filename).parent.mkdir(parents=True, exist_ok=True)
    driver = gdal.GetDriverByName("GTiff")  # 数据类型必须有，因为要计算需要多大内存空间
    dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

    dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
    dataset.SetProjection(im_proj)  # 写入投影
    if im_bands == 1:
        dataset.GetRasterBand(1).WriteArray(im_data)  # 写入数组数据
    else:
        for i in range(im_bands):
            dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
    del dataset


def read_img(filename):
    dataset = gdal.Open(filename)  # 打开文件
    if dataset == None:
        raise Exception(f"cant find/open {filename}")
    im_width = dataset.RasterXSize  # 栅格矩阵的列数
    im_height = dataset.RasterYSize  # 栅格矩阵的行数
    im_Band = dataset.RasterCount  # 栅格矩阵的行数

    im_geotrans = dataset.GetGeoTransform()
    im_proj = dataset.GetProjection()
    im_data = dataset.ReadAsArray(0, 0, im_width, im_height)  # 将数据写成数组，对应栅格矩阵

    del dataset  # 关闭对象，文件dataset
    return im_proj, im_geotrans, im_data, im_width, im_height, im_Band


def NDVI(redBand, nirBand):
    redBand = redBand.astype(np.float32)
    nirBand = nirBand.astype(np.float32)
    return (nirBand - redBand) / (nirBand + redBand)


if __name__ == '__main__':
    imgData = r"E:\pycharm\PycharmProjects\classify2\classify\Data\GF_test.tif"
    outPath = "NDVI.tif"
    im_proj, im_geotrans, im_data, im_width, im_height, im_Band = read_img(imgData)
    ndvi = NDVI(im_data[2, :, :], im_data[3, :, :])
    write_img(outPath, im_proj, im_geotrans, ndvi)

---