基于Python 中创建 Sentinel-2 RGB 合成图像

一、前言

下面的python代码将带您了解如何从原始 Sentinel-2 图像创建 RGB 合成图像的过程。

免费注册后,可以从 Open Access Hub 下载原始图像。 请注意,激活您的帐户可能需要 24 小时!

二、准备工作

(1)导入必要的库


import matplotlib.pyplot as plt
import numpy as np
import rasterio
import os
from osgeo import gdal

加载库后,我们必须定义输入和输出文件夹。 输入文件夹是您必须放置下载的 Sentinel-2 图像的地方。 输出文件夹是我们的脚本将保存 RGB 合成图像的地方。

(2)还必须为每个波段指定输入文件名。


fp_in='input/'
fp_out='output/'fn_blue='Tihany_T33TYN_A021798_20210509T094028_B02'
fn_green='Tihany_T33TYN_A021798_20210509T094028_B03'
fn_red='Tihany_T33TYN_A021798_20210509T094028_B04'

要使用 Sentinel-2 图像,首先我们必须将它们从 '.jp2' 文件格式转换为 '.tif' 文件格式,因为 Rasterio 库只能处理后者。


bandList = [band for band in os.listdir(fp_in) if band[-4:]=='.jp2']
for band in bandList:in_image = gdal.Open(fp_in+band)driver = gdal.GetDriverByName("GTiff")fp_tif = fp_in+band[:-4]+'.tif'out_image = driver.CreateCopy(fp_tif, in_image, 0)in_image = Noneout_image = None   

三、创建原始 RGB 合成

让我们为每个转换后的 Sentinel-2 图像定义文件路径,并使用 Rasterio 打开它们。


band_02=rasterio.open(fp_in+fn_blue+'.tif')
band_03=rasterio.open(fp_in+fn_green+'.tif')
band_04=rasterio.open(fp_in+fn_red+'.tif')

现在我们必须读入打开的文件。


red = band_04.read(1)
green = band_03.read(1)
blue = band_02.read(1)

(1) 如果我们查看红色波段栅格,我们将看到以下内容:


plt.imshow(red)

蓝色图像基本上是一个强度图,其中每个像素代表 Sentinel-2 传感器在红色波段中捕获的反射光量。 较亮的像素(较高的值)代表更多的红色内容,而较暗的像素(较低的值)代表较少的红色内容。

我们可以使用“cmap”命令更改蓝色表示。 在下面的示例中,我选择了“Reds”表示。

(请注意,还有许多其他选项。有关更多详细信息,请查看 Matplotlib 文档)。


plt.imshow(red, cmap='Reds')

现在让我们看看红色、绿色和蓝色通道图像是什么样子的。


fig = plt.figure(figsize=(20,6))
ax1 = fig.add_subplot(1,3,1)
ax1.imshow(red, cmap='Reds')
ax1 = fig.add_subplot(1,3,2)
ax1.imshow(green, cmap='Greens')
ax1 = fig.add_subplot(1,3,3)
ax1.imshow(blue, cmap='Blues')

您可以使用 shape 命令获取红色带图像的大小,如下所示。

如您所见,此图像是一个具有 582 行和 981 列的二维数组。

要制作 RGB 合成图,我们必须使用 np.dstack 命令将红色、绿色和蓝色波段图像堆叠在一起成为一个图像。

如果我们在新创建的 RGB 合成上再次调用形状命令,我们将看到,现在我们得到了一个包含红色、绿色和蓝色通道的 3D 数组。


rgb_composite_raw= np.dstack((red, green, blue))
rgb_composite_raw.shape

现在让我们看一下 RGB 图像...


plt.imshow(rgb_composite_raw)

这不是我们要找的,对吧? 问题的根源在于大多数图像的像素值范围为 0-255 或 0-1。 如果我们查看红色带的最大像素值,我们会得到超过 255。

(2) 我们现在能做什么? 此问题的解决方案是对 0..1 之间的所有像素值进行归一化。


def normalize(band):band_min, band_max = (band.min(), band.max())return ((band-band_min)/((band_max - band_min)))red_n = normalize(red)
green_n = normalize(green)
blue_n = normalize(blue)

在对我们的图像进行归一化处理后,一个波段的最大值和最小值应为 0 和 1。

现在让我们再次进行 RGB 堆栈,看看我们得到了什么结果。


rgb_composite_n= np.dstack((red_n, green_n, blue_n))
plt.imshow(rgb_composite_n)

最后我们可以看到我们感兴趣的区域,但是颜色似乎不太真实,整个图像有点暗。

四、基本图像处理技术

(1)波段运算

为了解决这个问题,我们必须先使每个波段变亮,然后将它们归一化并进行叠加。从数学的角度来看,增亮函数将每个像素值与“alpha”相乘,并在必要时添加“beta”值。如果完成此操作,我们必须将结果像素值裁剪在 0..255 之间。


def brighten(band):alpha=0.13beta=0return np.clip(alpha*band+beta, 0,255)red_b=brighten(red)
blue_b=brighten(blue)
green_b=brighten(green)red_bn = normalize(red_b)
green_bn = normalize(green_b)
blue_bn = normalize(blue_b)

现在让我们看一下对波段进行增亮和标准化后的新 RGB 合成图。


rgb_composite_bn= np.dstack((red_bn, green_bn, blue_bn))
plt.imshow(rgb_composite_bn)

现在我们的图像看起来非常逼真。 请注意,此图像并不代表真实的反射率值。

另一种图像处理技术是伽马校正。 它背后的数学原理是我们采用每个像素的强度值并将其提高到 (1/gamma) 的幂,其中 gamma 值由我们指定。

让我们使用我们的原始图像,进行伽马校正和归一化。


def gammacorr(band):gamma=2return np.power(band, 1/gamma)red_g=gammacorr(red)
blue_g=gammacorr(blue)
green_g=gammacorr(green)red_gn = normalize(red_g)
green_gn = normalize(green_g)
blue_gn = normalize(blue_g)

现在让我们看看结果。


rgb_composite_gn= np.dstack((red_gn, green_gn, blue_gn))
plt.imshow(rgb_composite_gn)

(2)将图像保存为PNG

此时大多数人想做的是将图像保存到文件中。这可以通过以下代码行来完成,将亮化和规范化的图像保存到 PNG 文件中。

请注意,给出了一种插值方法来平滑图像,并且还可以控制 dpi 值。有关详细信息,请访问 Pyplot 文档。


rgb_plot=plt.imshow(rgb_composite_bn, interpolation='lanczos')
plt.axis('off')
plt.savefig(fp_out+'tihany_rgb_composite.png',dpi=200,bbox_inches='tight')
plt.close('all')

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

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

相关文章

selenium的基础语法

📑打牌 : da pai ge的个人主页 🌤️个人专栏 : da pai ge的博客专栏 ☁️山水速疾来去易,襄樊镇固永难开 ☁️定位页面的元素 参数:抽象类By里…

springboot 开启和关闭kafka消费

关闭kafka自动消费 配置自定义容器工厂 import org.springframework.beans.factory.annotation.Autowired; import org.springframework.context.annotation.Bean; import org.springframework.context.annotation.Configuration; import org.springframework.kafka.config.C…

【从删库到跑路 | MySQL总结篇】数据库基础(增删改查的基本操作)

个人主页:兜里有颗棉花糖 欢迎 点赞👍 收藏✨ 留言✉ 加关注💓本文由 兜里有颗棉花糖 原创 收录于专栏【MySQL学习专栏】🎈 本专栏旨在分享学习MySQL的一点学习心得,欢迎大家在评论区讨论💌 重点放前面&am…

Android frameworks 开发总结之八

Quick Settings增加一項 XXX device要求在quick settings中增加一項touch panel. 在/frameworks/base/packages/SystemUI/res/values/config.xml文件中的quick_settings_tiles_default string 中增加touch panel。並在String resource文件中增加顯示的title <!-- The def…

网络数据结构skb_buff原理

skb_buff基本原理 内核中sk_buff结构体在各层协议之间传输不是用拷贝sk_buff结构体&#xff0c;而是通过增加协议头和移动指针来操作的。如果是从L4传输到L2&#xff0c;则是通过往sk_buff结构体中增加该层协议头来操作&#xff1b;如果是从L4到L2&#xff0c;则是通过移动sk_…

Kafka(五)消费者回调 +定时重试 + 理解Rebalance

文章目录 消费者回调如何抽象callBack消息&#xff1f;为什么要设置serverId?如何消费callBack消息&#xff1f; 定时重试消息失败表的设计重试逻辑设计 理解Rabalance通过日志来理解rebalance 参考资料结语示例源码仓库 消费者回调 有些邮件发送成功之后&#xff0c;需要执行…

【Linux】fork()

文章目录 一、fork()是什么&#xff1f;二、fork()干了什么&#xff1f;三、fork()怎么用&#xff1f; 一、fork()是什么&#xff1f; fork()函数其实是在Linux系统中用于创建一个新的进程。让我们看看Linux中是怎么描述的&#xff1f;运行man fork。 RETURN VALUE On success…

php站点伪静态配置(Apache+Linux)

404报错&#xff1a; 404 Not Found nginx/1.15.11 问题解决&#xff1a; 1、Linux location / { if (!-e $request_filename) { rewrite ^(.*)$ /index.php?s/$1 last; } } 2、Apache <IfModule mod_rewrite.c> RewriteEngine on RewriteBase / RewriteCond %{REQU…

英特尔和 ARM 将合作开发移动芯片技术,如何看待双方合作?

英特尔和 ARM 将合作开发移动芯片技术&#xff0c;如何看待双方合作&#xff1f; 最近市场传出Arm要自产芯片&#xff0c;供智能手机与笔电等使用后&#xff0c;外媒指Arm自产芯片将由英特尔晶圆代工部门打造&#xff0c;变成英特尔晶圆代工客户。将采用英特尔18A工艺&#xff…

利用Nginx与php处理方式不同绕过Nginx_host实现SQL注入

目录 首先需要搭建环境 nginxphpmysql环境&#xff1a; 搭建网站 FILTER_VALIDATE_EMAIL 绕过 方法1&#xff1a;冒号号分割host字段 方法2&#xff1a;冒号号分割host字段 方法3&#xff1a;SNI扩展绕过 首先需要搭建环境 nginxphpmysql环境&#xff1a; php安装包&a…

深入了解Spring Cloud中的分布式事务解决方案

引言 介绍分布式系统中事务管理的重要性&#xff0c;以及在云计算环境下分布式事务所面临的挑战。 传统事务和分布式事务 解释本地事务与分布式事务的区别&#xff0c;以及为什么在分布式环境中需要特殊的事务管理机制。 分布式事务的挑战 探讨在分布式系统中实现事务一致性所…

vite和webpack的区别和练习

Vite和Webpack都是现代化的前端构建工具&#xff0c;但它们之间存在一些区别&#xff1a; 构建性能&#xff1a;Vite使用ES Modules提高了构建性能&#xff0c;可以在构建时只构建需要的部分&#xff0c;而Webpack则需要在构建时处理整个应用程序。 开发体验&#xff1a;Vite具…

vue一个页面左边是el-table表格 当点击每条数据时可以在右边界面编辑表格参数,右边保存更新左边表格数据

实现思路&#xff1a; 1.点击当前行通过row拿到当前行数据。 2.将当前行数据传给子组件。 3.子组件监听父组件传过来的数据并映射在界面。 4.点击保存将修改的值传给父组件更新表格。 5.父组件收到修改过后的值&#xff0c;可以通过字段判断比如id&#xff0c;通过 findIn…

VR Interaction Framework2.0使用

1 按键 &#xff0c;比如按压下手柄的B键 if (InputBridge.Instance.BButtonDown){print("kkkkkkbbbbb456");} 2抓取某个物体&#xff0c;那么就在要抓取的那个物体上加一些组件&#xff0c;特别是Grabble Unity Events

cocos2dx DrawNode

cocos2dx 两种绘图方式 DrawPrimitivesDrawNode DrawPrimitives 3.x 已经弃用 绘制的图形可以是实心的&#xff0c;也可以是空心的。 DrawNode 在一个单独的批处理中绘制了所以元素&#xff0c;因此它绘制点、线段、多边形都要比“drawing primitives”快。 绘制的图形都…

【基础知识】AB软件RSLinx如何实现OPC通讯组态

哈喽&#xff0c;大家好&#xff0c;我是雷工。 在上一节了解了什么是RSLinx&#xff1f;以及RSLinx Lite、RSLinx Classice、RSLinx Professional、RSLinx Gateway几个版本的特点。 本节了解AB的RSLinx如何实现OPC组态。 一、创建RSLinx通讯&#xff1a; 1.1、【Communicati…

excel自己记录

1、清除换行符号 2、添加特殊符号&并清除换行符号 7日&15日&30日&60日 3、判断单元格最后一个字符是不是数字&#xff0c;不是就删掉 IF(ISNUMBER(--RIGHT(B2,1)),B2,SUBSTITUTE(B2,RIGHT(B2,1),"")) ISNUMBER(--RIGHT(B2,1))判断最右边的一个数是否…

Wireshark的捕获过滤器

Wireshark的过滤器&#xff0c;顾名思义&#xff0c;作用是对数据包进行过滤处理。具体过滤器包括捕获过滤器和显示过滤器。本文对捕获过滤器进行分析。 捕获过滤器&#xff1a;当进行数据包捕获时&#xff0c;只有那些满足给定的包含/排除表达式的数据包会被捕获。 捕获过滤器…

一起学docker系列之九docker运行mysql 碰到的各种坑及解决方法

目录 前言1 Docker 运行mysql命令2 坑一&#xff1a;无法读取/etc/mysql/conf.d目录的问题3 坑二&#xff1a;/tmp/ibnr0mis 文件无法创建/写入的问题4 坑三&#xff1a;Navicat 连接错误&#xff08;1045-access denied&#xff09;5 坑四&#xff1a;MySQL 登录失败问题结语 …

ros2文件package.xml与cmakelists.txt比较

每次在ros2里面添加文件以后&#xff0c;都要修改packages.xml,与cmakelists.txt文件。