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:
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))
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?
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')
No comments:
Post a Comment