1

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?

  • 1
    Just because you can plot two datasets as overlapping doesn't mean they share a coordinate reference -- It just means that they have correct metadata. Cropping an image requires a bounding box in the same CRS. GeoJSON is always geographic, so it is the wrong choice to clip a projected image. – Vince Nov 18 '18 at 13:15
3

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')

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.