Mengubah multipoligon besar menjadi poligon

9

Saya memiliki shapefile dengan beberapa multipoligon besar, dengan 100.000-an bagian. Apa cara termudah untuk membaginya menjadi poligon tunggal? Saya mencari fungsi seperti QGIS "Multipart to singlepart", tetapi file tersebut terlalu besar untuk ditangani QGIS. Saya menduga mungkin sudah ada beberapa modul Python yang dapat melakukannya untuk saya. Ada tips?

Leo
sumber
Apakah Anda memiliki pilihan untuk memuat layer ke PostGIS?
Dapatkan Spasial
Apa layer Anda sekarang? Bisa jadi mengubah format penyimpanan Anda akan memungkinkannya bekerja di QGIS, hanya memperhitungkan perbedaan efisiensi.
Dapatkan Spasial
1
Ada dalam shapefile, jadi solusi @ Barrett di bawah ini sempurna!
leo

Jawaban:

11

Shapefile tidak memiliki tipe MultiPolygon (type = Polygon), tetapi mereka tetap mendukungnya (semua cincin disimpan dalam satu poligon = daftar poligon, lihat GDAL: ESRI Shapefile )

Lebih mudah dengan Fiona dan Shapely :

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

dari milis GDAL menggunakan 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
sumber