0

I have a shapefile that contains multiple polygons. I would like to extract the extent of each of those polygons within that shapefile. The majority of answers that I have been able to find use ArcPy and I do not have ArcMap or ArcGIS Pro, so I cannot use it. So far I have code that will obtain the extent of the shapefile, but not each individual polygon within the shapefile. I am using fiona and shapely.

from shapely.geometry import shape
import fiona

path = 'D:/data/shapefile/'
file = '__932.shp'

c = fiona.open(path + file)
field = next(iter(c)
extent = shape(field['geometry']).bounds
print(extent)

# (465452.3213890702, 4924992.33678632, 465453.50800855865, 4924993.2240787465)

2 Answers 2

6

There are many solutions

With fiona

from shapely.geometry import shape
import fiona
layer = fiona.open("test.shp")
# now layer is a generator/iterator so you can use next()
shape(next(layer)['geometry']).bounds # gives the bound of the first feature
(204386.82131440676, 89543.4759701695, 204422.20968847454, 89578.86434423728)
shape(next(layer)['geometry']).bounds # gives the bound of the second feature
(204351.43294033897, 89508.0875961017, 204386.82131440676, 89543.47597016949)
.....

But it's easier with a for loop:

for geom in layer:
    print(shape(geom['geometry']).bounds)
(204386.82131440676, 89543.4759701695, 204422.20968847454, 89578.86434423728)
(204351.43294033897, 89508.0875961017, 204386.82131440676, 89543.47597016949)
(204386.82131440676, 89508.0875961017, 204422.20968847454, 89543.47597016949)
(204351.43294033897, 89472.69922203392, 204386.82131440676, 89508.0875961017)

With GeoPandas

import geopandas as gpd
df = gpd.read_file("test.shp")
print(df.geometry.bounds)
        minx          miny           maxx          maxy
204386.821314  89543.475970  204422.209688  89578.864344
204351.432940  89508.087596  204386.821314  89543.475970
204386.821314  89508.087596  204422.209688  89543.475970
204351.432940  89472.699222  204386.821314  89508.087596

print(df.total_bounds) # bounds of the whole shapefile
[204351.43294034  89472.69922203 204422.20968847  89578.86434424]

With PyShp (in pure Python)

import shapefile
reader = shapefile.Reader("test.shp")
for geom in reader.shapeRecords():
    print(shape(geom.shape.__geo_interface__).bounds)
(204386.82131440676, 89543.4759701695, 204422.20968847454, 89578.86434423728)
(204351.43294033897, 89508.0875961017, 204386.82131440676, 89543.47597016949)
(204386.82131440676, 89508.0875961017, 204422.20968847454, 89543.47597016949)
(204351.43294033897, 89472.69922203392, 204386.82131440676, 89508.0875961017)

Control with your solution

shapefile = ogr.Open("test.shp")
layer = shapefile.GetLayer(0)
for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    geometry = feature.GetGeometryRef()
    extent = geometry.GetEnvelope()
    print(extent)
(204386.82131440676, 204422.20968847454, 89543.4759701695, 89578.86434423728)
(204351.43294033897, 204386.82131440676, 89508.0875961017, 89543.47597016949)
(204386.82131440676, 204422.20968847454, 89508.0875961017, 89543.47597016949)
(204351.43294033897, 204386.82131440676, 89472.69922203392, 89508.0875961017)
0
2

I found a solution from Get feature extent using gdal/ogr at Stack Exchange using the osgeo library

from osgeo import ogr

path = 'D:/data/shapefile/'
file = '__932.shp'

shapefile = ogr.Open(path + file)
layer = shapefile.GetLayer(0)
for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    geometry = feature.GetGeometryRef()
    extent = geometry.GetEnvelope()
    print(extent)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.