如何用Rasterio将JP2格式R、G、B波段合成为RGB影像？

我需要将三个对应R、G、B波段的.jp2文件（Sentinel-2的B04=红、B03=绿、B02=蓝）合并为单个RGB文件，使用Rasterio库处理，编写的代码如下：

import os
import geopandas
import rasterio as rio
from PIL import Image

Image.MAX_IMAGE_PIXELS = 300000000

## RGB Transform

# Image paths: B04 = Red, B03 = Green, B02 = Blue
redPath = "assets/geodata/dump_r2/S2A_MSIL2A_20220512T133231_N0400_R081_T22JCP_20220512T202012.SAFE/GRANULE/L2A_T22JCP_A035971_20220512T133948/IMG_DATA/R10m/T22JCP_20220512T133231_B04_10m.jp2"
greenPath = "assets/geodata/dump_r2/S2A_MSIL2A_20220512T133231_N0400_R081_T22JCP_20220512T202012.SAFE/GRANULE/L2A_T22JCP_A035971_20220512T133948/IMG_DATA/R10m/T22JCP_20220512T133231_B03_10m.jp2"
bluePath = "assets/geodata/dump_r2/S2A_MSIL2A_20220512T133231_N0400_R081_T22JCP_20220512T202012.SAFE/GRANULE/L2A_T22JCP_A035971_20220512T133948/IMG_DATA/R10m/T22JCP_20220512T133231_B02_10m.jp2"

## Bands

"""
redConv = Image.open(redPath)
redConv.save("red.tiff", "TIFF")
greenConv = Image.open(greenPath)
greenConv.save("green.tiff", "TIFF")
blueConv = Image.open(bluePath)
blueConv.save("blue.tiff", "TIFF")
"""
red = rio.open(redPath)
green = rio.open(greenPath)
blue = rio.open(bluePath) 

# Creates RGB file.
rgb = rio.open("RGB.tiff", 'w+', driver="Gtiff", width=red.width, height=red.height, count=3, crs=red.crs, transform=red.transform, dtype=red.dtypes[0])
rgb.write(blue.read(1), 1)
rgb.write(green.read(1), 2)
rgb.write(red.read(1), 3)

print(red.shape)
print(green.shape)
print(blue.shape)

rgb.close()

但生成的RGB.tiff在GIMP中显示为灰度图像，无彩色通道。我的思考：

• 驱动可能存在问题，涉及JP2转TIFF、JP2读取或TIFF保存环节；
• 可能缺少影像处理步骤，目前未找到解决方向；
• rgb.read().shape为[3,10980,10980]，与JPEG的(1980,1080,3)形状不同，疑问TIFF的可视化特性。

GIMP打开的图像：

影像元数据（exiftool输出）：

ExifTool Version Number         : 12.16
File Name                       : RGB.tiff
Directory                       : .
File Size                       : 690 MiB
File Modification Date/Time     : 2022:08:24 11:56:03-03:00
File Access Date/Time           : 2022:08:24 11:55:18-03:00
File Inode Change Date/Time     : 2022:08:24 11:56:03-03:00
File Permissions                : rw-r--r--
File Type                       : TIFF
File Type Extension             : tif
MIME Type                       : image/tiff
Exif Byte Order                 : Little-endian (Intel, II)
Image Width                     : 10980
Image Height                    : 10980
Bits Per Sample                 : 16 16 16
Compression                     : Uncompressed
Photometric Interpretation      : BlackIsZero
Strip Offsets                   : (Binary data 108116 bytes, use -b option to extract)
Samples Per Pixel               : 3
Rows Per Strip                  : 1
Strip Byte Counts               : (Binary data 65879 bytes, use -b option to extract)
Planar Configuration            : Chunky
Extra Samples                   : Unknown (0 0)
Sample Format                   : Unsigned; Unsigned; Unsigned
Pixel Scale                     : 10 10 0
Model Tie Point                 : 0 0 0 300000 6900040 0
Geo Tiff Version                : 1.1.0
GT Model Type                   : Projected
GT Raster Type                  : Pixel Is Area
GT Citation                     : WGS 84 / UTM zone 22S
Geog Citation                   : WGS 84
Geog Angular Units              : Angular Degree
Projected CS Type               : WGS84 UTM zone 22S
Proj Linear Units               : Linear Meter
Image Size                      : 10980x10980
Megapixels                      : 120.6

问题根源 #

1. Photometric标识错误：从元数据可见Photometric Interpretation为BlackIsZero，这是单波段灰度图的默认标识，GIMP会以此判定图像为灰度图，而非RGB彩色图。
2. 16位数据未适配：Sentinel-2的L2A数据是16位无符号整数，多数图像查看工具（包括GIMP）默认不会自动将16位多波段图像渲染为彩色，需要转换为8位。
3. 波段顺序颠倒：原代码将蓝波段写入第1通道、红波段写入第3通道，输出的是BGR顺序，会导致颜色失真。

修改后的代码 #

import rasterio as rio
import numpy as np

# 波段路径
red_path = "assets/geodata/dump_r2/S2A_MSIL2A_20220512T133231_N0400_R081_T22JCP_20220512T202012.SAFE/GRANULE/L2A_T22JCP_A035971_20220512T133948/IMG_DATA/R10m/T22JCP_20220512T133231_B04_10m.jp2"
green_path = "assets/geodata/dump_r2/S2A_MSIL2A_20220512T133231_N0400_R081_T22JCP_20220512T202012.SAFE/GRANULE/L2A_T22JCP_A035971_20220512T133948/IMG_DATA/R10m/T22JCP_20220512T133231_B03_10m.jp2"
blue_path = "assets/geodata/dump_r2/S2A_MSIL2A_20220512T133231_N0400_R081_T22JCP_20220512T202012.SAFE/GRANULE/L2A_T22JCP_A035971_20220512T133948/IMG_DATA/R10m/T22JCP_20220512T133231_B02_10m.jp2"

def scale_to_8bit(arr):
    """将16位遥感数据缩放到8位范围"""
    min_val = arr.min()
    max_val = arr.max()
    # 避免除以0的情况
    if max_val == min_val:
        return np.zeros_like(arr, dtype=np.uint8)
    return ((arr - min_val) / (max_val - min_val) * 255).astype(np.uint8)

# 使用上下文管理器确保文件正确关闭
with rio.open(red_path) as red_src, rio.open(green_path) as green_src, rio.open(blue_path) as blue_src:
    # 读取波段数据
    red_data = red_src.read(1)
    green_data = green_src.read(1)
    blue_data = blue_src.read(1)
    
    # 转换为8位
    red_8bit = scale_to_8bit(red_data)
    green_8bit = scale_to_8bit(green_data)
    blue_8bit = scale_to_8bit(blue_data)
    
    # 创建输出文件，关键设置photometric='RGB'
    with rio.open(
        "RGB_color.tiff", "w",
        driver="GTiff",
        width=red_src.width,
        height=red_src.height,
        count=3,
        crs=red_src.crs,
        transform=red_src.transform,
        dtype=np.uint8,
        photometric="RGB"
    ) as dst:
        # 按RGB顺序写入通道1、2、3
        dst.write(red_8bit, 1)
        dst.write(green_8bit, 2)
        dst.write(blue_8bit, 3)

print("彩色RGB TIFF已生成")

关键修改说明 #

• 设置Photometric参数：显式指定photometric="RGB"，让图像工具识别这是彩色图像。
• 16位转8位：通过线性缩放将Sentinel-2的0-65535范围数据转换为0-255的8位数据，适配普通图像查看工具。
• 修正波段顺序：按红、绿、蓝的顺序写入通道1、2、3，确保颜色显示正常。
• 使用上下文管理器：替代手动close()，避免文件资源泄漏。

内容的提问来源于stack exchange，提问作者Oz_Magician