GIS: Merging MultiPolygons in Python

Question

I have a geojson file of the counties in the British Isles, and I'm trying to replace the individual counties of London with a single merged county, then save the result as geojson. (The London counties can be identified because their TYPE_2 attribute is set to London Borough.)

I thought I could execute this task with the following:

from shapely.geometry import Polygon, MultiPolygon, asShape
from shapely.ops import unary_union
import json, geojson

j = json.load(open('british-isles.geojson'))

# find the london counties
indices = [idx for idx, i in enumerate(j['features']) if \
    i['properties']['TYPE_2'] == 'London Borough']

# transform each london county into a shapely polygon
polygons = [asShape(j['features'][i]['geometry']) for i in indices]

# get the metadata for the first county
properties = j['features'][indices[0]]['properties']
properties['NAME_2'] = 'London'

# get the union of the polygons
joined = unary_union(polygons)

# delete the merged counties
d = j
for i in indices:
  del d['features'][i]

# add the new polygon to the features
feature = geojson.Feature(geometry=joined, properties=properties)
d['features'].append(feature)

# save the geojson
with open('geojson-british-isles-merged-london.geojson', 'w') as out:
  json.dump(d, out)

However, this does not properly merge the London counties--it results instead in a fragmented series of polygons where the London counties used to be.

Do others know how I can accomplish this task in Python? Any suggestions would be very helpful!

Answer 1

Well there were two problems with the above. The first was purely an oversight: when deleting from d['features'], I needed to delete array members in reverse order (deleting index 0 then 1 is different from deleting 1 then 0).

More fundamentally, the above geojson was already lossy. Coordinate values had limited decimal places to shave bytes off the JSON file size. But this makes merging the geometries only approximate, and results in small gaps between the merged polygons:

So my resulting workflow was to obtain a high-res topojson file, convert it to geojson, merge the geometries using the code below, then limit decimal precision (if necessary).

from shapely.geometry import Polygon, MultiPolygon, asShape
from shapely.ops import unary_union, cascaded_union
from geojson import Feature
import json

j = json.load(open('GBR_adm2.json'))

# find the london counties
indices = [idx for idx, i in enumerate(j['features']) if \
    'London Borough' in i['properties']['TYPE_2']]

# transform each london county into a shapely polygon
polygons = [asShape(j['features'][i]['geometry']) for i in indices]

# get the metadata for the first county
properties = j['features'][indices[0]]['properties']
properties['NAME_2'] = 'London'

# get the union of the polygons
joined = unary_union(polygons)

# delete the merged counties
d = j
for i in reversed(sorted(indices)):
  del d['features'][i]

# add the new polygon to the features
feature = Feature(geometry=joined, properties=properties)
d['features'].append(feature)

# save the geojson
with open('british-isles-merged-london.geojson', 'w') as out:
  json.dump(d, out)

Result: