2016-12-08 54 views
0

我在包含我做了兩個簡單的多邊形描述,並從http://ryanmullins.org/blog/2015/8/18/land-area-vectors-for-geographic-combatant-commandsPython的 - 身材勻稱大shape文件

6個複雜矢量我可以讀我自己的4-8點描述成身材勻稱多邊形A GeoJSON的文件中讀取。然而,從網站更復雜的描述上面給我以下錯誤:

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可以讀取的東西?

回答

0

那麼它會比在Shapely中一次性拋出所有座標複雜一點。

按照GeoJSON specShapely's documentation關於multipolygons:MultiPolygons由多邊形的陣列,和多邊形再次由該描繪的外部和內部區域/孔的LinearRing的。

這是我爲您的GeoJSON文件在MultiPolygon閱讀器上的嘗試,在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)))