Roger Almengor
Roger Almengor

Reputation: 464

Cannot create a virtual raster from a stack array in rasterio

I Need to know how to create a virtual raster. I iterated over a folder, which contains some binary rasters (1, 0). Those rasters I append it into a numpy array using numpy.concatenate. Then I would like to create a virtual raster using the number of raster concatenated as the number of bands that this raster will have. I receive the followgin message though:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-93-7a57711dc737> in <module>
     20     compress = 'lzw')
     21     with rasterio.open(path_scl + "/" + "scl_stack.vrt", "w", **profile) as dst:
---> 22         dst.write(final_array.astype(rasterio.uint8), dimension)

rasterio/_io.pyx in rasterio._io.DatasetWriterBase.write()

IndexError: band index out of range

I checked the number of rasters, and it corresponds to the variable "dimension", y just put when writing out my final virtual raster.

path_scl = r'I:\Sentinel-2\Central\2017\T32TNT'
files = [os.path.join(root, file) for root, directories, filenames in    os.walk(path_scl) for file in filenames]
scls = [file for file in files if    file.endswith("_01_cloud_mask_bin.tif")]
final_array = np.zeros((10980, 10980))
for scl in scls:
with rasterio.open(scl) as ds:

    profile = ds.profile 
    array = ds.read(1)
    np.concatenate((final_array, array), axis = 0)
    print(f"{scl} added")

dimension = len(scls)

with rasterio.Env():

    profile = profile 

    profile.update(
    dtype = rasterio.uint8,
    compress = 'lzw')
    with rasterio.open(path_scl + "/" + "scl_stack.vrt", "w", **profile) as dst:
        dst.write(final_array.astype(rasterio.uint8), dimension)

Does anyone know how to interprete this error message? Thanks

Upvotes: 0

Views: 879

Answers (1)

Christoph Rieke
Christoph Rieke

Reputation: 757

VRT is a read-only format, you can not write arrays to VRT. To create a VRT file you only reference the existing, on-disk raster files. E.g. using gdal.BuildVRT in Python with your list of raster files as in example 3 here.

Upvotes: 1

Related Questions