一、前言
下面的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 = None
out_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.13
beta=0
return 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=2
return 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')
