赞
踩
GDAL(Geospatial Data Abstraction Library)是一个在X/MIT许可协议下的开源栅格空间数据转换库。它利用抽象数据模型来表达所支持的各种文件格式。它还有一系列命令行工具来进行数据转换和处理。
Python的GDAL库作为栅格数据的处理转换库,其支持几百种栅格数据格式,如常见的TIFF、ENVI、HFA、HDF4等。因为遥感影像大部分都是栅格数据,所以GDAL库非常适合处理遥感影像、如光谱指数计算、波段合成、批量下载等。
因为博主的GDAL库的安装有些时候了,所以pip install 的方法能不能安装需要自己尝试。这里推荐大家下载对应包,使用本地安装。
1.打开GDAL库的下载链接,选择自己电脑Python对应的版本,别下载错了哦!
2.然后按住Win+R,输入cmd进入命令行。在cmd终端中,使用cd跳转至下载GDAL包的目录。如 cd /d G:\try\,G:\try\是你狭窄GDAL包的文件夹目录。然后在命令行里输入pip GDAL-3.4.3-cp38-cp38-win_amd64.whl,这是你下载的GDAL包的名称,每个人不一样,别傻乎乎地用我的!
3.可以在编译器里查看是否安装成功,博主使用的是PyCharm,在设置中可以看到是否安装成功。当然你也可以先编写程序,如果没安装成功,程序会报错=。=
1.导入计算所需库
numpy是数组处理的库,GDAL是我们安装的栅格数据处理库。遥感影像本质上是多维数组,如果这都不知道,建议先百度。
- import numpy as np
- from osgeo import gdal
2.查看栅格数据信息,如宽度、高度、波段数、投影等
当然这一步和我们的光谱指数计算没有太大的关系,但后续计算时需要这里面的相关参数,所以大家可以先了解一下。
- ds = gdal.Open(filepath) # 打开数据集dataset
- ds_width = ds.RasterXSize # 获取数据宽度
- ds_height = ds.RasterYSize # 获取数据高度
- ds_bands = ds.RasterCount # 获取波段数
- ds_geo = ds.GetGeoTransform() # 获取仿射地理变换参数
- ds_prj = ds.GetProjection() # 获取投影信息
- print("影像的宽度为:" + str(ds_width))
- print("影像的高度为:" + str(ds_height))
- print("仿射地理变换参数为:" + str(ds_geo))
- print("投影坐标系为:" + str(ds_prj))
3.以数组的形式打开红波段、近红外波段
每种卫星对应的红波段和近红外波段都不一样、博主这里red是3,nir是4。
- ds = gdal.Open(filepath) # 打开数据集dataset
- ds_width = ds.RasterXSize # 获取数据宽度
- ds_height = ds.RasterYSize # 获取数据高度
- ds_geo = ds.GetGeoTransform() # 获取仿射地理变换参数
- ds_prj = ds.GetProjection() # 获取投影信息
- array_red = ds.GetRasterBand(red).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float)
- # red是红波段对应的波段数
- array_nir = ds.GetRasterBand(nir).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float)
- # nir是近红波段对应的波段数
- # 以数组的形式读取红波段和近红外波段
4.计算NDVI(归一化植被指数)
其他光谱指数类似,指数更换red、nir对应的波段即可。
- b1 = array_nir-array_red
- b2 = array_nir+array_red
- NDVI_data = np.divide(b1, b2, out=np.zeros_like(b1), where=b2 != 0) # 计算NDVI
5.保存计算好的NDVI
这里需要先创建一个数据驱动,用于存储数组,类型可以替换(如可以保存为TIFF、ENVI、HFA等不同格式),这里我保存的TIFF格式,但需要注意的是GDAL中的类型命名不同,TIFF格式不能直接写TIFF,而是GTIFF,具体格式要求大家可以参考Python与开源GIS。
driver = gdal.GetDriverByName('GTiff') # 载入数据驱动,用于存储内存中的数组
随后我们需要创建一个空数组,给这个数组写入各类属性,如高度、宽度、波段数、投影等,一般我们都使用原始输入数据的属性。
- s_result = driver.Create(out_path, ds_width, ds_height, bands=1, eType=gdal.GDT_Float64)
- # 创建一个数组,宽高为原始尺寸
- ds_result.SetGeoTransform(ds_geo) # 导入仿射地理变换参数
- ds_result.SetProjection(ds_prj) # 导入投影信息
- ds_result.GetRasterBand(1).SetNoDataValue(9999) # 将无效值设为9999
最后,我们再将计算好的NDVI写入数组,删除缓存即可。注意缓存一定要删除,不然NDVI不会写入数组。
- ds_result.GetRasterBand(1).WriteArray(NDVI_data) # 将NDVI的计算结果写入数组
- del ds_result
- # 删除内存中的结果,否则结果不会写入图像中
- # -*- coding: utf-8 -*-
- """
- @Time : 2023/3/30 9:11
- @Auth : RS迷途小书童
- @File :NDVI计算.py
- @IDE :PyCharm
- """
- import numpy as np
- from osgeo import gdal
-
-
- def Get_data(filepath):
- ds = gdal.Open(filepath) # 打开数据集dataset
- ds_width = ds.RasterXSize # 获取数据宽度
- ds_height = ds.RasterYSize # 获取数据高度
- ds_bands = ds.RasterCount # 获取波段数
- ds_geo = ds.GetGeoTransform() # 获取仿射地理变换参数
- ds_prj = ds.GetProjection() # 获取投影信息
- print("影像的宽度为:" + str(ds_width))
- print("影像的高度为:" + str(ds_height))
- print("仿射地理变换参数为:" + str(ds_geo))
- print("投影坐标系为:" + str(ds_prj))
- # data = ds.ReadAsArray(0, 0, ds_width, ds_height) # 以数组的形式读取整个数据集
-
-
- def Write_NDVI(filepath, out_path, red, nir):
- print("-----正在进行归一化植被指数NDVI计算-----")
- ds = gdal.Open(filepath) # 打开数据集dataset
- ds_width = ds.RasterXSize # 获取数据宽度
- ds_height = ds.RasterYSize # 获取数据高度
- ds_geo = ds.GetGeoTransform() # 获取仿射地理变换参数
- ds_prj = ds.GetProjection() # 获取投影信息
- array_red = ds.GetRasterBand(red).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float)
- # red是红波段对应的波段数
- array_nir = ds.GetRasterBand(nir).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float)
- # nir是近红波段对应的波段数
- # 以数组的形式读取红波段和近红外波段
- b1 = array_nir-array_red
- b2 = array_nir+array_red
- NDVI_data = np.divide(b1, b2, out=np.zeros_like(b1), where=b2 != 0) # 计算NDVI
- driver = gdal.GetDriverByName('GTiff') # 载入数据驱动,用于存储内存中的数组
- ds_result = driver.Create(out_path, ds_width, ds_height, bands=1, eType=gdal.GDT_Float64)
- # 创建一个数组,宽高为原始尺寸
- ds_result.SetGeoTransform(ds_geo) # 导入仿射地理变换参数
- ds_result.SetProjection(ds_prj) # 导入投影信息
- ds_result.GetRasterBand(1).SetNoDataValue(9999) # 将无效值设为9999
- ds_result.GetRasterBand(1).WriteArray(NDVI_data) # 将NDVI的计算结果写入数组
- del ds_result
- # 删除内存中的结果,否则结果不会写入图像中
- print("计算完成")
-
-
- if __name__ == "__main__":
- file_path = r"G:/Shanghai_TreeClassicfication/Sen2Cor_Shanghai/Shanghai.dat"
- # 输入的栅格数据路径
- out_path1 = r"G:/GDAL_try/Shanghai-1.tif"
- # 导出的文件路径
- Get_data(file_path) # 执行函数,获取影像基本信息
- print("\n")
- print("--------------NDVI计算--------------")
- red1 = int(input("请输入红波段:"))
- nir1 = int(input("请输入近红外波段:"))
- Write_NDVI(file_path, out_path1, red1, nir1) # 执行函数,计算NDVI
使用Python去处理遥感影像可以大大减少我们的工作时间,提升工作效率。可能在一幅遥感影像处理时感觉还没有ENVI计算NDVI的速度快,但实际上如果你需要研究长时间序列的NDVI或者大量卫星影像,你就会发现编程给你带来的便利不是一点半点。针对上面分享的程序,只需加入os.listdir和for循环即可遍历文件夹中所有的栅格数据,成百上千幅影像在后台一会就处理好了。同时如果有需要,还可以直接将程序作为函数调用,用在完整的遥感影像预处理、定量分析中。本文章仅作为学习使用,如有侵权请联系作者删除。
如果大家在学习Python或者RS时有什么问题,可以随时留言交流!如果大家对批量处理有兴趣同样可以留言给博主,博主会分享相关代码以供学习!
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。