基于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')