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?