Merging Sentinel 2 RGB bands with rasterio

Question

I'm trying to merge Sentinel 2 bands RGB bands with rasterio like this:

import rasterio
import numpy as np
from rasterio import plot
from rasterio.plot import show

band2=rasterio.open("B02.jp2")
band3=rasterio.open("B03.jp2")
band4=rasterio.open("B04.jp2")

rgb=rasterio.open('rgb.tiff', 'w', driver='Gtiff',
                          width=band2.width, height=band2.height,
                          count=3,
                          crs=band2.crs,
                          transform=band2.transform,
                          dtype='uint16')
rgb.write(band4.read(1),1)
rgb.write(band3.read(1),2)
rgb.write(band2.read(1),3)
rgb.close()

This works as expected and creates a RGB-GeoTIFF that is properly displayed in QGIS. However if I try to visualize it with rasterio, I get an almost black image:

show(plot.adjust_band(rasterio.open('rgb.tiff').read([1,2,3])))

when I inspect the normalized values the mean for each band is around 0.05. Is there something wrong how I create the RGB?

Answer 1

I experienced the same "issue" when working on Sentinel-2 images. Your code looks fine but incomplete. QGIS automatically rescales image intensity, which explains why your images is displayed correctly. However, rasterio does not proceed likewise, so you need to rescale image intensity before plotting. I suggest you to use skimage for achieving it. See the following code (worked well for me):

import rasterio
import numpy as np
from rasterio import plot
from rasterio.plot import show
from skimage import exposure

band2=rasterio.open("B02.jp2")
band3=rasterio.open("B03.jp2")
band4=rasterio.open("B04.jp2")

band2_geo = band2.profile
band2_geo.update({"count": 3})

with rasterio.open('rgb.tiff', 'w', **band2_geo) as dest:
# I rearanged the band order writting to 2→3→4 instead of 4→3→2
    dest.write(band2.read(1),1)
    dest.write(band3.read(1),2)
    dest.write(band4.read(1),3)

# Rescale the image (divide by 10000 to convert to [0:1] reflectance
img = rasterio.open('rgb.tiff')
image = np.array([img.read(3), img.read(2), img.read(1)]).transpose(1,2,0)
p2, p98 = np.percentile(image, (2,98))
image = exposure.rescale_intensity(image, in_range=(p2, p98)) / 100000

# Plot the RGB image
show(image.transpose(2,0,1), transform=img.transform)