Masking GeoTIFF – Solve ‘Input Shapes Do Not Overlap Raster’ with GeoJSON and Rasterio

geojsongeotiff-tiffmaskingplanetrasterio

So my intention is quite simple:
I have a 4-Band TIFF file I want to crop after a polygon shape. However, I'm getting the following two exceptions

  • WindowError: windows do not intersect
  • ValueError: Input shapes do not overlap raster.

after running this code

with rasterio.open(raw_filename) as dataset:
   out_image, out_transform = rasterio.mask.mask(dataset, [polygon], crop=True)

Now this seems really strange to me because the same polygon was used to query my data. Moreover, I made sure there's a 100% overlap between them
prior to downloading and I can see it that the areas overlap.

My polygon is the blue box below:

AOI

with the values:

print(polygon)
# {'type': 'Polygon', 'coordinates': [[[36.044253, 5.166209], [36.044253, 5.180148], [36.061434, 5.180148], [36.061434, 5.166209], [36.044253, 5.166209]]]}

And my file looks something like this:

image = rasterio.open(raw_filename)
rasterio.plot.show((image, 1))

Full shape

which might not mean much to you, but I can guarantee that's in there.

Digging a little bit on my own I came to the conclusion that there ought to be something in the configuration of the GeoTIFF file that's not compatible with / adjusted to my polygon coordinates, most probably through the dataset.trasnform affine transformation matrix as seen here.

My file's aforementioned matrix, bounding box and coordinate reference system are as follows:

print(dataset.bounds)
# BoundingBox(left=160803.0, bottom=567978.0, right=187356.0, top=581073.0)

print(dataset.transform)
# | 3.00, 0.00, 160803.00|
# | 0.00,-3.00, 581073.00|
# | 0.00, 0.00, 1.00|

print(dataset.crs)
# +init=epsg:32637

Does anyone have any insight into why this might happen?

Best Answer

I managed to crop the GeoTiff file as intended.

The issue was that my GeoJSON / shape / polygon had its points in "normal" coordinate reference systems (which actually is epsg:4326 or WSG 84) while my geoTiff data is expressed in epsg:32637.

I used the following function to project my polygon:

def project_wsg_shape_to_csr(shape, csr):
  project = lambda x, y: pyproj.transform(
    pyproj.Proj(init='epsg:4326'),
    pyproj.Proj(init=csr),
    x,
    y
  )
  return shapely.ops.transform(project, shape)


projected_shape = project_wsg_shape_to_csr(shapely.geo.shape(polygon), 'epsg:32637')
Related Question