Konvertieren eines riesigen Multipolygons in Polygone

8

Ich habe ein Shapefile mit einigen riesigen Multipolygonen mit 100.000 Teilen. Was wäre der einfachste Weg, sie in einteilige Polygone aufzuteilen? Ich suche nach etwas wie der QGIS-Funktion "Multipart to Singlepart", aber die Datei ist zu groß, als dass QGIS sie verarbeiten könnte. Ich vermute, dass es wahrscheinlich bereits ein Python-Modul gibt, das dies für mich tun kann. Irgendwelche Tipps?

qgis python splitting multipart Löwe
quelle
Haben Sie die Möglichkeit, den Layer in PostGIS zu laden?
Holen Sie sich Spatial
In was ist deine Ebene jetzt gespeichert? Es könnte sein, dass das Ändern Ihres Speicherformats es ermöglicht, in QGIS zu arbeiten, wobei nur Effizienzunterschiede berücksichtigt werden.
Holen Sie sich Spatial
1
Es befindet sich in einem Shapefile, daher ist die unten stehende Lösung von @ Barrett perfekt!
Leo

Antworten:

11

Shapefiles haben keinen Typ MultiPolygon (Typ = Polygon), aber sie unterstützen sie trotzdem (alle Ringe werden in einem Polygon gespeichert = Liste der Polygone, siehe GDAL: ESRI Shapefile )

Mit Fiona und Shapely ist es einfacher :

import fiona
from shapely.geometry import shape, mapping

# open the original MultiPolygon file
with fiona.open('multipolygons.shp') as source:
    # create the new file: the driver and crs are the same
    # for the schema the geometry type is "Polygon" instead
    output_schema = dict(source.schema)  # make an independant copy
    output_schema['geometry'] = "Polygon"

    with fiona.open('output.shp', 'w', 
                    driver=source.driver,
                    crs=source.crs,
                    schema=output_schema) as output:

        # read the input file
        for multi in source:

           # extract each Polygon feature
           for poly in shape(multi['geometry']):

              # write the Polygon feature
              output.write({
                  'properties': multi['properties'],
                  'geometry': mapping(poly)
              })
Gen
quelle
11

von der GDAL-Mailingliste mit Python

import os
from osgeo import ogr

def multipoly2poly(in_lyr, out_lyr):
    for in_feat in in_lyr:
        geom = in_feat.GetGeometryRef()
        if geom.GetGeometryName() == 'MULTIPOLYGON':
            for geom_part in geom:
                addPolygon(geom_part.ExportToWkb(), out_lyr)
        else:
            addPolygon(geom.ExportToWkb(), out_lyr)

def addPolygon(simplePolygon, out_lyr):
    featureDefn = out_lyr.GetLayerDefn()
    polygon = ogr.CreateGeometryFromWkb(simplePolygon)
    out_feat = ogr.Feature(featureDefn)
    out_feat.SetGeometry(polygon)
    out_lyr.CreateFeature(out_feat)
    print 'Polygon added.'

from osgeo import gdal
gdal.UseExceptions()
driver = ogr.GetDriverByName('ESRI Shapefile')
in_ds = driver.Open('data/multipoly.shp', 0)
in_lyr = in_ds.GetLayer()
outputshp = 'data/poly.shp'
if os.path.exists(outputshp):
    driver.DeleteDataSource(outputshp)
out_ds = driver.CreateDataSource(outputshp)
out_lyr = out_ds.CreateLayer('poly', geom_type=ogr.wkbPolygon)
multipoly2poly(in_lyr, out_lyr)
Barrett
quelle