代码之家 › 专栏 › 技术社区 › duhaime

GIS:在Python中合并多多边形

shapely geojson gis json python

duhaime · 技术社区 · 6 年前

我有一个 geojson file 我试图用一个合并的郡替换伦敦的各个郡,然后将结果保存为geojson。(伦敦各县之所以能够被识别,是因为它们 TYPE_2 属性设置为 London Borough .)

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)

然而,这并不能很好地合并伦敦县——而是导致了一系列支离破碎的多边形,而伦敦县过去是这样的。

其他人知道我如何用Python完成这个任务吗?任何建议都会很有帮助!

1 回复 | 直到 6 年前

duhaime 6 年前

上面有两个问题。第一个纯粹是疏忽:从 d['features'] ,我需要按相反的顺序删除数组成员(先删除索引0,然后删除1与先删除1,然后删除0不同)。

更重要的是,上面的geojson已经是有损的了。坐标值的小数位数有限,可以减少JSON文件大小的字节数。但这使得合并几何体只是近似的,并导致合并的多边形之间的小间隙:

所以我的工作流程是获得高分辨率 topojson file ,将其转换为geojson,使用下面的代码合并几何体, 然后限制小数精度(如有必要)。

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)

结果: