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的替换者)