Dissoudre des polygones basés sur des attributs avec Python (galbe, fiona)?

J'ai essayé de créer une fonction qui fait essentiellement la même chose que la fonction "dissoudre" de QGIS. Je pensais que ce serait super facile, mais apparemment non. Donc, d'après ce que j'ai rassemblé, l'utilisation de fiona avec galbe devrait être la meilleure option ici. Je viens de commencer à jouer avec les fichiers vectoriels, donc ce monde est assez nouveau pour moi et python aussi.

Pour ces exemples, je travaille avec un fichier de formes de comté fondé ici http://tinyurl.com/odfbanu Voici donc un morceau de code que j'ai rassemblé mais ne trouve pas de moyen de les faire fonctionner ensemble

Pour l'instant ma meilleure méthode est la suivante basée sur: https://sgillies.net/2009/01/27/a-more-perfect-union-continued.html . Cela fonctionne très bien et j'obtiens une liste des 52 états en tant que géométrie galbée. N'hésitez pas à commenter s'il existe un moyen plus simple de faire cette partie.

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.

Donc, l'écriture n'est vraiment pas simple de ce que j'ai vu, je veux vraiment le même fichier de forme uniquement avec le pays se dissolvant en états, je n'ai même pas besoin de beaucoup de la table d'attributs mais je suis curieux de voir comment vous pouvez passer depuis la source vers le nouveau fichier de formes créé.

J'ai trouvé de nombreux morceaux de code pour écrire avec fiona mais je ne suis jamais en mesure de le faire fonctionner avec mes données. Exemple de Comment écrire des géométries galbées dans des fichiers de formes? :

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

Le problème ici est de savoir comment faire de même avec une liste de géométrie et comment recréer les mêmes propriétés que la source.

La question concerne Fiona et Shapely et l'autre réponse utilisant GeoPandas nécessite de connaître également les Pandas . De plus, GeoPandas utilise Fiona pour lire / écrire des fichiers de formes.

Je ne remets pas en question ici l'utilité de GeoPandas, mais vous pouvez le faire directement avec Fiona en utilisant le module standard itertools , spécialement avec la commande groupby ("En bref, groupby prend un itérateur et le décompose en sous-itérateurs en fonction des changements dans la "clé" de l'itérateur principal. Cela se fait bien sûr sans lire l'intégralité de l'itérateur source en mémoire ", itertools.groupby ).

Shapefile original coloré par le champ 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]})

Résultat

Je recommande fortement GeoPandas pour gérer de nombreux assortiments de fonctionnalités et effectuer des opérations en masse.

Il étend les cadres de données Pandas et utilise une forme galbée sous le capot.

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")

En complément de la réponse de @ gene , j'avais besoin de dissoudre par plus d'un champ, j'ai donc modifié son code pour gérer plusieurs champs. Le code ci-dessous permet operator.itemgetterde regrouper par plusieurs champs:

# 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"))