[GIS] How to return WKT of entire polygon with shapely

pythonshapefileshapelywell-known-text

Im learning to use shapely and have two shapefiles–forest, and basins. The two shapefiles clearly intersect (basins is completely within forest). But the way I've gone about it only gets the feature at index 0 of each shapefile, which don't intersect. I know where this is taking place in the script, but how do I get the whole polygon shape as wkt instead of just one record? Do I build a list of the wkt geometry one record at a time?

# Basins shapefile to wkt
basins = ogr.Open("data/Basins.shp")
basins_layer = basins.GetLayer(0)
basins_feature = basins_layer.GetFeature(0) # This is the line
# I'm thinking there should probably be a loop here to get the geometry of the whole shapefile
basins_geometry = basins_feature.GetGeometryRef()
basins_wkt = shapely.wkt.loads(basins_geometry.ExportToWkt())

How do I loop through the whole shapefile and build the outer polygon as wkt or is there another function to get the outer boundary of the shapefile? I guess this is called the outer ring. I want wkt of the whole shapefile exterior, rather than the individual polygons that make it up.

Best Answer

To get wkt outer ring of a polygon with shapely, you can use next code with your particular paths:

import fiona
from shapely.geometry import shape, LineString

path = '/home/zeito/pyqgis_data/polygon1.shp' #polygon with only one feature

c = fiona.open(path)

collection = [ shape(item['geometry']) for item in c ]

rings = [ LineString(pol.coords).wkt for pol in collection ]

print rings[0] #0 because polygon has only one feature (one outer ring)

After running above code at Python Console, for my particular path, I got:

LINESTRING (389535.8208391897 4448641.082016046, 397951.9055779715 4459595.351041127, 418925.3231016023 4456522.812168239, 425070.4008473795 4441160.117803795, 412513.0680625305 4427801.253139063, 392341.1824187837 4435148.628704665, 389535.8208391897 4448641.082016046)

By using QuickWKT plugin of QGIS, it can be observed that it works as expected:

enter image description here

Editing Note:

If you have issues to install fiona, an alternative code, by using ogr python module, is the following:

from osgeo import ogr
from shapely.wkt import loads
from shapely.geometry import LineString

path = '/home/zeito/pyqgis_data/polygon1.shp'

basins = ogr.Open(path)

layer = basins.GetLayer()

geoms = []

for feature in layer:
    geom = feature.GetGeometryRef()
    geoms.append(geom.ExportToWkt())

pol = loads(geoms[0])

print LineString(pol.exterior.coords).wkt

It produces same result that first code.

Related Question