1

I have a TIFF file (Landsat8 tile) and a polygon of a burned area, I would like to plot only the polygon appropriately within the bounds of the TIFF and save that as an image however I am struggling.

The polygon has been uploaded as a GitHub gist here. The image TIFF can be found here.

Plotting the TIFF with rasterio and matplotlib:

with rasterio.open(f'{input_dir}/{band2}', 'r', height=1024, width=1024) as t:
    f, ax = plt.subplots(figsize=(10,10))
    show(t.read(1), ax=ax)
    print(t.crs)

Landsat 8 tile in TIFF format

The CRS of this image is: EPSG:32610

Plotting the polygon produces the following:

f, ax = plt.subplots()
gpd.GeoSeries(boundary[0], crs='EPSG:32610').plot(ax=ax, color='red')

Polygon plot

The CRS doesn't seem to be correct because this polygon is now way larger than the image it should be in the bounds of. When plotting the polygon directly onto the tile (just to see what that would look like) produces the following:

with rasterio.open(f'{input_dir}/{band2}', 'r', height=1024, width=1024) as t:
    f, ax = plt.subplots(figsize=(10,10))
    show(t.read(1), ax=ax)
    gpd.GeoSeries(boundary[0], crs='EPSG:32610').plot(ax=ax, color='red')

Polygon plotted over the tile

The end result I am hoping for is just the polygon where it should exist in the bounds of the tile on a blank background.

Edit The original CRS of the polygon was EPSG:4326. However if I then try the following:

gpd.GeoSeries(boundary[0], crs='EPSG:4326').to_crs(epsg=32610)

I get the polygon as seen above.

0    POLYGON ((515747.044 4334707.851, 515863.285 4...
dtype: geometry

I've uploaded the original polygon in the EPSG:4326 here.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.