Melarutkan poligon berdasarkan atribut dengan Python (rupawan, fiona)?

Saya telah mencoba untuk membuat fungsi yang pada dasarnya melakukan hal yang sama dengan fungsi "larut" QGIS. Saya pikir itu akan sangat mudah tetapi ternyata tidak. Jadi dari apa yang telah saya kumpulkan, penggunaan fiona dengan rupawan harus menjadi pilihan terbaik di sini. Saya baru saja mulai dipusingkan dengan file vektor sehingga dunia ini cukup baru bagi saya dan juga python.

Untuk contoh ini, saya bekerja dengan county shapefile yang didirikan di sini http://tinyurl.com/odfbanu Jadi di sini ada beberapa kode yang saya kumpulkan tetapi tidak dapat menemukan cara untuk membuatnya bekerja bersama

Untuk saat ini metode terbaik saya adalah sebagai berikut berdasarkan: https://sgillies.net/2009/01/27/a-more-perfect-union-continued.html . Ini berfungsi dengan baik dan saya mendapatkan daftar 52 negara sebagai geometri Shapely. Silakan berkomentar jika ada cara yang lebih mudah untuk melakukan bagian ini.

from osgeo import ogr
from shapely.wkb import loads
from numpy import asarray
from shapely.ops import cascaded_union

ds = ogr.Open('counties.shp')
layer = ds.GetLayer(0)

#create a list of unique states identifier to be able
#to loop through them later
STATEFP_list = []
for i in range(0 , layer.GetFeatureCount()) :
    feature = layer.GetFeature(i)
    statefp = feature.GetField('STATEFP')
    STATEFP_list.append(statefp)

STATEFP_list = set(STATEFP_list)

#Create a list of merged polygons = states 
#to be written to file
polygons = []

#do the actual dissolving based on STATEFP
#and append polygons
for i in STATEFP_list : 
    county_to_merge = []
    layer.SetAttributeFilter("STATEFP = '%s'" %i ) 
    #I am not too sure why "while 1" but it works 
    while 1:
        f = layer.GetNextFeature()
        if f is None: break
        g = f.geometry()
        county_to_merge.append(loads(g.ExportToWkb()))

    u = cascaded_union(county_to_merge)
    polygons.append(u)

#And now I am totally stuck, I have no idea how to write 
#this list of shapely geometry into a shapefile using the
#same properties that my source.

Jadi tulisannya benar-benar tidak lurus ke depan dari apa yang saya lihat, saya benar-benar hanya ingin bentuk yang sama hanya dengan negara larut menjadi negara, saya bahkan tidak perlu banyak tabel atribut tetapi saya ingin tahu bagaimana Anda dapat lulus pada dari sumber ke shapefile yang baru dibuat.

Saya menemukan banyak potongan kode untuk menulis dengan fiona tetapi saya tidak pernah bisa membuatnya berfungsi dengan data saya. Contoh dari Cara menulis geometri Shapely ke shapefile? :

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

Masalahnya di sini adalah bagaimana melakukan hal yang sama dengan daftar geometri dan cara membuat ulang properti yang sama dari sumbernya.

Pertanyaannya adalah tentang Fiona dan Shapely dan jawaban lain yang menggunakan GeoPandas juga perlu diketahui Panda . Apalagi GeoPandas menggunakan Fiona untuk membaca / menulis shapefile.

Saya tidak mempertanyakan utilitas GeoPandas di sini, tetapi Anda dapat melakukannya langsung dengan Fiona menggunakan modul standar itertools , khususnya dengan perintah groupby ("Singkatnya, groupby mengambil iterator dan memecahnya menjadi sub-iterator berdasarkan perubahan di "kunci" dari iterator utama. Ini tentu saja dilakukan tanpa membaca seluruh iterator sumber ke dalam memori ", itertools.groupby ).

Shapefile asli diwarnai oleh bidang STATEFP

from shapely.geometry import shape, mapping
from shapely.ops import unary_union
import fiona
import itertools
with fiona.open('cb_2013_us_county_20m.shp') as input:
    # preserve the schema of the original shapefile, including the crs
    meta = input.meta
    with fiona.open('dissolve.shp', 'w', **meta) as output:
        # groupby clusters consecutive elements of an iterable which have the same key so you must first sort the features by the 'STATEFP' field
        e = sorted(input, key=lambda k: k['properties']['STATEFP'])
        # group by the 'STATEFP' field 
        for key, group in itertools.groupby(e, key=lambda x:x['properties']['STATEFP']):
            properties, geom = zip(*[(feature['properties'],shape(feature['geometry'])) for feature in group])
            # write the feature, computing the unary_union of the elements in the group with the properties of the first element in the group
            output.write({'geometry': mapping(unary_union(geom)), 'properties': properties[0]})

Hasil

Saya sangat merekomendasikan GeoPandas untuk berurusan dengan berbagai macam fitur dan melakukan operasi massal.

Itu memperluas bingkai data Pandas, dan menggunakan Shapely di bawah tenda.

from geopandas import GeoSeries, GeoDataFrame

# define your directories and file names
dir_input = '/path/to/your/file/'
name_in = 'cb_2013_us_county_20m.shp'
dir_output = '/path/to/your/file/'
name_out = 'states.shp'

# create a dictionary
states = {}
# open your file with geopandas
counties = GeoDataFrame.from_file(dir_input + name_in)

for i in range(len(counties)):
    state_id = counties.at[i, 'STATEFP']
    county_geometry = counties.at[i, 'geometry']
    # if the feature's state doesn't yet exist, create it and assign a list
    if state_id not in states:
        states[state_id] = []
    # append the feature to the list of features
    states[state_id].append(county_geometry)

# create a geopandas geodataframe, with columns for state and geometry
states_dissolved = GeoDataFrame(columns=['state', 'geometry'], crs=counties.crs)

# iterate your dictionary
for state, county_list in states.items():
    # create a geoseries from the list of features
    geometry = GeoSeries(county_list)
    # use unary_union to join them, thus returning polygon or multi-polygon
    geometry = geometry.unary_union
    # set your state and geometry values
    states_dissolved.set_value(state, 'state', state)
    states_dissolved.set_value(state, 'geometry', geometry)

# save to file
states_dissolved.to_file(dir_output + name_out, driver="ESRI Shapefile")

Sebagai tambahan untuk jawaban gen @ , saya perlu membubarkan lebih dari satu bidang sehingga memodifikasi kodenya untuk menangani beberapa bidang. Kode di bawah ini digunakan operator.itemgetteruntuk mengelompokkan berdasarkan beberapa bidang:

# Modified from /gis//a/150001/2856
from shapely.geometry import shape, mapping
from shapely.ops import unary_union
import fiona
import itertools
from operator import itemgetter


def dissolve(input, output, fields):
    with fiona.open(input) as input:
        with fiona.open(output, 'w', **input.meta) as output:
            grouper = itemgetter(*fields)
            key = lambda k: grouper(k['properties'])
            for k, group in itertools.groupby(sorted(input, key=key), key):
                properties, geom = zip(*[(feature['properties'], shape(feature['geometry'])) for feature in group])
                output.write({'geometry': mapping(unary_union(geom)), 'properties': properties[0]})


if __name__ == '__main__':
    dissolve('input.shp', 'input_dissolved.shp', ["FIELD1", "FIELD2", "FIELDN"))