¿Disolver polígonos basándose en atributos con Python (shapely, fiona)?

Question

He estado intentando crear una función que haga básicamente lo mismo que la función "disolver" de QGIS. Pensé que sería super fácil pero bueno al parecer no. Así que de lo que he reunido alrededor, el uso de fiona con shapely debe ser la mejor opción aquí. Acabo de empezar a trastear con archivos vectoriales así que este mundo es bastante nuevo para mi y python también.

Para este ejemplo, estoy trabajando con un shapefile del condado fundado aquí http://tinyurl.com/odfbanu Así que aquí están algunas piezas de código que reuní, pero no puedo encontrar una manera de hacer que funcionen juntos

Por ahora mi mejor método es el siguiente basado en : https://sgillies.net/2009/01/27/a-more-perfect-union-continued.html . Funciona bien y obtengo una lista de los 52 estados como geometría Shapely. Por favor, no dude en comentar si hay una manera más directa de hacer 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.

Así que la escritura no es realmente sencillo de lo que he visto, en realidad sólo quiero el mismo shapefile sólo con el país disolver en los estados, ni siquiera necesito mucho de la tabla de atributos, pero tengo curiosidad por ver cómo se puede pasar de la fuente a la nueva shapefile creado.

He encontrado muchos trozos de código para escribir con fiona pero nunca soy capaz de hacerlo funcionar con mis datos. Ejemplo de ¿Cómo escribir geometrías Shapely en 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},
    })

El problema aquí es cómo hacer lo mismo con una lista de geometría y cómo recrear las mismas propiedades que la fuente.

Answer 1

La pregunta es sobre Fiona y Shapely y la otra respuesta utilizando GeoPandas requiere conocer también Pandas . Además, GeoPandas utiliza Fiona para leer y escribir archivos shape.

No cuestiono aquí la utilidad de GeoPandas, pero puedes hacerlo directamente con Fiona utilizando el módulo estándar itertools especialmente con el comando agruparpor ("En pocas palabras, groupby toma un iterador y lo divide en sub-iteradores basados en cambios en la "clave" del iterador principal. Esto se hace, por supuesto, sin leer todo el iterador fuente en memoria", itertools.groupby ).

Shapefile original coloreado por el 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

Answer 2

Recomiendo encarecidamente GeoPandas para tratar con grandes surtidos de características y realizar operaciones masivas.

Extiende Pandas dataframes, y utiliza shapely bajo el 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")

Answer 3

Como adición a @gene's responder En mi caso, necesitaba disolver por más de un campo, así que modifiqué su código para manejar múltiples campos. El siguiente código utiliza operator.itemgetter para agrupar por varios campos:

# Modified from https://gis.stackexchange.com/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"))