Rasterio - Extract raster pixel values based on another raster's pixel values

Question

I have two single band rasters, both with the same spatial resolution and CRS. A land cover raster (named clc), and another one, that is for a smaller area (named patch) and falls within the extent of clc. Here are some details of the two rasters:

	clc	patch
Number of Bands	1	1
CRS	EPSG:32632	EPSG:32632
Spatial resolution	(60.0, 60.0)	(60.0, 60.0)
Rows	14452	135
Columns	12067	134

I want to extract a new raster with the extent and pixel values of patch on the locations where clc has values defined in a list [10, 13, 22].

I am new to rasterio and only managed to extract the clc values instead of the corresponding patch values.

Answer 1

You can use numpy.isin to create a boolean array where clc values are not in [10, 13, 22] and use that to mask patch with nodata. This means you will retain only patch values where clc values are in [10, 13, 22].

As your clc raster is larger, you'll need to read only from a window matching the bounds (extent) of the patch raster:

E.g.

import numpy as np
import rasterio as rio
from rasterio.windows import from_bounds

with rio.open('patch.tif') as patch:
    patchdata = patch.read()
    profile = patch.profile.copy()
    bounds = patch.bounds

with rio.open('clc.tif') as clc:
    window = from_bounds(*bounds, clc.transform)
    clcdata = clc.read(window=window)

# create bool array where clc is not in list of values
mask = np.isin(clcdata, [10, 13, 22], invert=True)

# replace None with lowest value in dtype range
if profile["nodata"] is None:
    profile["nodata"] = np.iinfo(patchdata.dtype).min


# set everything where mask is True to NoData
patchdata[mask] = profile["nodata"]

with rio.open('output.tif', "w", **profile) as output:
    output.write(patchdata)