1

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.

2 Answers 2

3

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
4
  • 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?
    – GCGM
    Oct 3, 2019 at 10:29
  • 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.
    – user2856
    Oct 3, 2019 at 10:38
  • Coul you please ellaborate on how to mask the array with the shapefile with some example?
    – GCGM
    Oct 3, 2019 at 10:50
  • 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.
    – user2856
    Oct 3, 2019 at 10:50
1

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)

Your Answer

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

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