Python读取遥感影像并计算NDVI指数

Python读取遥感影像并计算NDVI指数

该篇计算NDVI存在问题，查看GDAL计算NDVI及常用函数介绍
先附上源码
读取的遥感影像数据需要先进行预处理

from osgeo import gdal
import numpy as np
import pandas as pd
import os
np.seterr(divide='ignore', invalid='ignore')

class Grid:
    # =============================================================================
    #     1. 函数1 read_tiff读取影像
    # =============================================================================
    def read_tiff(self, filename):
        # 读取影像，获取影像位置
        dataset = gdal.Open(filename)
        # 获取影像波段数,获取影像长宽
        im_band=dataset.GetRasterBand   #波段数
        im_width = dataset.RasterXSize  # 宽，列数
        im_height = dataset.RasterYSize  # 高,行数
        # 仿射矩阵
        im_GeoTransform = dataset.GetGeoTransform()
        # 地图投影信息
        img_proj = dataset.GetProjection()

        im_data = dataset.ReadAsArray(0, 0, im_width, im_height)
        del dataset
        return img_proj, im_GeoTransform, im_data

    # =============================================================================
    #     函数2 write_tiff，存储为GTIFF
    # =============================================================================
    def write_tiff(self, filename, im_proj, im_GeoTransform, im_data):
        # 判断栅格数据类型
        if 'int8' in im_data.dtype.name:  # int8 uint8
            datatype = gdal.GDT_Byte
        elif 'int16' in im_data.dtype.name:  # int16 uint16~DN值/反射率
            datatype = gdal.GDT_UInt16
        elif 'int32' in im_data.dtype.name:  # int16 uint16~DN值/反射率
            datatype = gdal.GDT_UInt32
        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
            # 创建文件
        driver = gdal.GetDriverByName("GTiff")
        dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
        dataset.SetGeoTransform(im_GeoTransform)
        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

    # =============================================================================
    #     函数3 calndvi，输入波段B G R NIR数组，返回单波段数组
    # =============================================================================
    def calndvi(self, im_data):
        # print(im_data.shape)
        band_r = im_data[2, :, :]#band3 R
        band_nir = im_data[3, :, :]#band4 NIR
        ndvi = (band_nir - band_r) / (band_nir + band_r)
        # print(ndvi.dtype.name)
        return ndvi

if __name__ == "__main__":


    filename = r'D:\浏览器下载位置\GF1_WFV4_E116.1_N38.5_20230621_L1A0007353828.tiff'
    run=Grid()
    img_proj, im_GeoTransform, im_data=run.read_tiff(filename)
    ndvi=run.calndvi(im_data)
    print(ndvi)
    run.write_tiff('1.tiff', img_proj, im_GeoTransform, ndvi)

代码讲解

1、GDAL库的open函数可以获取.tiff格式的影像，将影像存储在dataset的数据集对象中，该对象包括了波段数、宽高、仿射矩阵、地图投影信息信息。

分别用GetRasterBand、RasterXSize、RasterYSize、GetGeoTransform()、GetProjection()来获取这些信息。

2、ReadAsArray() 方法接受四个参数：起始点的横坐标和纵坐标 ( 0, 0)，以及要读取的宽度和高度 (im_width, im_height)。这意味着它会从数据集中的 (0, 0) 点开始，读取宽度为 im_width，高度为 im_height 的数组数据。

返回的结果是一个包含数组数据的变量 im_data。你可以在后续的代码中使用 im_data 来处理或分析数据集中的数组信息。

3、通常情况下，使用del可以释放变量所占用的内存空间，或者删除列表中的某个元素，或者删除对象等。

4、后续的操作只需要对im_data进行就行。im_data[0, :, :]，im_data[1, :, :]，im_data[2, :, :]，im_data[3, :, :]分别获取4个波段的值

5、然后计算NDVI为（近红-红)/(近红+红)

6、创建一个Grid类的子对象run。调用Grid类的读数据函数，获取im_data，使用calndvi计算NDVI，最后调用写函数，记得名字要带后缀.tiff.

验证

看看计算的NDVI的tiff图像在Arcgis中运行怎么样