3

I have a satellite image raster file and a Shapefile consisting of several discrete polygons. I would like to individually extract (or mask) the Numpy image array for each of these polygons from the raster file.

Currently, I import the Shapefile using Fiona and create a list of the polygons. I have no problem masking the raster file using all the polygons at once. However when I try to use an individual polygon from the list I get an error: "ValueError: Input shapes do not overlap raster.". This is despite getting a successful result previously, and both files using the same CRS.

import rasterio
import shapefile
import fiona
import numpy as np

with fiona.open("test.shp", "r") as shapefile:

    features = [feature["geometry"] for feature in shapefile]

features = [x for x in features if x is not None]

This works:

with rasterio.open('sat_img_B01.jp2') as src:

    out_image, out_transform = rasterio.mask.mask(src, features,
                                                        crop=True)

This doesn't work (WindowError: windows do not intersect and ValueError: Input shapes do not overlap raster):

index = 0

with rasterio.open('sat_img_B01.jp2') as src:

    out_image, out_transform = rasterio.mask.mask(src, [features[index]], crop=True)

I guess I'm missing something fundamental! Is there an elegant way of extracting each polygon in the list 'features' as an individual image file/numpy array from the raster image?

Thanks!!

Will000
  • 31
  • 1
  • 4

1 Answers1

4

You have to make sure that both the shapefile and the raster image are using the same crs (Coordinate Reference System): try printing them.

Firas Omrane
  • 894
  • 1
  • 14
  • 21
  • can you elaborate this point? I'm experiencing the same problem – Claus Jan 07 '20 at 22:19
  • 2
    @Claus I remember that my problem was that the shapefile and the raster image are not loaded with the same crs (coordinate reference system). I think that you already loaded both of them (shapefile using fiona) and (raster image using resterio). For the shapefile to print the crs using `.crs` property and for the resterio try to load the raster image using `src = rasterio.open('image')` and them print its crs `src.crs`. If the two crs are not the same the mask won't work. The solution is to change the crs of one of them to the other. I know this is late. – Firas Omrane May 30 '20 at 23:57