39

This question is one that would likely be closed quickly if asked today because it effectively asks for code without making a code attempt. Since it has many upvotes the community has assessed it to be useful, and so it is likely to be left open, but should not be used as an example of a question that should be asked today.

Can someone demonstrate a simple way to write geometry data-structures from shapely into shapefiles?

I am particularly interested in polygons with holes and linestrings. It would also be beneficial to stay away from ArcPy (so osgeo, pyshp, etc. would all be better).

5 Answers 5

65

Well-known binary is a good binary exchange format that can be exchanged with plenty of GIS software, including Shapely and GDAL/OGR.

This is a tiny example of the workflow with osgeo.ogr:

from osgeo import ogr
from shapely.geometry import Polygon

# Here's an example Shapely geometry
poly = Polygon([(0, 0), (0, 1), (1, 1), (0, 0)])

# Now convert it to a shapefile with OGR    
driver = ogr.GetDriverByName('Esri Shapefile')
ds = driver.CreateDataSource('my.shp')
layer = ds.CreateLayer('', None, ogr.wkbPolygon)
# Add one attribute
layer.CreateField(ogr.FieldDefn('id', ogr.OFTInteger))
defn = layer.GetLayerDefn()

## If there are multiple geometries, put the "for" loop here

# Create a new feature (attribute and geometry)
feat = ogr.Feature(defn)
feat.SetField('id', 123)

# Make a geometry, from Shapely object
geom = ogr.CreateGeometryFromWkb(poly.wkb)
feat.SetGeometry(geom)

layer.CreateFeature(feat)
feat = geom = None  # destroy these

# Save and close everything
ds = layer = feat = geom = None

Although the poster has accepted the GDAL/OGR answer, here is a Fiona equivalent:

from shapely.geometry import mapping, Polygon
import fiona

# Here's an example Shapely geometry
poly = Polygon([(0, 0), (0, 1), (1, 1), (0, 0)])

# Define a polygon feature geometry with one attribute
schema = {
    'geometry': 'Polygon',
    'properties': {'id': 'int'},
}

# Write a new Shapefile
with fiona.open('my_shp2.shp', 'w', 'ESRI Shapefile', schema) as c:
    ## If there are multiple geometries, put the "for" loop here
    c.write({
        'geometry': mapping(poly),
        'properties': {'id': 123},
    })
0
35

I've designed Fiona to work well with Shapely. Here is a very simple example of using them together to "clean" shapefile features:

import logging
import sys

from shapely.geometry import mapping, shape

import fiona

logging.basicConfig(stream=sys.stderr, level=logging.INFO)

with fiona.open('docs/data/test_uk.shp', 'r') as source:

    # **source.meta is a shortcut to get the crs, driver, and schema
    # keyword arguments from the source Collection.
    with fiona.open(
            'with-shapely.shp', 'w',
            **source.meta) as sink:

        for f in source:

            try:
                geom = shape(f['geometry'])
                if not geom.is_valid:
                    clean = geom.buffer(0.0)
                    assert clean.is_valid
                    assert clean.geom_type == 'Polygon'
                    geom = clean
                f['geometry'] = mapping(geom)
                sink.write(f)

            except Exception, e:
                # Writing uncleanable features to a different shapefile
                # is another option.
                logging.exception("Error cleaning feature %s:", f['id'])

From https://github.com/Toblerity/Fiona/blob/master/examples/with-shapely.py.

10

You can also write Shapely geometries by using PyShp (since the original poster also asked about PyShp).

One way would be to convert your shapely geometry to geojson (with the shapely.geometry.mapping method) and then use my modified fork of PyShp which provides a Writer method that accepts geojson geometry dictionaries when writing to a shapefile.

If you would rather rely on the main PyShp version, I have also provided a conversion function below:

# THIS FUNCTION CONVERTS A GEOJSON GEOMETRY DICTIONARY TO A PYSHP SHAPE OBJECT
def shapely_to_pyshp(shapelygeom):
    # first convert shapely to geojson
    try:
        shapelytogeojson = shapely.geometry.mapping
    except:
        import shapely.geometry
        shapelytogeojson = shapely.geometry.mapping
    geoj = shapelytogeojson(shapelygeom)
    # create empty pyshp shape
    record = shapefile._Shape()
    # set shapetype
    if geoj["type"] == "Null":
        pyshptype = 0
    elif geoj["type"] == "Point":
        pyshptype = 1
    elif geoj["type"] == "LineString":
        pyshptype = 3
    elif geoj["type"] == "Polygon":
        pyshptype = 5
    elif geoj["type"] == "MultiPoint":
        pyshptype = 8
    elif geoj["type"] == "MultiLineString":
        pyshptype = 3
    elif geoj["type"] == "MultiPolygon":
        pyshptype = 5
    record.shapeType = pyshptype
    # set points and parts
    if geoj["type"] == "Point":
        record.points = geoj["coordinates"]
        record.parts = [0]
    elif geoj["type"] in ("MultiPoint","Linestring"):
        record.points = geoj["coordinates"]
        record.parts = [0]
    elif geoj["type"] in ("Polygon"):
        record.points = geoj["coordinates"][0]
        record.parts = [0]
    elif geoj["type"] in ("MultiPolygon","MultiLineString"):
        index = 0
        points = []
        parts = []
        for eachmulti in geoj["coordinates"]:
            points.extend(eachmulti[0])
            parts.append(index)
            index += len(eachmulti[0])
        record.points = points
        record.parts = parts
    return record

Simply copy and paste the function to your own script and call on it to convert any of your shapely geometries to a pyshp compatible shape. To save them you then simply append each resulting pyshp shape to the shapefile.Writer instance's ._shapes list (for an example see the test script at the bottom of this post).

Note however: the function will NOT handle any interior polygon holes if there are any, it simply ignores them. It is certainly possible to add that functionality to the function but I have simply not bothered yet. Suggestions or edits to improve the function are welcome :)

Here is a full standalone test script:

### HOW TO SAVE SHAPEFILE FROM SHAPELY GEOMETRY USING PYSHP

# IMPORT STUFF
import shapefile
import shapely, shapely.geometry

# CREATE YOUR SHAPELY TEST INPUT
TEST_SHAPELYSHAPE = shapely.geometry.Polygon([(133,822),(422,644),(223,445),(921,154)])

#########################################################
################## END OF USER INPUT ####################
#########################################################

# DEFINE/COPY-PASTE THE SHAPELY-PYSHP CONVERSION FUNCTION
def shapely_to_pyshp(shapelygeom):
    # first convert shapely to geojson
    try:
        shapelytogeojson = shapely.geometry.mapping
    except:
        import shapely.geometry
        shapelytogeojson = shapely.geometry.mapping
    geoj = shapelytogeojson(shapelygeom)
    # create empty pyshp shape
    record = shapefile._Shape()
    # set shapetype
    if geoj["type"] == "Null":
        pyshptype = 0
    elif geoj["type"] == "Point":
        pyshptype = 1
    elif geoj["type"] == "LineString":
        pyshptype = 3
    elif geoj["type"] == "Polygon":
        pyshptype = 5
    elif geoj["type"] == "MultiPoint":
        pyshptype = 8
    elif geoj["type"] == "MultiLineString":
        pyshptype = 3
    elif geoj["type"] == "MultiPolygon":
        pyshptype = 5
    record.shapeType = pyshptype
    # set points and parts
    if geoj["type"] == "Point":
        record.points = geoj["coordinates"]
        record.parts = [0]
    elif geoj["type"] in ("MultiPoint","Linestring"):
        record.points = geoj["coordinates"]
        record.parts = [0]
    elif geoj["type"] in ("Polygon"):
        record.points = geoj["coordinates"][0]
        record.parts = [0]
    elif geoj["type"] in ("MultiPolygon","MultiLineString"):
        index = 0
        points = []
        parts = []
        for eachmulti in geoj["coordinates"]:
            points.extend(eachmulti[0])
            parts.append(index)
            index += len(eachmulti[0])
        record.points = points
        record.parts = parts
    return record

# WRITE TO SHAPEFILE USING PYSHP
shapewriter = shapefile.Writer()
shapewriter.field("field1")
# step1: convert shapely to pyshp using the function above
converted_shape = shapely_to_pyshp(TEST_SHAPELYSHAPE)
# step2: tell the writer to add the converted shape
shapewriter._shapes.append(converted_shape)
# add a list of attributes to go along with the shape
shapewriter.record(["empty record"])
# save it
shapewriter.save("test_shapelytopyshp.shp")
0
7

Karim's answer is pretty old but I have used his code. One minor thing I figured out using his code: If the shape type is Polygon or Multipolygon, it might still have multiple parts (holes inside). Therefore part of his code should be changed to

elif geoj["type"] == "Polygon":
    index = 0
    points = []
    parts = []
    for eachmulti in geoj["coordinates"]:
        points.extend(eachmulti)
        parts.append(index)
        index += len(eachmulti)
    record.points = points
    record.parts = parts
elif geoj["type"] in ("MultiPolygon", "MultiLineString"):
    index = 0
    points = []
    parts = []
    for polygon in geoj["coordinates"]:
        for part in polygon:
            points.extend(part)
            parts.append(index)
            index += len(part)
    record.points = points
    record.parts = parts
5

You can also use the geopandas lib for this purpose:

import geopandas as gpd
from shapely.geometry import Polygon

# Here's an example Shapely geometry
poly = Polygon([(0, 0), (0, 1), (1, 1), (0, 0)])

# use the feature loop in case you polygon is a multipolygon
features = [i for i in range(len(poly))]
# add crs using wkt or EPSG to have a .prj file
gdr = gpd.GeoDataFrame({'feature': features, 'geometry': poly}) #, crs='EPSG:4326)

gdr.to_file("Poly.shp")
3
  • Unfortunately this results in error. I found a fix those who are trying it out using geopandas at gis.stackexchange.com/a/395327 Oct 1, 2021 at 10:41
  • @AmanBagrecha I would be happy to edit my answer (that works fine on my side) if you can point me what is not working. Oct 3, 2021 at 16:07
  • The index for geodataframe was missing. The answer in the previous comment is what worked for me. Currently the "feature" key and value are not required I guess Oct 4, 2021 at 7:32

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.