Python - Mask raster array using shapefile or rasterio

Question

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.

Answer 1

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

• GCGM: Thanks. So, no option to concatenate a resample + mask operatiosn using rasterio. If I want to avoid the writting the resampled product, is it possible to mask the array using the shapefile? I am seeing some posts about using rasterio.warp Would that be an option?
• user2856: Yes, you just need to rasterize the shapefile with the same dimensions & transform as the resampled array. I don't know about warp, don't use it.
• GCGM: Coul you please ellaborate on how to mask the array with the shapefile with some example?
• user2856: Sorry, don't have one handy. See the doc rasterio.readthedocs.io/en/stable/api/…. You still need to adjust the transform per my answer.

Answer 2

You might be interested in rioxarray.

import rioxarray
xds = rioxarray.open_rasterio(...)

https://corteva.github.io/rioxarray/stable/examples/reproject.html

xds_lonlat = xds.rio.reproject("epsg:4326")

https://corteva.github.io/rioxarray/stable/examples/clip_geom.html

clipped = xds.rio.clip(features, features_crs)