[GIS] Python – Mask raster array using shapefile or rasterio

arraypythonrasterrasterioshapefile

I have a raster that I have resample using rasterio, where img.shape is (5490, 5490)

if img.res[1] == 20:
      img = img.read(
      out_shape=(img.count, img.height * 2, img.width * 2),
      resampling=Resampling.nearest)

The Resampling function of rasterio is giving me the result as an array

img 
array([[[2014, 2014, 1999, ..., 1705, 2136, 2136],
        [2014, 2014, 1999, ..., 1705, 2136, 2136],
        [1997, 1997, 1947, ..., 1709, 1769, 1769],
        ...,
        [1871, 1871, 1975, ..., 1868, 1966, 1966],
        [1817, 1817, 1892, ..., 1824, 1870, 1870],
        [1817, 1817, 1892, ..., 1824, 1870, 1870]]], dtype=uint16)

However, in my next step I want to use rasterio again to mask my raster with a shapefile using the following code:

out_image, out_transform = rasterio.mask.mask(img, features, crop=True)

This is not working as the input for mask (img) is not a rasterio object but an array. 

AttributeError: 'numpy.ndarray' object has no attribute 'nodata'

What should be the procedure in this case. From other posts I have seen that an option would be to rasterize my shapefile and read it as an array but I would like avoid work-around solutions.

Best Answer

It's a bit of a pain, but you need to write the resampled numpy array to a rasterio Dataset (either on file or in memory) and adjust the transform to match the resampling..

Here's an example of both:

# Example licensed under cc by-sa 4.0 with attribution required

from contextlib import contextmanager

import rasterio
from rasterio import Affine, MemoryFile
from rasterio.enums import Resampling


@contextmanager
def resample_raster(raster, out_path=None, scale=2):
    """ Resample a raster
        multiply the pixel size by the scale factor
        divide the dimensions by the scale factor
        i.e
        given a pixel size of 250m, dimensions of (1024, 1024) and a scale of 2,
        the resampled raster would have an output pixel size of 500m and dimensions of (512, 512)
        given a pixel size of 250m, dimensions of (1024, 1024) and a scale of 0.5,
        the resampled raster would have an output pixel size of 125m and dimensions of (2048, 2048)
        returns a DatasetReader instance from either a filesystem raster or MemoryFile (if out_path is None)
    """
    t = raster.transform

    # rescale the metadata
    transform = Affine(t.a * scale, t.b, t.c, t.d, t.e * scale, t.f)
    height = int(raster.height / scale)
    width = int(raster.width / scale)

    profile = src.profile
    profile.update(transform=transform, driver='GTiff', height=height, width=width)

    data = raster.read(
            out_shape=(raster.count, height, width),
            resampling=Resampling.bilinear,
        )

    if out_path is None:
        with write_mem_raster(data, **profile) as dataset:
            del data
            yield dataset

    else:
        with write_raster(out_path, data, **profile) as dataset:
            del data
            yield dataset


@contextmanager
def write_mem_raster(data, **profile):
    with MemoryFile() as memfile:
        with memfile.open(**profile) as dataset:  # Open as DatasetWriter
            dataset.write(data)

        with memfile.open() as dataset:  # Reopen as DatasetReader
            yield dataset  # Note yield not return


@contextmanager
def write_raster(path, data, **profile):

    with rasterio.open(path, 'w', **profile) as dataset:  # Open as DatasetWriter
        dataset.write(data)

    with rasterio.open(path) as dataset:  # Reopen as DatasetReader
        yield dataset
Related Question