Python读取遥感影像并计算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)