16

I have a raster that I opened with rasterio. When I do raster.bounds I get BoundingBox(left=399960.0, bottom=-3309780.0, right=609780.0, top=-3099960.0). My question is now, if there is any simple way to make a shapefile of it?

Vince
  • 20,017
  • 15
  • 45
  • 64
Robin Kohrs
  • 675
  • 1
  • 8
  • 21
  • the thing is I have a raster with a bounding box. I don't want to create a bounding box using coordinates – Robin Kohrs Feb 29 '20 at 13:34
  • Please don't comment on your own question. Instead, [Edit] the Question to contain this information. One in mind that raster boundaries are sometimes computed from center-center of the corner pixels, so your polygon may be a half-pixel too small in each corner. – Vince Feb 29 '20 at 13:58
  • thank you very much for the advices:)! – Robin Kohrs Feb 29 '20 at 14:17

4 Answers4

26

With shapely box and GeoPandas

import rasterio as rio
ra = rio.open("raster.tif")
bounds  = ra.bounds

Convert bounds to shapely geometry

from shapely.geometry import box
geom = box(*bounds)
print(geom.wkt)
POLYGON ((205625.4144070312 88411.04799999669, 205625.4144070312 90534.35044407193, 202086.5770002506 90534.35044407193, 202086.5770002506 88411.04799999669, 205625.4144070312 88411.04799999669))

Create a shapefile with GeoPandas

df = gpd.GeoDataFrame({"id":1,"geometry":[geom]})
df.to_file("boundary.shp")

In one line:

df = gpd.GeoDataFrame({"id":1,"geometry":[box(*bounds)]})

Result:

enter image description here

gene
  • 54,868
  • 3
  • 110
  • 187
6

For scripting, other answers does the job. If you need generating shp from command line, I would recommand to use rasterio command line utility rio with

rio bounds /tmp/RGBA.byte.tif > tmp.geojson
ogr2ogr RGBA.shp tmp.geojson
rm tmp.geojson

You could also use gdaltindex e.g

gdaltindex RGBA.shp RGBA.byte.tif
ThomasG77
  • 30,725
  • 1
  • 53
  • 93
5

Here is a solution using shapely and fiona that takes into account the spatial reference of the raster and makes sure the new shapefile has it:

# imports
from shapely.geometry import box, mapping
import fiona
import rasterio

# create a Polygon from the raster bounds
bbox = box(*raster.bounds)

# create a schema with no properties
schema = {'geometry': 'Polygon', 'properties': {}}

# create shapefile
with fiona.open('bbox.shp', 'w', driver='ESRI Shapefile',
                crs=raster.crs.to_dict(), schema=schema) as c:
    c.write({'geometry': mapping(bbox), 'properties': {}})
Marcelo Villa
  • 5,928
  • 2
  • 19
  • 38
3

There is also a GDAL program gdaltindex https://gdal.org/programs/gdaltindex.html for this purpose.

Description

This program builds a shapefile with a record for each input raster file, an attribute containing the filename, and a polygon geometry outlining the raster.

user30184
  • 65,331
  • 4
  • 65
  • 118