I have a raster which I want to crop multiple times using different shapefiles. My code uses mask method which returns numpy array and its metadata:

out_img, out_transform = mask(raster=data, shapes=coords, crop=True)

I run different kinds of thresholding on out_img and I want to crop it again using another shapefile but I'm not sure how to do that since out_img is simply a numpy array.

Do I need to save out_img and out_transform as TIFF, load it and use mask method again? Is there any more straightforward method of achieving it (i.e. without saving the intermediate steps)?


You don't necessarily need to write your crop to disk, you can write it to a rasterio in-memory file.

Your example would look something like this:

with rasterio.open("your_file.tif") as src:
    out_img, out_transform = mask(raster=src, shapes=coords, crop=True)

    # Update the `src`'s profile with the crop's attributes
    profile = src.profile
    profile["height"] = out_image.shape[1]
    profile["width"] = out_image.shape[2]
    profile["transform"] = out_transform

# Write the crop to a memfile
with rasterio.MemoryFile() as memfile:
    with memfile.open(**profile) as dst:

        # Crop the memfile with the new shape `coords2`
        out_image2, out_transform2 = mask(dst, shapes=coords2, crop=True)

However just a note: if the new shapefile you are using to crop your out_image is independent of the thresholding step you describe, you could simply crop with the list of shapefiles straight away, since mask() takes an iterator of shapes.

| improve this answer | |

When using rasterio to work on rasters in memory as numpy arrays, and doing multiple masking steps together with similar operations, I usually rasterize the masks to numpy arrays. Then it is possible to do everything without using rasterio.MemoryFile. See rasterio.features.rasterize.

So say I've done your first step, then want to load and apply another mask:

out_img, out_transform = mask(raster=data, shapes=coords, crop=True)
new_img = out_img * 2  # Some intermediate operation
new_mask = rasterio.features.rasterize(new_mask_geometries, out_shape=new_img.shape,
                                       transform=out_transform, default_value=1)
new_img[new_mask==1] = np.nan

Note: this doesn't do any cropping, just masking of values.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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