Python Shapely – How to Get Polygon Breadth

pythonshapely

Is it possible to get polygon breadth in Shapely, similar to how it provides polygon length?
http://toblerity.org/shapely/shapely.geometry.html

I am also ok with a solution that uses fiona or any non-arcgis open source python solution.

Best Answer

# create example polygon
poly = Polygon([(0, 0), (4, 0), (5, 2), (7, 5), (3, 2), (1, 3)])

# get minimum bounding box around polygon
box = poly.minimum_rotated_rectangle

# get coordinates of polygon vertices
x, y = box.exterior.coords.xy

# get length of bounding box edges
edge_length = (Point(x[0], y[0]).distance(Point(x[1], y[1])), Point(x[1], y[1]).distance(Point(x[2], y[2])))

# get length of polygon as the longest edge of the bounding box
length = max(edge_length)

# get width of polygon as the shortest edge of the bounding box
width = min(edge_length)

Related Solutions

Shapely Python – Looping Through a Number of Polygons

Some solutions according to the position of an element in a list:

polA = Polygon([(0,0), (3,0), (3,3), (0,3)])
polB = Polygon([(2,-1), (5,-1), (5,2), (2,2)])
polC = Polygon([(5,2), (8,2), (8,5), (5,5)]) 
collection = [polA, polB, polC]

Iterating by index:

for i in range(len(collection)-1):
    print collection[i], collection[i+1], collection[i].touches(collection[i+1])
POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) False
POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) True

Iterating with itertools (standard module), all combinations without repeated elements

import itertools
for pol in  itertools.combinations(collection, 2):
    print pol[0],pol[1], pol[0].touches(pol[1])
 POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) False
 POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) False
 POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) True

or with the suggestion of Mike T

for pol1, pol2 in  itertools.combinations(collection, 2):
     print pol1, pol2, pol1.touches(pol2)
 POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) False
 POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) False
 POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) True

Compare with

for pol in  itertools.permutations(collection, 2):
    print pol[0],pol[1], pol[0].touches(pol[1]) 
POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) False
POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) False
POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) False
POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) True
POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) False
POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) True

Look also at Python - Previous and next values inside a loop

from itertools import tee, islice, chain, izip
def previous_and_next(some_iterable):
    prevs, items, nexts = tee(some_iterable, 3)
    prevs = chain([None], prevs)
    nexts = chain(islice(nexts, 1, None), [None])
    return izip(prevs, items, nexts)
 for previous, item, nxt in previous_and_next(collection):
      print "Item is now", item, "next is", nxt, "previous is", previous

 Item is now POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) next is POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) previous is None
 Item is now POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) next is POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) previous is POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0))
 Item is now POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) next is None previous is POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1))

[GIS] Convert shapely polygon coordinates

To change projections with Fiona, use the pyproj module.

Example with a point shapefile (you can simplify the algorithm):

from pyproj import Proj, transform
import fiona
from fiona.crs import from_epsg
shape = fiona.open('sample.shp')
original = Proj(shape.crs) # EPSG:4326 in your case
destination = Proj(init='EPSG:...') # your new EPSG
with fiona.open('new.shp', 'w', 'ESRI Shapefile', shape.schema.copy(), crs=from_epsg(...)) as output:
    for feat in shape:
        long,lat =  feat['geometry']['coordinates']
        x,y = transform(original, destination,long,lat)
        # change only the coordinates of the feature
        feat['geometry']['coordinates'] = (x,y)
        output.write(feat)

If you want to work with the GeoJSON format:

with fiona.open('sample.shp') as source:
    records = list(source)
geojson = {"type": "FeatureCollection","features": records}

And use the same method with the GeoJSON format (with pyproj)

For a polygon, do the same thing with the points of the polygon

for feat in shape:
   for point in feat['geometry']['coordinates'][0]: # = LinearRing of the polygon
       x,y = pt
       ....

new

In your script, why counter?
- feat['geometry']['coordinates'] is a nested list (Polygon)
- feat['geometry']['coordinates'][0] also (LinearRing of the Polygon)
- feat['geometry']['coordinates'][0][x] is a Point of the LinearRing

So:

with fiona.open('new.shp', 'w', 'ESRI Shapefile', shape.schema.copy(), crs=from_epsg(...)) as output:
    for feat in shape: # feat = one polygon of the shapefile
        out_linearRing = [] # empty list for the LinearRing of transformed coordinates
        for point in feat['geometry']['coordinates'][0]: # LinearRing of the Polygon
            long,lat =  point  # one point of the LinearRing
            x,y = transform(original, destination,long,lat) # transform the point
            out_linearRing.append((x,y)) # add all the points to the new LinearRing
        # transform the resulting LinearRing to  a Polygon and write it
        feat['geometry']['coordinates'] = [out_linearRing]
        output.write(feat)

Best Answer

Related Solutions

Shapely Python – Looping Through a Number of Polygons

[GIS] Convert shapely polygon coordinates

new

Related Question