一些简单的遥感中相关光谱指数计算

　　本节主要利用python 和 GDAL来计算一些常用的遥感指数，如，NDVI（归一化植被指数）、RVI（比值植被指数）、NDWI （归一化水体指数）、SAVI（土壤调节植被指数）。

1.读取波段数据

import numpy as np
from osgeo import gdal
in_ds = gdal.Open(r"D:\data\test.tif")  # 打开样本文件
xsize = in_ds.RasterXSize  # 获取行列数
ysize = in_ds.RasterYSize
bands = in_ds.RasterCount  # 获取波段数
geotransform = in_ds.GetGeoTransform() # 获取投影信息
projection = in_ds.GetProjectionRef()
block_data = in_ds.ReadAsArray(0,0,xsize,ysize).astype(np.float32)# 获取影像信息
B = block_data[0, :, :]
G = block_data[1,:, :]
R = block_data[2,:, :]
NIR = block_data[3,:, :]
    

2.表达式计算

#归一化植被指数
NDVI=(NIR-R)/(NIR+R)
#归一化水体指数
NDWI=(G-NIR)/(G+NIR)
#比值植被指数
RVI=NIR/R
#土壤调节植被指数
L=0.5 #L是随着植被密度变化的参数，取值范围从0-1，当植被覆盖度很高时为0，很低时为1。一般取0.5
SAVI=(NIR-R)*(1+L)/(NIR+R+L)

3.保存成果

#以NDVI为例
driver = gdal.GetDriverByName('GTiff')
out_dataset=driver.Create("NDVI.tif",xsize,ysize,1,gdal.GDT_Float32)
out_band1 = out_dataset.GetRasterBand(1)
out_band1.WriteArray(NDVI)
out_dataset.SetGeoTransform(geotransform)  # 写入仿射变换
out_dataset.SetProjection(projection)