基于python的多光谱影像植被指数计算

利用GDAL库对多光谱影像进行分析，光谱5波段，顺序蓝-绿-红-红边-近红外

"""
author: Shuai-jie Shen 沈帅杰
CSDN: https://blog.csdn.net/weixin_45452300
公众号: AgBioIT
"""
from osgeo import gdal, osr


def export_tif(out_tif_name, var_lat, var_lon, NDVI):
    # 判断数组维度
    if len(NDVI.shape) > 2:
        im_bands, im_height, im_width = NDVI.shape
    else:
        im_bands, (im_height, im_width) = 1, NDVI.shape
    LonMin, LatMax, LonMax, LatMin = [min(var_lon), max(var_lat), max(var_lon), min(var_lat)]
    # 分辨率计算
    Lon_Res = (LonMax - LonMin) / (float(im_width) - 1)
    Lat_Res = (LatMax - LatMin) / (float(im_height) - 1)
    driver = gdal.GetDriverByName('GTiff')
    out_tif = driver.Create(out_tif_name, im_width, im_height, im_bands, gdal.GDT_Float32)  # 创建框架

    # 设置影像的显示范围
    # Lat_Res一定要是-的
    geotransform = (LonMin, Lon_Res, 0, LatMax, 0, -Lat_Res)
    out_tif.SetGeoTransform(geotransform)

    # 获取地理坐标系统信息，用于选取需要的地理坐标系统
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)  # 定义输出的坐标系为"WGS 84"，AUTHORITY["EPSG","4326"]
    out_tif.SetProjection(srs.ExportToWkt())  # 给新建图层赋予投影信息

    # 数据写出
    if im_bands == 1:
        out_tif.GetRasterBand(1).WriteArray(NDVI)
    else:
        for bands in range(im_bands):
            out_tif.GetRasterBand(bands + 1).WriteArray(NDVI[bands])
    del out_tif


filePath = r'E:\orthomosaic_rotation_exp20211013_ref.tif'  # 五波段影像
dataset = gdal.Open(filePath)  # 读取数据
adfGeoTransform = dataset.GetGeoTransform()
nXSize = dataset.RasterXSize  # 列数
nYSize = dataset.RasterYSize  # 行数
im_data = dataset.ReadAsArray(0, 0, nXSize, nYSize)  # 读取数据
var_lat = []  # 纬度
var_lon = []  # 经度
for j in range(nYSize):
    lat = adfGeoTransform[3] + j * adfGeoTransform[5]
    var_lat.append(lat)
for i in range(nXSize):
    lon = adfGeoTransform[0] + i * adfGeoTransform[1]
    var_lon.append(lon)
# 不同的指标结合成三维以上的array也可以输出为多通道光谱数据
GNDVI = (im_data[4] - im_data[1])/(im_data[4] + im_data[1])
NDRE = (im_data[4] - im_data[3])/(im_data[4] + im_data[3])
LCI = (im_data[4] - im_data[3])/(im_data[4] + im_data[2])
OSAVI = (im_data[4] - im_data[2])/(im_data[4] + im_data[2]+0.16)
out_tif_name = r'E:\GNDVI_py.tif'
out_tif_name1 = r'E:\NDRE_py.tif'
out_tif_name2 = r'E:\LCI_py.tif'
out_tif_name3 = r'E:\OSAVI_py.tif'
export_tif(out_tif_name, var_lat, var_lon, GNDVI)
export_tif(out_tif_name1, var_lat, var_lon, NDRE)
export_tif(out_tif_name2, var_lat, var_lon, LCI)
export_tif(out_tif_name3, var_lat, var_lon, OSAVI)