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.

Improve this question

1 Answer 1

Reset to default
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)
Improve this answer
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? 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. 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? 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. Commented Feb 6, 2023 at 20:08

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.