1

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
2
  • 1
    What happens with dst.write(array) or dst.write(array, [1,2,3]) (or 0,1,2?) instead? Aug 3, 2023 at 18:43
  • 1
    All bands got the value zero, or I got the Source shape error. Thanks anyway
    – BERA
    Aug 4, 2023 at 5:30

3 Answers 3

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)
0

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)
0
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

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.