Friday, 15 December 2017

Masking GeoTIFF file after GeoJSON through rasterio - "Input shapes do not overlap raster"


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?



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

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...