在 Python 中合并两个 GEOJSON 多边形
Merging two GEOJSON polygons in Python
有没有办法在 python 中合并两个重叠的 GEOJSON 多边形,返回一个合并的 GEOJSON 对象?
这就是我使用 packages/modules json、geojson、shapely、pyproj 和来自 functools 的 partial 来做到这一点的方法:
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)
here 中的这个示例似乎更加简洁:
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)
GeoPandas 中的 dissolve() 函数使这变得非常简单。例如,我有一个包含美国县(及其多边形)及其相应的天主教教区的 GeoDataFrame。我想创建多边形来显示每个县每个教区的轮廓。为此,我使用了以下代码:
diocese_boundaries = df_dioceses.dissolve(by = 'Diocese')
这返回了一个新的 DataFrame,其中包含每个教区的一行。每行都有一个新的几何列,其中包含每个教区的轮廓。
有没有办法在 python 中合并两个重叠的 GEOJSON 多边形,返回一个合并的 GEOJSON 对象?
这就是我使用 packages/modules json、geojson、shapely、pyproj 和来自 functools 的 partial 来做到这一点的方法:
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)
here 中的这个示例似乎更加简洁:
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)
GeoPandas 中的 dissolve() 函数使这变得非常简单。例如,我有一个包含美国县(及其多边形)及其相应的天主教教区的 GeoDataFrame。我想创建多边形来显示每个县每个教区的轮廓。为此,我使用了以下代码:
diocese_boundaries = df_dioceses.dissolve(by = 'Diocese')
这返回了一个新的 DataFrame,其中包含每个教区的一行。每行都有一个新的几何列,其中包含每个教区的轮廓。