【PythonRS】植被显示增强(多光谱、正射、照片等)

        很多时候我们需要某个区域的正射图,虽然正射图一般都运用了匀色的算法,整体色彩比较均衡。但如果研究区内有大量的植被,这个时候植被突出显示就很有必要了。所以今天给大家分享一下使用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.8+NIR*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, bands=ds_bands+1, eType=gdal.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_bands+1).WriteArray(array_new_green)  # 将计算好的新绿波段写入del ds_result# 删除内存中的结果,否则结果不会写入图像中

三、RGB图像植被增强(正射影像)

1.算法核心

        由于RGB影像是没有近红外波段的,那么我们该怎么才能即增强植被的显示效果,又不改变其他地物的显示效果呢?这里我通过计算两种无需近红外波段就能计算的植被指数作为阈值计算的对象。很多植被指数我都试验过了,就这两种效果比较好:

1)绿叶指数 (GLI)

        GLI 指数在识别所有这些绿色和深色时非常敏感。负值将显示裸露的土壤、水体和人造基础设施。相反,正值就是类似于 NDVI 的颜色渐变来识别植被覆盖。如果你所在的区域有阴影、绿色潮湿的区域或裸露、潮湿和深色的土壤,应注意不要将植物覆盖分类错误。

GLI = ((绿 - 红) + (绿 - 蓝)) / ((2 * 绿) + 红 + 蓝)

2)红绿蓝植被指数 (RGBVI)

        利用可见光的三个波段。需要对值进行无监督的分类,以将植物与其他土地覆盖区分开来。

RGBVI = ((绿 * 绿) - (红 * 蓝)) / ((绿 * 绿) + (红 + 蓝))

        我们得到两种植被指数后,就需要建立一个合适的算法,去近似代替绿波段。这里我自己提出了两种算法,都是我实验之后效果比较好的:

Green = Green+Green*GLI

Green = Green*0.8+RGBVI*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(divide='ignore', invalid='ignore')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, bands=ds_bands+1, eType=gdal.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_red+array_blue)RGBVI = ((array_green*array_green)-(array_red*array_blue))/((array_green*array_green)+(array_red+array_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.8+RGBVI*0.2*255else:print('选择错误!')for i in range(1, ds_bands+1):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_bands+1).SetNoDataValue(0)ds_result.GetRasterBand(ds_bands+1).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, bands=ds_bands+1, eType=gdal.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_bands+1).WriteArray(array_new_green)  # 将计算好的新绿波段写入del ds_result# 删除内存中的结果,否则结果不会写入图像中def Image_Compress(path_image, path_out_image):""":param path_image: 输入需要压缩的影像路径:param path_out_image: 输出压缩后的影像路径:return: None"""ds = gdal.Open(path_image)# 打开影像数据driver = gdal.GetDriverByName('GTiff')# 创建输出的数据驱动driver.CreateCopy(path_out_image, ds, strict=1, options=["TILED=YES", "COMPRESS=PACKBITS", "BIGTIFF=YES"])# 设置压缩参数"""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(divide='ignore', invalid='ignore')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, bands=ds_bands+1, eType=gdal.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_red+array_blue)RGBVI = ((array_green*array_green)-(array_red*array_blue))/((array_green*array_green)+(array_red+array_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.8+RGBVI*0.2*255else:print('选择错误!')for i in range(1, ds_bands+1):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_bands+1).SetNoDataValue(0)ds_result.GetRasterBand(ds_bands+1).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, green=2, nir=4)  # 执行多光谱植被增强函数Enhance_Veg_RGB(fn_input.name, fn_output, red=1, green=2, blue=3)  # 执行RGB植被增强函数# Image_Compress(fn_output, fn_output)

五、效果图

        两图均是2%线性拉伸显示的效果。植被增强的效果很明显。

图1   原始图像

 图2   植被增强后

        需要注意的是,我为了保护原始文件,将新计算的绿波段加在了最后一个波段,所以要想查看增强效果,就用最后一个波段当成绿波段去真彩色显示!!!

        总的来说,本次博文分享的算法效果还不错。但代码并没有集成好,勉勉强强能用就行,感兴趣的话可以自己将它封装成程序。

        由于我代码中已给详细的解释,所以就不单独加以文字说明了。本文章主要是分享个人在学习Python过程中写过的一些代码。有些部分借鉴了前人以及官网的教程,如有侵权请联系作者删除,大家有问题可以随时留言交流,博主会及时回复。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/32367.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

【学习FreeRTOS】第2章——FreeRTOS基础知识

1.任务调度 1.1.任务调度简介 调度器就是使用相关的调度算法来决定当前需要执行的哪个任务FreeRTOS 一共支持三种任务调度方式: 抢占式调度:针对优先级不同的任务,每个任务都有一个优先级,优先级高的任务可以抢占优先级低的任务…

(3)原神角色数据分析-3

绘图类 在名为“WRITEPHOT.py”的文件中,定义如下绘图方式,则在主页面(app.py)文件中,可通过如下方式调用: from WRITEPHOTO import WriteScatter,WriteFunnel,WriteBarData,WritePie,WriteLineBar 代码如下: "…

【前端】CSS垂直居中的7种方法

文章目录 line-height绝对定位margin:autoflex绝对定位margin:负值定位transformvertical-align:middledisplay:table-cell思维导图 前文:【前端】CSS水平居中的6种方法_karshey的博客-CSDN博客 有很多相似的部分。 line-height 适用于单行的行内元素设置line-he…

N皇后算法的思路是什么?

N皇后问题是经典的回溯算法问题,其目标是在一个NxN的棋盘上放置N个皇后,使得它们彼此之间不能互相攻击,即任意两个皇后都不在同一行、同一列或同一斜线上。 算法的思路如下: 使用递归函数来解决问题。递归函数的参数通常包括当前…

模仿火星科技 基于cesium+角度测量+高度测量+可编辑

1. 创建提示窗: 启动Cesium应用,地图场景将打开,欢迎您进入编辑模式。 在屏幕的一角,一个友好的提示窗将呈现,随着您的操作,它会为您提供有用的信息和指导。 2. 绘制面积: 轻轻点击鼠标左键&a…

Spring之Aop切面---日志收集(环绕处理、前置处理方式)--使用/教程/实例

Spring之Aop切面---日志收集(环绕处理、前置处理方式)--使用/教程/实例 简介系统登录日志类LoginLogEntity .java 一、环绕处理方式1、自定义注解类LoginLogAop.class2、切面处理类LogoutLogAspect.java 二、前置处理方式:1、自定义注解类Log…

MySQL之 show profile 相关总结

MySQL之 show profile 相关总结 MySQL官网show profile介绍:https://dev.mysql.com/doc/refman/8.0/en/show-profile.html 1. 简介 show profile 和 show profiles 命令用于展示SQL语句的资源使用情况,包括CPU的使用,CPU上下文切换&#xf…

【实战】 九、深入React 状态管理与Redux机制(五) —— React17+React Hook+TS4 最佳实践,仿 Jira 企业级项目(二十)

文章目录 一、项目起航:项目初始化与配置二、React 与 Hook 应用:实现项目列表三、TS 应用:JS神助攻 - 强类型四、JWT、用户认证与异步请求五、CSS 其实很简单 - 用 CSS-in-JS 添加样式六、用户体验优化 - 加载中和错误状态处理七、Hook&…

Docker 数据管理

文章目录 前言1、Dcoker 文件体系2、volume挂载案例2.1、挂载运行一个容器实例方法1方法2 3、volumes-from 案例4、备份/恢复数据卷5、删除数据卷 前言 为什么要有数据管理? 因为: Docker 是不提供持久化的 ,容器是不稳定的;一个…

mac ssh连接另一台window虚拟机vm

vmware配置端口映射 编辑(E) > 虚拟网络编辑器(N)... > NAT设置(S)... window防火墙,入站规则添加5555端口 控制面板 > 系统和安全 > Windows 防火墙>高级设置>入站规则>新建规则... tips windows查看端口命令:netstat -ano | f…

【从零学习python 】14.Python字符串常见操作(二)

文章目录 count替换内容分隔splitrsplitsplitlinespartitionrpartition 修改大小写capitalizetitlelowerupper 空格处理1. ljust2. rjust3. center4. lstrip5. rstrip6. strip 字符串拼接字符串运算符进阶案例 count 返回 str 在 start 和 end 之间在 mystr 里面出现的次数。 …

EthereumArbitrumOptimism节点搭建

Ethereum节点搭建 Ethereum节点搭建可基于源码或基于docker。这里介绍基于docker的方式搭建。 拉取镜像 docker pull ethereum/client-go 建议将镜像源配置为阿里云的镜像源 启动容器 docker run -itd -p 8545:8545 -p 30303:30303 --name eth -v /data/eth:/root/.ethere…

java-IDEA MAVEN查看依赖树,解决jar包重复和冲突

如果这里面的依赖关系有红线,就说明有包冲突,一般都是版本不一致,可以在idea里下一个插件Maven Helper,点击install并重启IDEA 打开pom.xml文件,在下方会出现Dependency Analyzer,选择它会出现重复依赖列表,选择对应的依赖,右键红…

BI技巧丨利用Index计算半累计

在实际的业务场景中,特别是财务模块和库存管理模块,经常需要我们针对每个月的期初期末进行相关指标计算,这也是我们之前曾经提到的Calculate基础应用——半累计计算。 现在我们也可以通过微软新推出的Index开窗函数来解决这一问题。 INDEX函…

Mapbox加载天地图CGCS2000矢量瓦片地图

1.背景 最近在做天地图的项目,要基于MapBox添加CGCS2000矢量切片数据,但是 Mapbox 只支持web 墨卡托(3857)坐标系的数据。Github有专业用户修改了mapbox-gl的相关代码,支持CGCS2000的切片数据加载,并且修改…

dji uav建图导航系列(一)建图

文章目录 1、uav + rplidir雷达1.2、思岚激光雷达1.3、dji uav的launch文件2、cartographer激光建图2.1、启动文件2.2、config修改2.3、建图过程3、融合odom+laser建图1、uav + rplidir雷达 思岚激光雷达frame为base_laser_link, 无人机frame为base_footprint。 文件uav_lid…

【Python】同步、异步、堵塞、非堵塞、回调

原文作者:我辈李想 版权声明:文章原创,转载时请务必加上原文超链接、作者信息和本声明。 文章目录 一、Python中的同步异步二、Python中的堵塞非堵塞三、Python中的回调四、异步编程 一、Python中的同步异步 在 Python 中,同步和…

【OS】操作系统中,页框/页帧见解

文章目录 页框大小的取决因素页框大小的查询不同大小的页框,适用于哪些场景? 页框大小的取决因素 页框(Page Frame)是指操作系统中用于存储页面的物理内存块,其大小通常由操作系统的设计和硬件架构决定。在现代计算机…

手动实现 Spring 底层机制 实现任务阶段一编写自己 Spring 容器-准备篇【2】

😀前言 手动实现 Spring 底层机制的第2篇 实现了任务阶段一编写自己 Spring 容器-准备篇【2】 🏠个人主页:尘觉主页 🧑个人简介:大家好,我是尘觉,希望我的文章可以帮助到大家,您的…

docker 删除镜像文件

docker 容器里面太多镜像,D盘满了 四 查看和移除镜像 1 查看镜像 docker images 2 移除镜像命令 docker rmi 镜像名称 # 只输入前四位即可 五 实际有效操作 清除所有不使用的资源 docker system prune 这个命令将会删除所有不使用的镜像、容器和数据卷等资…