Write three band RGB to file

Question

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

Answer 1

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)

Answer 2

You should take a look a the docs for the write function.
The second argument:

indexes (int or list, optional) – Which bands of the dataset to write to. The default is all.

Here when you pass the value 3, you're asking rasterio to write only the third band.
dst.write(array) should work for a shape == (3, 10980, 10980).
For 2d arrays (single band) the right way to do it is with dst.write_band(1, b4).

Edit: just saw the comment above. If write(array) produces only zeros, you may have a problem elsewhere. Possibly because of the reading part, and the JP2 driver. You can try this:

with rasterio.open(bluefile) as ds2:
    b2 = ds2.read(1)
    profile = ds2.profile.copy()
with rasterio.open(greenfile) as ds3:
    b3 = ds3.read(1)
with rasterio.open(redfile) as ds4:
    b4 = ds4.read(1)

array = np.stack((b4,b3,b2))
print(array.shape)

print(profile)
profile["count"] = 3
profile["driver"] = "GTiff"

with rasterio.open(outfile, "w", **profile) as dst:
    dst.write(array)

Answer 3

from osgeo import gdal, gdal_array
import numpy as np

# input

band2 = gdal.Open(r'My inputpath')
band3 = gdal.Open(r'My inputpath')
band4 = gdal.Open(r'My inputpath')

array2 = band2.ReadAsArray()
array3 = band3.ReadAsArray()
array4 = band4.ReadAsArray()


stack = np.array([array2,array3,array4])

gdal_array.SaveArray(stack, r'My outputpath", prototype=band2)
                     
banda2 = None