Python - 匀称的大型 shapefile

Question

我正在读取一个 GeoJSON 文件，其中包含我制作的两个简单的多边形描述和来自 http://ryanmullins.org/blog/2015/8/18/land-area-vectors-for-geographic-combatant-commands

的六个复杂向量

我可以将自己的 4-8 点描述读入 Shapely Polygons。但是，上面网站的更复杂的描述给我以下错误：

from shapely.geometry import Polygon

jsonFile="path/to/file.json"

with open(jsonFile) as f:
    data=json.load(f)

    for feature in data['features']:
        #This is not how I'm saving the polygons, and is only for testing purposes:
        myPoly=Polygon(feature['geometry']['coordinates'])

错误信息：

File "/.../anaconda2/lib/python2.7/site-packages/shapely/geometry/polygon.py", line 229, in __init__
self._geom, self._ndim = geos_polygon_from_py(shell, holes)

File "/.../anaconda2/lib/python2.7/site-packages/shapely/geometry/polygon.py", line 508, in geos_polygon_from_py
geos_shell, ndim = geos_linearring_from_py(shell)

File "/.../anaconda2/lib/python2.7/site-packages/shapely/geometry/polygon.py", line 454, in geos_linearring_from_py
assert (n == 2 or n == 3)

AssertionError

它们被读取为 list，其中 USAFRICOM 的长度为 113。

有没有办法把这些很长的向量读成匀称的？我试过多边形、多点、asMultiPointIf。如果没有，您能否建议如何将此矢量描述简化为 Shapely 可以读取的内容？

Answer 1

好吧，它比一次性将所有坐标都扔到 Shapely 上要复杂一些。

根据关于多边形的GeoJSON spec and Shapely's documentation：多边形由多边形数组组成，多边形又由描述外部和内部区域/孔的线性环组成。

这是我为您的 GeoJSON 文件尝试的 MultiPolygon reader，在 QGIS 中打开输出显示正确，如果您有问题请告诉我。

import json
from shapely.geometry import mapping
from shapely.geometry import Polygon
from shapely.geometry import LinearRing
from shapely.geometry import MultiPolygon

jsonFile = "USCENTCOM.json"
polygons = []

with open(jsonFile) as f:
    data = json.load(f)
    for feature in data['features']:
        for multi_polygon in feature['geometry']['coordinates']:
            # collect coordinates (LinearRing coordinates) for the Polygon
            tmp_poly = [] 
            for polygon in multi_polygon:
                tmp_poly.append(polygon)
            if len(tmp_poly) > 1:
                # exterior LinearRing at [0], all following interior/"holes"
                polygons.append(Polygon(tmp_poly[0], tmp_poly[1:]))
            else:
                # there's just the exterior LinearRing
                polygons.append(Polygon(tmp_poly[0]))

# finally generate the MultiPolygon from our Polygons
mp = MultiPolygon(polygons)
# print GeoJSON string
print(json.dumps(mapping(mp)))

Python - 匀称的大型 shapefile

Python - Shapely large shapefiles

python-2.7

gis

shapely