14

Is there a way to merge two overlapping GEOJSON polygons in python, returning a single merged GEOJSON object?

traggatmot
  • 1,423
  • 5
  • 26
  • 51
  • Pure python, or are you looking for a library, or what? – Erica Dec 17 '15 at 02:05
  • Prefer pure python, but open to anything. I'm thinking about trying pygeoproseccing package, but I can't get it to install, because GDAL won't install....so, I am open to suggestions while I grapple with that. – traggatmot Dec 17 '15 at 02:10
  • [shapely](https://pypi.python.org/pypi/Shapely) is great for geospatial operations like that. I have no experience with GeoJSON though. – Matt Hall Dec 17 '15 at 03:00
  • I'm reading through the shapely manual now and it doesn't seem like it can easily digest GeoJSON objects. – traggatmot Dec 17 '15 at 03:34

3 Answers3

11

This is how I was able to do it using the packages/modules json, geojson, shapely, pyproj, and partial from functools:

import json
import geojson
from functools import partial
import pyproj
import shapely.geometry
import shapely.ops

# reading into two geojson objects, in a GCS (WGS84)
with open('file1.json') as geojson1:
    poly1_geojson = json.load(geojson1)

with open('file2.json') as geojson2:
    poly2_geojson = json.load(geojson2)


# pulling out the polygons
poly1 = shapely.geometry.asShape(poly1_geojson['features'][2]['geometry'])
poly2 = shapely.geometry.asShape(poly2_geojson['features'][2]['geometry'])

# checking to make sure they registered as polygons
print poly1.geom_type
print poly2.geom_type

# merging the polygons - they are feature collections, containing a point, a polyline, and a polygon - I extract the polygon
# for my purposes, they overlap, so merging produces a single polygon rather than a list of polygons
mergedPolygon = poly1.union(poly2)

# using geojson module to convert from WKT back into GeoJSON format
geojson_out = geojson.Feature(geometry=mergedPolygon, properties={})

# outputting the updated geojson file - for mapping/storage in its GCS format
with open('Merged_Polygon.json', 'w') as outfile:
    json.dump(geojson_out.geometry, outfile, indent=3, encoding="utf-8")
outfile.close()

# reprojecting the merged polygon to determine the correct area
# it is a polygon covering much of the US, and dervied form USGS data, so using Albers Equal Area
project = partial(
    pyproj.transform,
    pyproj.Proj(init='epsg:4326'),
    pyproj.Proj(init='epsg:5070'))

mergedPolygon_proj = shapely.ops.transform(project,mergedPolygon)
traggatmot
  • 1,423
  • 5
  • 26
  • 51
11

This example from here seems to be a lot more concise:

from shapely.geometry import Polygon
from shapely.ops import cascaded_union

polygon1 = Polygon([(0, 0), (5, 3), (5, 0)])
polygon2 = Polygon([(0, 0), (3, 10), (3, 0)])

polygons = [polygon1, polygon2]

u = cascaded_union(polygons)
Luke Madhanga
  • 6,871
  • 2
  • 43
  • 47
0

The dissolve() function within GeoPandas makes this very simple. For example, I had a GeoDataFrame containing US counties (along with their polygons) and their corresponding Catholic diocese. I wanted to create polygons showing the outline of each diocese for each county. To do so, I used the following code:

diocese_boundaries = df_dioceses.dissolve(by = 'Diocese')

This returned a new DataFrame containing one row for each diocese. Each row had a new geometry column that contained the outline of each diocese.

KBurchfiel
  • 635
  • 6
  • 15