Summing four rasters to one suitability raster using Rasterio & GeoPandas in Python

Question

I am trying to create a basic suitability raster from three vector sources and one raster source using Rasterio & GeoPandas in Python. I have already clipped and rasterized the vector sources to the raster extent.

All rasters contain nodata values. How do I combine all four rasters into one raster? I want to:

• Keep nodata values only if all 4 rasters have nodata on the cell.
• Replace nodata if 1 or more rasters contain data on the cell.

My current implementation does almost what I want, but replaces all nodata by 0 and will not write the result to a new raster.

import rasterio
import rasterio.plot
import numpy as np

# Merge the rasterized vectors

# Read raster bands directly to Numpy arrays.

with rasterio.open('water.tif') as src:
    water = src.read()
    print(water.dtype)

with rasterio.open('road.tif') as src:
    roads = src.read()

with rasterio.open('natura2000.tif') as src:
    natura2000 = src.read()

with rasterio.open('elevation.tif') as src:
    elevation = src.read()

suit_raster = np.zeros(elevation.shape, dtype=rasterio.int16)
for data in water, roads, natura2000, elevation:
    suit_raster += data

# Update metadata
kwargs = src.meta
kwargs.update(
    dtype=rasterio.uint16,
    count=1,
    compress='lzw',
    nodata=0
)

# Show resulting raster.
rasterio.plot.show(suit_raster)

# Write result to file fails?
with rasterio.open('a5_suit_raster.tif', 'w', **kwargs) as dst:
    dst.write_band(1, suit_raster.astype(rasterio.uint16))

This results in the following raster which is almost what I want as the nodata values are all replaced by 0 (or is this because it is not a raster yet and I am plotting an numpy array?):

Additionaly, I cant write the result to a new raster as I get the error message:

Source shape (1, 1, 1437, 1628) is inconsistent with given indexes 1

What am I missing here? How do I update the 'shape' in the new rasters? And how do I keep the nodata values?

Answer 1

After re-reading the rasterio docs I changed the writing to file step to make it work:

• Use write() rather than write_band() in writing to a raster file.

The resulting 'raster' shown above was a plot of a plain numpy array, therefore nodata values were not shown. The new resulting raster now looks like this as was expected:

import rasterio
import rasterio.plot
import numpy as np

# Merge the rasterized vectors

# Read raster bands directly to Numpy arrays.

with rasterio.open('water.tif') as src:
    water = src.read()
    print(water.dtype)

with rasterio.open('road.tif') as src:
    roads = src.read()

with rasterio.open('natura2000.tif') as src:
    natura2000 = src.read()

with rasterio.open('elevation.tif') as src:
    elevation = src.read()

suit_raster = np.zeros(elevation.shape, dtype=rasterio.int16)
for data in water, roads, natura2000, elevation:
    suit_raster += data

# Update metadata
kwargs = src.meta
kwargs.update(
    dtype=rasterio.uint16,
    count=1,
    compress='lzw',
    nodata=0
)

with rasterio.open('data/a5_suit_raster.tif', 'w+', **kwargs) as dst:
    # Don't use write_band, but use write!
    dst.write(suit_raster.astype(rasterio.uint16))