python - GIS:在 Python 中合并多面体
问题描述
我有一个不列颠群岛县的geojson文件,我试图用一个合并的县替换伦敦的各个县,然后将结果保存为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 中完成这项任务吗?任何建议都会非常有帮助!
解决方案
那么上面有两个问题。第一个纯粹是一个疏忽:从 删除时d['features']
,我需要以相反的顺序删除数组成员(删除索引 0 然后 1 与删除 1 然后 0 不同)。
更根本的是,上面的 geojson 已经是有损的了。坐标值的小数位数有限,以减少 JSON 文件大小的字节数。但这使得合并几何图形只是近似的,并导致合并多边形之间的小间隙:
所以我得到的工作流程是获取一个高分辨率的topojson 文件,将其转换为 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)
推荐阅读
- c# - 在 C# Visual Studio 中,一旦应用程序使用了设定的内存量,是否可以拍摄内存快照?
- azure-devops - 无法在 Azure DevOps 中获得单个平台的徽章
- security - 攻击者可以从 AES-GCM 密文中学到什么?
- python - Subset pandas data frame based on tuple
- python - “set”类型的对象不是 JSON 可序列化的
- sql - 查询以根据 UI 中的搜索条件从多个表中选择记录
- search - 对大于或等于给定键的元素进行二分查找
- python - Function to return the index at which a specific letter occurs always returns nothing regardless of input in Python?
- git - 如何忽略“git存档”中的文件/目录,只创建子目录的存档?
- android - 如何为字符串中的文本添加下划线?