Crafting Shapefiles with Well-Structured Geometries

A complete test script is available below. To achieve the same outcome as the feature dict mentioned by the Fiona user manual, refer to the link provided: https://fiona.readthedocs.io/en/latest/manual.html#writing-vector-data. Even though the user has already accepted the GDAL/OGR response, a Fiona equivalent is also provided. It is important to note that Windows users have no reason not to use it.


Solution 1:

A reliable exchange format for binary data is the well-known binary format, which is compatible with various GIS software such as Shapely and GDAL/OGR.

Below is a small instance of the process involving

osgeo.ogr

, which is diminutive in size.

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


The GDAL/OGR answer has been accepted by the poster, but for the sake of providing an alternative, 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},
    })

(Attention Windows users: there are no reasons for your inaction)


Solution 2:


Fiona has been developed to seamlessly integrate with Shapely, enabling the user to easily enhance shapefile features. As a demonstration, here’s a basic example of how to use these tools in combination to “tidy up” 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'])

This code snippet is sourced from the GitHub repository of Fiona, located at https://github.com/Toblerity/Fiona/blob/master/examples/with-shapely.py.


Solution 3:


PyShp can be utilized to write Shapely geometries, which was also inquired by the original poster.

To write shapely geometry to a shapefile, you can first convert it to geojson using the shapely.geometry.mapping method. Then, use my PyShp fork that has a Writer method which can accept the geojson geometry dictionaries.

Below, I have included a conversion function for those who prefer to depend on the primary PyShp version.

# 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

To convert your shapely geometries into a pyshp compatible format, copy and paste the function into your script and call it as needed. After conversion, append each resulting pyshp shape to the ._shapes list of the shapefile.Writer instance to save them. For an example, refer to the test script provided at the end of this post.

Please note that the function does not address the presence of interior polygon holes; it disregards them completely. While it is feasible to incorporate this feature into the function, I have not yet taken the initiative to do so. Any recommendations or modifications to enhance the function are gladly accepted 🙂

Below is a comprehensive test script that can be used independently.

### 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")

Frequently Asked Questions