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).
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},
})
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.
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")
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")
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
