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?