2023.7.6 汇报（自用

利用遥感图像计算植被覆盖度

1.计算NDVI（归一化植被指数）

from osgeo import gdal
 
filename='D:\\影像练习\\001.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(3)
data_red=band_red.ReadAsArray(0,0,cols,rows) #获取像素值
#提取近红外波段
band_nir=dataset.GetRasterBand(4)
data_nir=band_nir.ReadAsArray(0,0,cols,rows)
 
ndvi=(data_nir-data_red)/(data_nir+data_red)
#创建栅格数据集
filename_out='D:\\影像练习\\NDVI输出\\001_ndvi.tif'
'''
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

2.计算植被覆盖度（FVC）

采用像元二分模型来反演研究区域的植被覆盖度

 通常情况下：NDVIsoil取值为与NDVI累计百分比5%最接近的值；NDVIveg取值为与NDVI累计百分比95%最接近的值。

代码明天写，今天把数据中的红光和近红光光谱提取了出来，然后计算ndvi，并用arcmap验证代码计算的是否正确。