Python Rasterio – Write Three Band RGB to File

pythonrasterio

I'm trying to create a 3 band RGB image from Sentinel 2 bands but I seem to be trying to write an array of inconsistent shape:

For example:

ValueError: Source shape (1, 10980, 10980, 3) is inconsistent with
given indexes 1

I've tried reshaping it but I still get errors.

import numpy as np
import rasterio, os

folder = r"/home/bera/data/Sentinel2/S2B_MSIL2A_20200625T101559_N0214_R065_T33VWD_20200625T134347.SAFE/GRANULE/L2A_T33VWD_A017251_20200625T101903/IMG_DATA/R10m/"
bluefile = os.path.join(folder, "T33VWD_20200625T101559_B02_10m.jp2")
greenfile = os.path.join(folder, "T33VWD_20200625T101559_B03_10m.jp2")
redfile = os.path.join(folder, "T33VWD_20200625T101559_B04_10m.jp2")

#Read as datasets
dsblue = rasterio.open(bluefile)
dsgreen = rasterio.open(greenfile)
dsred = rasterio.open(redfile)

#And to arrays
blue = dsblue.read(1)
green = dsgreen.read(1)
red = dsred.read(1)

#Composite 
rgb = np.dstack((red, green, blue))
print(rgb.shape)
#(10980, 10980, 3)

def writetofile(output, array):
    with rasterio.open(fp=outfile, mode="w", driver="GTiff", 
                       width=rgb.shape[1], height=rgb.shape[0],
                       count=3, crs=dsblue.crs, 
                       transform=dsblue.transform, dtype=rgb.dtype) as dst:
        dst.write(array, 3)

outfile = r"/home/bera/Desktop/GIStest/rgbtest.tif"
writetofile(output=outfile, array=rgb)
#ValueError: Source shape (1, 10980, 10980, 3) is inconsistent with given indexes 1

from rasterio.plot import reshape_as_raster, reshape_as_image

raster = reshape_as_raster(rgb)
print(raster.shape)
#(3, 10980, 10980)
writetofile(output=outfile, array=raster)
#ValueError: Source shape (1, 3, 10980, 10980) is inconsistent with given indexes 1

image = reshape_as_image(rgb)
print(image.shape)
#(10980, 3, 10980)
writetofile(output=outfile, array=image)
#ValueError: Source shape (1, 10980, 3, 10980) is inconsistent with given indexes 1

Best Answer

I found that it is possible to write the bands one by one to the same output file, you dont need to stack them.

This is working:

import rasterio, os

folder = r"C:\S2_230731\S2A_MSIL2A_20230611T103631_N0509_R008_T33VUH_20230611T173458\S2A_MSIL2A_20230611T103631_N0509_R008_T33VUH_20230611T173458.SAFE\GRANULE\L2A_T33VUH_A041618_20230611T103626\IMG_DATA\R10m"
bluefile = os.path.join(folder, "T33VUH_20230611T103631_B02_10m.jp2")
greenfile = os.path.join(folder, "T33VUH_20230611T103631_B03_10m.jp2")
redfile = os.path.join(folder, "T33VUH_20230611T103631_B04_10m.jp2")

outfile = r"C:\folder\rgbtest.tif"

ds2=rasterio.open(bluefile)
ds3=rasterio.open(greenfile)
ds4=rasterio.open(redfile)
b2 = ds2.read(1)
b3 = ds3.read(1)
b4 = ds4.read(1)

with rasterio.open(fp=outfile, mode="w", driver="GTiff", width=ds2.shape[1], height=ds2.shape[0], 
                   count=3, crs=ds2.crs, transform=ds2.transform, dtype=b2.dtype) as dst:
    dst.write(b2, 3)
    dst.write(b3, 2)
    dst.write(b4, 1)
Related Question