用阳寿做交易的网站,北京工信部网站,湖南常德石门县,中国中建设计集团有限公司很多时候我们需要某个区域的正射图#xff0c;虽然正射图一般都运用了匀色的算法#xff0c;整体色彩比较均衡。但如果研究区内有大量的植被#xff0c;这个时候植被突出显示就很有必要了。所以今天给大家分享一下使用Python对多光谱、正射影像进行植被显示增强的算法。 一、… 很多时候我们需要某个区域的正射图虽然正射图一般都运用了匀色的算法整体色彩比较均衡。但如果研究区内有大量的植被这个时候植被突出显示就很有必要了。所以今天给大家分享一下使用Python对多光谱、正射影像进行植被显示增强的算法。 一、加载需要的库 这里的主要库是numpy用来读取图片的数组gdal库用来读取多光谱数据。
import numpy as np
from osgeo import gdal
import tkinter.filedialog
# 创建窗口打开图片模块
二、多光谱数据增强
1.获取影像的基本信息 这一步没啥用可以不放进代码。主要让你先了解一下自己的数据。
def Get_data(filepath)::param filepath: 输入数据路径:return: 输出影像的基本信息ds gdal.Open(filepath) # 打开数据集datasetds_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_bands))print(仿射地理变换参数为 str(ds_geo))print(投影坐标系为 str(ds_prj))# data ds.ReadAsArray(0, 0, ds_width, ds_height) # 以数组的形式读取整个数据集
2.算法核心 由于多光谱是包含近红外波段的同时呢我们知道植被对于近红外波段有着很强的反射其栅格值也会很大。所以这里我们使用近红外波段近似代替绿波段这样就可以达到植被显示增强的目的了。但如果只是单纯的替代那么其他对绿波段有反射的地物的颜色就会失真你就会发现增强后整个图变得非常奇怪。所以这里我们设定一个阈值将绿波段和近红外波段综合起来再去代替绿波段。
Green Green*0.8NIR*0.2
3.代码实现 这里我就不过多的解释了因为里面的一些函数我之前发布的博文里都已经解释过了同时代码我也做了详细的注释所以大家直接看吧
def Enhance_Veg(filepath, out_path, green, nir)::param filepath: 输入需要增强的图片路径:param out_path: 输入保存文件的路径:param green: 输入绿波段:param nir: 输入近红外波段:return:ds gdal.Open(filepath) # 打开数据集datasetds_width ds.RasterXSize # 获取数据宽度ds_height ds.RasterYSize # 获取数据高度ds_geo ds.GetGeoTransform() # 获取仿射地理变换参数ds_prj ds.GetProjection() # 获取投影信息ds_bands ds.RasterCount # 获取波段数driver gdal.GetDriverByName(GTiff) # 载入数据驱动用于存储内存中的数组ds_result driver.Create(out_path, ds_width, ds_height, bandsds_bands1, eTypegdal.GDT_Float64)# 创建一个数组宽高为原始尺寸ds_result.SetGeoTransform(ds_geo) # 导入仿射地理变换参数ds_result.SetProjection(ds_prj) # 导入投影信息array_green ds.GetRasterBand(green).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float64)array_nir ds.GetRasterBand(nir).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float64)array_new_green array_green * 0.8 array_nir * 0.2for i in range(1, ds_bands):array_band ds.GetRasterBand(i).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float64)ds_result.GetRasterBand(i).SetNoDataValue(0) # 将无效值设为0ds_result.GetRasterBand(i).WriteArray(array_band) # 将每个波段写入新的文件中ds_result.GetRasterBand(ds_bands1).WriteArray(array_new_green) # 将计算好的新绿波段写入del ds_result# 删除内存中的结果否则结果不会写入图像中
三、RGB图像植被增强正射影像
1.算法核心 由于RGB影像是没有近红外波段的那么我们该怎么才能即增强植被的显示效果又不改变其他地物的显示效果呢这里我通过计算两种无需近红外波段就能计算的植被指数作为阈值计算的对象。很多植被指数我都试验过了就这两种效果比较好
1绿叶指数 (GLI) GLI 指数在识别所有这些绿色和深色时非常敏感。负值将显示裸露的土壤、水体和人造基础设施。相反正值就是类似于 NDVI 的颜色渐变来识别植被覆盖。如果你所在的区域有阴影、绿色潮湿的区域或裸露、潮湿和深色的土壤应注意不要将植物覆盖分类错误。
GLI ((绿 - 红) (绿 - 蓝)) / ((2 * 绿) 红 蓝)
2红绿蓝植被指数 (RGBVI) 利用可见光的三个波段。需要对值进行无监督的分类以将植物与其他土地覆盖区分开来。
RGBVI ((绿 * 绿) - (红 * 蓝)) / ((绿 * 绿) (红 蓝)) 我们得到两种植被指数后就需要建立一个合适的算法去近似代替绿波段。这里我自己提出了两种算法都是我实验之后效果比较好的
Green GreenGreen*GLI
Green Green*0.8RGBVI*255*0.2
2.代码实现 这个函数在运行时会询问你需要使用哪种算法再程序询问时输入1或者2回车即可。
def Enhance_Veg_RGB(filepath, out_path, red, green, blue)::param filepath: 输入需要增强的图片路径:param out_path: 输入保存文件的路径:param red: 输入红波段:param green: 输入绿波段:param blue: 输入蓝波段:return:np.seterr(divideignore, invalidignore)ds gdal.Open(filepath) # 打开数据集datasetds_width ds.RasterXSize # 获取数据宽度ds_height ds.RasterYSize # 获取数据高度ds_geo ds.GetGeoTransform() # 获取仿射地理变换参数ds_prj ds.GetProjection() # 获取投影信息ds_bands ds.RasterCount # 获取波段数driver gdal.GetDriverByName(GTiff) # 载入数据驱动用于存储内存中的数组ds_result driver.Create(out_path, ds_width, ds_height, bandsds_bands1, eTypegdal.GDT_Float64)# 创建一个数组宽高为原始尺寸ds_result.SetGeoTransform(ds_geo) # 导入仿射地理变换参数ds_result.SetProjection(ds_prj) # 导入投影信息array_red ds.GetRasterBand(red).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float64)array_green ds.GetRasterBand(green).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float64)array_blue ds.GetRasterBand(blue).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float64)GLI ((array_green-array_red) (array_green-array_blue))/((array_green*2)array_redarray_blue)RGBVI ((array_green*array_green)-(array_red*array_blue))/((array_green*array_green)(array_redarray_blue))array_new_green 0chance int(input(请选择增强算法1:GLI, 2:RGBVI\n))if chance 1:array_new_green array_green array_green*GLIelif chance 2:array_new_green array_green*0.8RGBVI*0.2*255else:print(选择错误!)for i in range(1, ds_bands1):array_band ds.GetRasterBand(i).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float64)ds_result.GetRasterBand(i).SetNoDataValue(0) # 将无效值设为0ds_result.GetRasterBand(i).WriteArray(array_band) # 将每个波段写入新的文件中ds_result.GetRasterBand(ds_bands1).SetNoDataValue(0)ds_result.GetRasterBand(ds_bands1).WriteArray(array_new_green) # 将计算好的新绿波段写入del ds_result# 删除内存中的结果否则结果不会写入图像中
四、完整代码 我这里使用了tkinter.filedialog调用了GUI窗口以便于更方便地选取需要增强的图片。本来还想封装成程序的但我太懒了。本来还想再出个ENVI实现的教程的但我太懒了。所以现在还需要手动输入波段索引同时还需要自己确定使用多光谱还是RGB增强函数。 我这里还写了一个压缩函数可以不要但不能没有。代码中已经注释掉了不会用就不用。
# -*- coding: utf-8 -*-Time 2023/8/10 17:45
Auth RS迷途小书童
File Vegetation Enhancement.py
IDE PyCharm
Purpose 影像中植被增强显示import numpy as np
from osgeo import gdal
import tkinter.filedialog
# 创建窗口打开图片模块def Get_data(filepath)::param filepath: 输入数据路径:return: 输出影像的基本信息ds gdal.Open(filepath) # 打开数据集datasetds_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_bands))print(仿射地理变换参数为 str(ds_geo))print(投影坐标系为 str(ds_prj))# data ds.ReadAsArray(0, 0, ds_width, ds_height) # 以数组的形式读取整个数据集def Enhance_Veg(filepath, out_path, green, nir)::param filepath: 输入需要增强的图片路径:param out_path: 输入保存文件的路径:param green: 输入绿波段:param nir: 输入近红外波段:return:ds gdal.Open(filepath) # 打开数据集datasetds_width ds.RasterXSize # 获取数据宽度ds_height ds.RasterYSize # 获取数据高度ds_geo ds.GetGeoTransform() # 获取仿射地理变换参数ds_prj ds.GetProjection() # 获取投影信息ds_bands ds.RasterCount # 获取波段数driver gdal.GetDriverByName(GTiff) # 载入数据驱动用于存储内存中的数组ds_result driver.Create(out_path, ds_width, ds_height, bandsds_bands1, eTypegdal.GDT_Float64)# 创建一个数组宽高为原始尺寸ds_result.SetGeoTransform(ds_geo) # 导入仿射地理变换参数ds_result.SetProjection(ds_prj) # 导入投影信息array_green ds.GetRasterBand(green).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float64)array_nir ds.GetRasterBand(nir).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float64)array_new_green array_green * 0.8 array_nir * 0.2for i in range(1, ds_bands):array_band ds.GetRasterBand(i).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float64)ds_result.GetRasterBand(i).SetNoDataValue(0) # 将无效值设为0ds_result.GetRasterBand(i).WriteArray(array_band) # 将每个波段写入新的文件中ds_result.GetRasterBand(ds_bands1).WriteArray(array_new_green) # 将计算好的新绿波段写入del ds_result# 删除内存中的结果否则结果不会写入图像中def Image_Compress(path_image, path_out_image)::param path_image: 输入需要压缩的影像路径:param path_out_image: 输出压缩后的影像路径:return: Noneds gdal.Open(path_image)# 打开影像数据driver gdal.GetDriverByName(GTiff)# 创建输出的数据驱动driver.CreateCopy(path_out_image, ds, strict1, options[TILEDYES, COMPRESSPACKBITS, BIGTIFFYES])# 设置压缩参数PACKBITS连续字节压缩快速无损压缩LZW所有信息全部保留可逆以某一数值代替字符串快速无损压缩del dsdef Enhance_Veg_RGB(filepath, out_path, red, green, blue)::param filepath: 输入需要增强的图片路径:param out_path: 输入保存文件的路径:param red: 输入红波段:param green: 输入绿波段:param blue: 输入蓝波段:return:np.seterr(divideignore, invalidignore)ds gdal.Open(filepath) # 打开数据集datasetds_width ds.RasterXSize # 获取数据宽度ds_height ds.RasterYSize # 获取数据高度ds_geo ds.GetGeoTransform() # 获取仿射地理变换参数ds_prj ds.GetProjection() # 获取投影信息ds_bands ds.RasterCount # 获取波段数driver gdal.GetDriverByName(GTiff) # 载入数据驱动用于存储内存中的数组ds_result driver.Create(out_path, ds_width, ds_height, bandsds_bands1, eTypegdal.GDT_Float64)# 创建一个数组宽高为原始尺寸ds_result.SetGeoTransform(ds_geo) # 导入仿射地理变换参数ds_result.SetProjection(ds_prj) # 导入投影信息array_red ds.GetRasterBand(red).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float64)array_green ds.GetRasterBand(green).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float64)array_blue ds.GetRasterBand(blue).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float64)GLI ((array_green-array_red) (array_green-array_blue))/((array_green*2)array_redarray_blue)RGBVI ((array_green*array_green)-(array_red*array_blue))/((array_green*array_green)(array_redarray_blue))array_new_green 0chance int(input(请选择增强算法1:GLI, 2:RGBVI\n))if chance 1:array_new_green array_green array_green*GLIelif chance 2:array_new_green array_green*0.8RGBVI*0.2*255else:print(选择错误!)for i in range(1, ds_bands1):array_band ds.GetRasterBand(i).ReadAsArray(0, 0, ds_width, ds_height).astype(np.float64)ds_result.GetRasterBand(i).SetNoDataValue(0) # 将无效值设为0ds_result.GetRasterBand(i).WriteArray(array_band) # 将每个波段写入新的文件中ds_result.GetRasterBand(ds_bands1).SetNoDataValue(0)ds_result.GetRasterBand(ds_bands1).WriteArray(array_new_green) # 将计算好的新绿波段写入del ds_result# 删除内存中的结果否则结果不会写入图像中if __name__ __main__:fn_input open(tkinter.filedialog.askopenfilename(title选择图片, filetypes[(所有文件, .*), (JPG, .jpg), (JPG, .jpeg),(TIFF, .tif), (DAT, .dat), (PNG, .png)]), rb)# 输入的栅格数据路径fn_output tkinter.filedialog.asksaveasfilename(defaultextension.tif)# 导出的文件路径Get_data(fn_input.name)# Enhance_Veg(fn_input.name, fn_output, green2, nir4) # 执行多光谱植被增强函数Enhance_Veg_RGB(fn_input.name, fn_output, red1, green2, blue3) # 执行RGB植被增强函数# Image_Compress(fn_output, fn_output)
五、效果图 两图均是2%线性拉伸显示的效果。植被增强的效果很明显。 图1 原始图像 图2 植被增强后 需要注意的是我为了保护原始文件将新计算的绿波段加在了最后一个波段所以要想查看增强效果就用最后一个波段当成绿波段去真彩色显示 总的来说本次博文分享的算法效果还不错。但代码并没有集成好勉勉强强能用就行感兴趣的话可以自己将它封装成程序。 由于我代码中已给详细的解释所以就不单独加以文字说明了。本文章主要是分享个人在学习Python过程中写过的一些代码。有些部分借鉴了前人以及官网的教程如有侵权请联系作者删除大家有问题可以随时留言交流博主会及时回复。
