如何将Shapely几何体写入shapefile？

有人可以演示一种简单的方法将几何数据结构从形状写入shapefile吗？我对带有孔和线串的多边形特别感兴趣。远离arcpy也将是有益的（因此osgeo，pyshp等都会更好）。

#1 楼

众所周知的二进制文件是一种很好的二进制交换格式，可以与许多GIS软件（包括Shapely和GDAL / OGR）进行交换。

这是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

更新：尽管张贴者接受了GDAL / OGR的答案，但这是Fiona等效项：

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

（请注意Windows用户：您没有任何借口）

有趣的是为什么选择了此方法而不是Fiona库。

–内森·W
13年2月24日在7:21

好吧，发布者正在寻找osgeo.ogr示例，并且比较很有趣。

–傻瓜
13年2月24日在17:08

添加了@sgillies显式比较

– Mike T
13年2月24日在18:33

好吧，说实话这主要是实用主义的。我很感谢为回答我的问题而演示代码的工作，并且我已经对osgeo产生了兴趣。从那以后，我尝试了两种方法，它们都是足够的答案。我感谢响应者所做的努力，既准确又快速。

–terra_matics
13年2月24日在22:06

@Mike T关于osgeo.ogr方法，我在QGIS的Python插件中使用它。考虑写入的shapefile是Line（Shapely中的LineString）。在定义了“ poly”变量的地方，我定义了一个“ line”变量，其坐标来自Qgs.Rectangle。我使用了确切的代码，没有错误，但是它没有添加功能，并且给了我一个没有功能的shapefile。

–阿克希尔
2015年1月9日在11:06

#2 楼

我设计的Fiona与Shapely可以很好地协作。这是一个非常简单的示例，将它们一起使用以“清理” shapefile功能：

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'])

来自https://github.com/Toblerity/Fiona/blob/master /examples/with-shapely.py。

#3 楼

您也可以使用PyShp（因为原始海报也询问了PyShp）来编写Shapely几何图形。然后使用我修改过的PyShp分支，该分支提供了Writer方法，该方法在写入shapefile时接受geojson几何字典。

如果您想依赖主要的PyShp版本，我还在下面提供了转换功能：函数执行您自己的脚本，然后调用它以将任何形状的几何体转换为pyshp兼容的形状。要保存它们，您只需将每个生成的pyshp形状附加到shapefile.Writer实例的._shapes列表中（有关示例，请参见本文底部的测试脚本）。

但是请注意：该函数将不会处理任何内部多边形孔（如果有的话），它只会忽略它们。当然可以将该功能添加到功能中，但是我还没有打扰。欢迎提出建议或进行修改以改善功能：)

以下是完整的独立测试脚本：

# 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

#4 楼

Karim的答案很老，但我使用了他的代码，并想感谢他。我用他的代码弄清楚了一件小事：如果形状类型是Polygon或Multipolygon，它可能仍然具有多个部分（内部有孔）。因此，他的部分代码应更改为

 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 楼

我在这里发布（针对任何有兴趣的人），您也可以出于以下目的使用geopandas库：

编程黑洞网