3

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.

  • 1
    If you're using the ogr module you don't need shapely. All of the functionality of shapely is present in ogr. – sgillies Jul 24 '17 at 15:09
4

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.exterior.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.

  • This is very close to what I suggest: iterate over all records in the shapefile as above, get the individual polygon geometries, and make a multipolygon geometry from them. You can then compare this multipolygon to other geometries or dissolve it into a single polygon using the unary_union function. – sgillies Jul 24 '17 at 15:13
  • Thanks. This looks like a good solution and I've been trying all day to get Fiona to no avail. I've been able to successfully install it using pip, but can't import it in python. – Alex R Jul 24 '17 at 18:56
  • As you have issues to install fiona, I modified my first code to obtain same result (see my Editing Note) but, by using org python module. – xunilk Jul 24 '17 at 22:42
0

You can just use .wkt as wkt_polygon = ShapelyPolygon(..).wkt

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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