0

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.

1 Answer 1

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)
4
  • Thaanks for your response. I receive an error "IndexError: boolean index did not match indexed array along dimension 0; dimension is 135 but corresponding boolean dimension is 14452" Does that mean that I have to limit the operation to the extent of patch?
    – KosV
    Commented Feb 6, 2023 at 12:06
  • Yes, both arrays have to be the same shape. You neglected to mention that they differed in your question. Please edit it to specify the no of rows, cols and bands in each raster.
    – user2856
    Commented Feb 6, 2023 at 12:23
  • Thank you, I've edited my original question. Do you happen to know a method to do this?
    – KosV
    Commented Feb 6, 2023 at 13:10
  • That's great, I did something similar in the end but your way has fewer lines. I also added a condition to replace None values with 0, because otherwise it throws an error. Thank you so much for your time.
    – KosV
    Commented Feb 6, 2023 at 20:08

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.