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?
Best Answer
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):