Bagaimana cara menulis geometri Shapely ke shapefile?

26

Dapatkah seseorang mendemonstrasikan cara sederhana untuk menulis struktur data geometri dari shapely menjadi shapefile? Saya khususnya tertarik pada poligon dengan lubang dan linestrings. Ini juga akan bermanfaat untuk menjauh dari arcpy (jadi osgeo, pyshp, dll. Semuanya akan lebih baik).

terra_matics
sumber

Jawaban:

44

Biner terkenal adalah format pertukaran biner yang baik yang dapat ditukar dengan banyak perangkat lunak GIS, termasuk Shapely dan GDAL / OGR.

Ini adalah contoh kecil dari alur kerja dengan 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

Pembaruan : Meskipun poster telah menerima jawaban GDAL / OGR, berikut ini adalah padanan 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},
    })

(Catatan pengguna Windows: Anda tidak punya alasan )

Mike T
sumber
Tertarik mengapa Anda memilih metode ini daripada perpustakaan Fiona.
Nathan W
1
Nah, poster itu sedang mencari contoh osgeo.ogr, dan perbandingannya menarik.
sgillies
@sgillies perbandingan eksplisit ditambahkan
Mike T
3
Yah, sejujurnya itu sebagian besar pragmatik. Saya menghargai upaya menunjukkan kode dalam menanggapi pertanyaan saya dan saya sudah mucking dengan osgeo. Saya telah mencoba kedua metode dan keduanya merupakan jawaban yang cukup. Saya menghargai upaya dari para responden yang akurat dan cepat.
terra_matics
@ Mike T Mengenai pendekatan osgeo.ogr, saya menggunakannya dalam plugin Python untuk QGIS. Shapefile yang dipertimbangkan untuk ditulis adalah Garis (LineString in Shapely). Di mana Anda telah mendefinisikan variabel "poli", saya telah mendefinisikan variabel "baris" dengan koordinat dari Qgs.Rectangle. Saya telah menggunakan kode yang tepat, tidak ada kesalahan, tetapi tidak menambah fitur, dan memberi saya sebuah shapefile tanpa fitur.
Akhil
28

Saya telah merancang Fiona agar bekerja dengan baik dengan Shapely. Berikut adalah contoh yang sangat sederhana untuk menggunakannya bersama-sama untuk "membersihkan" fitur 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'])

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

sgillies
sumber
6

Anda juga dapat menulis geometri Shapely dengan menggunakan PyShp (karena poster asli juga bertanya tentang PyShp).

Salah satu caranya adalah dengan mengubah geometri bentuk Anda menjadi geojson (dengan metode shapely.geometry.mapping) dan kemudian menggunakan garpu PyShp saya yang dimodifikasi yang menyediakan metode Writer yang menerima kamus geometri geojson ketika menulis ke shapefile.

Jika Anda lebih suka mengandalkan versi PyShp utama, saya juga menyediakan fungsi konversi di bawah:

# 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

Cukup salin dan tempel fungsi ke skrip Anda sendiri dan panggil untuk mengonversi salah satu geometri bentuk Anda ke bentuk pyshp yang kompatibel. Untuk menyimpannya, Anda cukup menambahkan setiap bentuk pyshp yang dihasilkan ke daftar ._shapes dari instance shapefile.Writer (sebagai contoh lihat skrip uji di bagian bawah posting ini).

Namun, perhatikan: fungsinya TIDAK akan menangani lubang poligon interior jika ada, hanya mengabaikannya. Sangat mungkin untuk menambahkan fungsionalitas itu ke fungsi tetapi saya belum merasa terganggu. Saran atau suntingan untuk meningkatkan fungsi dipersilahkan :)

Berikut ini skrip uji mandiri lengkap:

### 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")
Karim Bahgat
sumber
5

Jawaban Karim cukup lama tetapi saya telah menggunakan kodenya dan ingin mengucapkan terima kasih kepadanya. Satu hal kecil yang saya tahu menggunakan kode-nya: Jika tipe bentuknya Polygon atau Multipolygon, mungkin masih ada banyak bagian (lubang di dalamnya). Oleh karena itu bagian dari kodenya harus diubah menjadi

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
lebenlechzer
sumber