Tuesday 16 April 2019

shapefile - Dissolving polygons based on multiple attributes with Python (shapely and fiona)?


I would like to dissolve a polygon shapefile based on two different fields, using Python and open source libraries.


I'm just able to dissolve the shapefile using one single field, with Gene's answer to Dissolving polygons based on attributes with Python (shapely, fiona)?


Here's the code he suggests, that should be the starting point for what I need:


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

Answer



Just change two key=lambda k: k['properties']['FIELD'] parts into key=lambda k: (k['properties']['FIELD_1'], k['properties']['FIELD_2']).


Shortly, you have to use ('field_1', 'field_2') tuple instead of 'field'.


from shapely.geometry import shape, mapping
from shapely.ops import unary_union
import fiona
import itertools

with fiona.open('source.shp') as input:

meta = input.meta
with fiona.open('target.shp', 'w', **meta) as output:
# First, sort data by field_1 then field_2
e = sorted(input, key=lambda k: (k['properties']['FIELD_1'], k['properties']['FIELD_2']) )
# group by field_1 then filed_2
for key, group in itertools.groupby(e, key=lambda k: (k['properties']['FIELD_1'], k['properties']['FIELD_2'])):
properties, geom = zip(*[(feature['properties'], shape(feature['geometry'])) for feature in group])
output.write({'geometry': mapping(unary_union(geom)), 'properties': properties[0]})

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...