¿Cómo escribir geometrías de Shapely en shapefiles?

22

¿Puede alguien demostrar una forma sencilla de escribir estructuras de datos de geometría desde bien formados en shapefiles? Estoy particularmente interesado en polígonos con agujeros y cadenas. También sería beneficioso mantenerse alejado de arcpy (por lo que osgeo, pyshp, etc. serían mejores).

    
pregunta terra_matics 23.02.2013 - 05:24

4 respuestas

35

Binario conocido es un buen formato de intercambio binario que se puede intercambiar con un montón de software GIS, incluido Shapely y GDAL / OGR.

Este es un ejemplo pequeño del flujo de trabajo con 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

Actualización : Aunque el póster ha aceptado la respuesta GDAL / OGR, aquí hay un Fiona equivalente :

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

(Tenga en cuenta a los usuarios de Windows: no tiene excusa )

    
respondido por el Mike T 23.02.2013 - 07:35
27

He diseñado Fiona para que funcione bien con Shapely. Este es un ejemplo muy simple de usarlos juntos para "limpiar" las características del 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'])

De enlace .

    
respondido por el sgillies 23.02.2013 - 06:05
5

También puedes escribir geometrías de Shapely usando PyShp (ya que el póster original también preguntó sobre PyShp).

Una forma sería convertir su geometría bien proporcionada a geojson (con el método shapely.geometry.mapping) y luego usar mi bifurcación modificada de PyShp que proporciona un método Writer que acepta los diccionarios de geometría geojson al escribir en un shapefile.

Si prefiere confiar en la versión principal de PyShp, también proporcioné una función de conversión a continuación:

# 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

Simplemente copie y pegue la función en su propio script y solicítelo para convertir cualquiera de sus geometrías bien definidas a una forma compatible con pyshp. Para guardarlos, simplemente agregue cada forma de pyshp resultante a la lista de formas de forma del archivo shape.Writer (para ver un ejemplo, vea el script de prueba al final de esta publicación).

Sin embargo, tenga en cuenta que la función NO manejará ningún agujero interior de polígono si hay alguno, simplemente los ignora. Ciertamente es posible agregar esa funcionalidad a la función, pero simplemente no me he molestado todavía. Las sugerencias o ediciones para mejorar la función son bienvenidas :)

Aquí hay un script de prueba completamente independiente:

### 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")
    
respondido por el Karim Bahgat 17.02.2014 - 14:29
3

La respuesta de Karim es bastante antigua, pero he usado su código y quería agradecerle por ello. Una cosa menor que descubrí usando su código: si el tipo de forma es Polígono o Multipolígono, es posible que aún tenga varias partes (agujeros en el interior). Por lo tanto parte de su código debe ser cambiado a

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
    
respondido por el lebenlechzer 13.12.2017 - 16:06

Lea otras preguntas en las etiquetas