Shapely – How to Convert Polygons to Lines

shapely

Given a shapefile with polygons, how do I convert the polygons to individual lines instead? I know how to do this in QGIS, but I'm looking for a Shapely function which does the same thing.

Best Answer

bugmenot123 is ok but I find easier use the boundary of the polygons. If you have a multipolygon or a polygon with holes the boundary returns a multilinestrig, if you have a polygon without holes the boundary returns a linestring.

Here is a simple example so you can see how it works:

import shapely
from shapely.geometry import MultiPolygon, Point

pol1 = MultiPolygon([Point(0, 0).buffer(2.0), Point(1, 1).buffer(2.0)])
pol2 = Point(7, 8).buffer(1.0)
pols = [pol1, pol2]

lines = []
for pol in pols:
    boundary = pol.boundary
    if boundary.type == 'MultiLineString':
        for line in boundary:
            lines.append(line)
    else:
        lines.append(boundary)

for line in lines:
    print line.wkt

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