Dissolvendo polígonos com base em atributos com Python (bem torneado, fiona)?

Eu tenho tentado criar uma função que faz basicamente a mesma coisa que o QGIS "dissolve" a função. Eu pensei que seria super fácil, mas aparentemente não. Então, pelo que entendi, o uso de fiona com bem torneado deve ser a melhor opção aqui. Eu comecei a mexer com arquivos vetoriais, então esse mundo é bem novo para mim e também para python.

Para este exemplo, estou trabalhando com um shapefile do condado fundado aqui http://tinyurl.com/odfbanu Então, aqui estão alguns trechos de código que reuni, mas não consigo encontrar uma maneira de fazê-los trabalhar juntos

Por enquanto, meu melhor método é o seguinte, com base em: https://sgillies.net/2009/01/27/a-more-perfect-union-continued.html . Funciona bem e recebo uma lista dos 52 estados como geometria Shapely. Por favor, não hesite em comentar se existe uma maneira mais direta de fazer esta parte.

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.

Portanto, a escrita não é realmente direta do que vi, quero apenas o mesmo arquivo de forma, apenas com o país se dissolvendo em estados, nem preciso muito da tabela de atributos, mas estou curioso para ver como você pode passar desde a fonte até o novo shapefile criado.

Encontrei muitos códigos para escrever com fiona, mas nunca consigo fazê-lo funcionar com meus dados. Exemplo de Como escrever geometrias bem torneadas em shapefiles? :

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

O problema aqui é como fazer o mesmo com uma lista de geometria e como recriar as mesmas propriedades que a origem.

A pergunta é sobre Fiona e Shapely e a outra resposta usando o GeoPandas também exige conhecer Pandas . Além disso, o GeoPandas usa o Fiona para ler / escrever arquivos de forma.

Não questiono aqui o utilitário do GeoPandas, mas você pode fazê-lo diretamente com o Fiona usando o módulo padrão itertools , especialmente com o comando groupby ("Em poucas palavras, o groupby pega um iterador e o divide em sub-iteradores com base nas alterações na "chave" do iterador principal. É claro que isso é feito sem a leitura de todo o iterador de origem na memória ", itertools.groupby ).

Shapefile original colorido pelo campo 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]})

Resultado

Eu recomendo o GeoPandas para lidar com uma grande variedade de recursos e executar operações em massa.

Ele estende os quadros de dados do Pandas e usa bem torneados sob o capô.

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

Como um adendo à resposta do @ gene , eu precisava me dissolver em mais de um campo para modificar seu código para lidar com vários campos. O código abaixo utiliza operator.itemgetterpara agrupar por vários campos:

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