[GIS] Writing a 3-band image with rasterio

geotiff-tiffpythonrasterio

I would like to save a 3-band raster to file using rasterio.

The official docs have examples for saving a single band image in both the quickstart and in the advanced topics.

What I would like to do is something like..

with rasterio.open('image.tif', 'w', **meta) as dst:
    dst.write(array)

This raises the following error for an array with a shape of (256, 256, 3).

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

Best Answer

The numpy array needs to be in bands, rows, cols order (z, y, x) not cols, rows, bands (x, y, z).

You can rearrange the axes with numpy.moveaxis

e.g.

import numpy as np
import rasterio

with rasterio.open('image.tif', 'w', **meta) as dst:
    # If array is in (x, y, z) order (cols, rows, bands)
    dst.write(np.moveaxis(array, [0, 1, 2], [2, 1, 0]))

    # If array is in (y, x, z) order (rows, cols, bands)
    dst.write(np.moveaxis(array, [0, 1, 2], [2, 0, 1]))  # same as np.rollaxis(array, axis=2)

Related Solutions

[GIS] Windowed writing in rasterio

You can directly write the result blocks to the big result tiff file via windowed writing.

import rasterio

src_RP1 = rasterio.open(r'RP1_filename.tif')
src_RP2 = rasterio.open(r'RP2_filename.tif')

out_profile = src_RP1.profile.copy()
dst = rasterio.open(r'result.tif', 'w', **out_profile)

for block_index, window in src_RP1.block_windows(1):
    RP1_block = src_RP1.read(window=window, masked=True)
    RP2_block = src_RP2.read(window=window, masked=True)

    result_block = RP2_block - RP1_block
    dst.write(result_block, window=window)

src_RP1.close()
src_RP2.close()
dst.close()

If indeed you want to assemble one big results array in memory you would have to initialize an empty array with the dimensions of RP1 and iteratively assign the block arrays to the respective rows/columns (you can get these from the window).

[GIS] Reading a 4 band image with rasterio (python)

Reading and plotting works just fine, but the problem is in the warning you get from matplotlib:

If an image has integer values matplotlib will clip them to interval [0, 255]
If an image has float values matplotlib will clip them to interval [0, 1.0]

In your case all values were clipped to either 255 or to 1.0 and that is why you see a white image. Also note that this happens only when you try to plot a RGB image. When plotting an image with a single channel matplotlib won't do any clipping.

Solution - check dtype of your image values and decrease / change them accordingly. Example:

image = dataset.read([1,2,3])
print(image.dtype, np.max(image))
image = (255 * image / np.max(image)).astype(np.uint8)
plot.show(image)

This transforms image values to integer values from [0, 255]. (Instead by multiplying images by 255 / np.max(image) you might want to try some other constant to get nicer image. I used this one because I don't know anything about values of your image.)

Best Answer

Related Solutions

[GIS] Windowed writing in rasterio

[GIS] Reading a 4 band image with rasterio (python)

Related Question