简单的计算NDVI（以landsat 8为例)

# -*- coding: UTF-8 -*-
from osgeo import gdal

filename = 'C:\\Users\\pzt\\Desktop\\NDVI\\2023年NDVI1.tif'
dataset = gdal.Open(filename)
if dataset is None:
    print('could not open' + filename)

cols = dataset.RasterXSize  # 列数
rows = dataset.RasterYSize  # 行数  空间坐标 可以理解为行列号
print(cols, rows)
bands = dataset.RasterCount  # 波段数

# dataset.GetRasterBand():获取指定波段  dataset.ReadAsArray(xoffset,yoffset,1,1):根据像素位置获取像素值
band1 = dataset.GetRasterBand(1)  # 读取第1个波段
proj = dataset.GetProjection()  # 投影

geotrans = dataset.GetGeoTransform()  # 左上角x坐标、水平分辨率、旋转参数、左上角y坐标、旋转参数、竖直分辨率
print(geotrans)

data_all = dataset.ReadAsArray(0, 0, cols, rows)  # 读取从（0，0）开始，大小为（xsize,ysize）的矩阵
print(data_all)

# 计算 NDVI=（NIR-RED）/（NIR+RED）
# 提取红波段
band_red = dataset.GetRasterBand(4)
data_red = band_red.ReadAsArray(0, 0, cols, rows)  # 获取像素值
# 提取近红外波段
band_nir = dataset.GetRasterBand(5)
data_nir = band_nir.ReadAsArray(0, 0, cols, rows)

ndvi = (data_nir - data_red) / (data_nir + data_red)
# 创建栅格数据集
filename_out = 'C:\\Users\\pzt\\Desktop\\NDVI'
'''
GetDriverByName(driver_name)函数中的driver_name是文件类型的名称，
如TIFF/GeoTIFF文件的名称为GTiff
Erdas Imagine文件的名称为HFA
ENVI文件的名称为ENVI
详细的栅格文件名称可参考https://gdal.org/drivers/raster/index.html。
'''
driver = gdal.GetDriverByName('GTiff')  # 将GTiff图纸实体化为driver，括号内为文件类型

dataset_out = driver.Create(filename_out, cols, rows, 1, gdal.GDT_Float32)  # 这里的1指的是创建的栅格数据文件只能容纳一个波段
# 定义仿射矩阵和投影信息
dataset_out.SetGeoTransform(geotrans)
dataset_out.SetProjection(proj)
# 写入像素值
band_out = dataset_out.GetRasterBand(1)
band_out.WriteArray(ndvi)  # 把数组写入到栅格数据相应的波段中
# 清空缓存，关闭数据集
dataset.FlushCache()
dataset_out.FlushCache()  # 把内存中的数据写入到文件中
del dataset
del dataset_out