Python - Shapely большие шейп-файлы

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

Они читаются как список, с USAFRICOM длиной 113.

Есть ли способ считать эти очень длинные векторы красивыми? Я пробовал Polygon, MultiPoint, asMultiPointIf. Если нет, то не могли бы вы предложить, как упростить это векторное описание во что-то, что может быть прочитано Shapely?

Answer 1

Ну, это немного сложнее, чем просто бросить все координаты в Shapely за один раз.

Согласно спецификации GeoJSON и документации Shapely относительно мультиполигонов: мультиполигоны состоят из массивов полигонов, а полигоны снова состоят из линейных колец, которые изображают внешние и внутренние области / отверстия.

Вот моя попытка чтения MultiPolygon для ваших файлов GeoJSON, открытие вывода обратно в 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)))