4

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?):

plot of the input rasters and the resulting output 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?

1 Answer 1

2

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: resulting suitability 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
)

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))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.