python 影像读取和计算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

geopandas —— 矢量数据读取与基本操作

rasterio       —— 栅格数据读写与基本操作

matplotlib   —— 基本绘图控件

cartopy       —— 地理空间绘图的基础、关键库（Basemap的替换者）