4

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?

2
  • please check the nodata value of your bands. maybe 65535 it is considered as a value so that the adjustment is not done on the correct range. rasterio plot options are unfortunately quite limited, while QGIS does it better.This is not a merging but a display isse as far as I can tell from your explanations.
    – radouxju
    Commented Nov 14, 2019 at 10:40
  • The pixel values before the normalization are between 1 and around 28000. The relevant output of gdalinfo doesn't yield infos about no_data: Band 1 Block=10980x1 Type=UInt16, ColorInterp=Gray Band 2 Block=10980x1 Type=UInt16, ColorInterp=Undefined Band 3 Block=10980x1 Type=UInt16, ColorInterp=Undefined
    – user153583
    Commented Nov 14, 2019 at 10:52

1 Answer 1

7

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)
2
  • Awesome, I've been looking for this. Do you know this issue exists? C
    – Geosphere
    Commented Jan 3, 2022 at 16:30
  • Can the data range for the bands be transformed in values between 0-255?
    – Geosphere
    Commented Jan 3, 2022 at 16:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.