1

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?

  • 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 Nov 14 '19 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 Nov 14 '19 at 10:52
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)
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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