2013-04-12 54 views
6

Tôi đang cố vẽ đồ thị đa giác chứa đầy các quốc gia trên bản đồ thế giới bằng matplotlib trong python.shapefile và matplotlib: tập hợp đa giác lô tọa độ shapefile

Tôi có một shapefile với các tọa độ ranh giới quốc gia của mỗi quốc gia. Bây giờ, tôi muốn chuyển đổi các tọa độ này (cho mỗi quốc gia) thành một đa giác với matplotlib. Không sử dụng Sơ đồ trang web. Thật không may, các bộ phận đang vượt qua hoặc chồng chéo. Có một workarund, có thể sử dụng khoảng cách từ điểm đến điểm .. hoặc sắp xếp lại chúng? enter image description here

+0

Bạn đã tạo hình ảnh này bằng cách nào? Tôi luôn sử dụng các đường dẫn cho điều này, như trong ví dụ này: http://matplotlib.org/examples/api/donut_demo.html –

Trả lời

11

Ha! Tôi phát hiện ra, làm thế nào .. Tôi hoàn toàn bị bỏ quên, các sf.shapes [i] .parts thông tin! Sau đó, nó đi xuống đến:

# -- import -- 
import shapefile 
import matplotlib.pyplot as plt 
import matplotlib.patches as patches 
from matplotlib.patches import Polygon 
from matplotlib.collections import PatchCollection 
# -- input -- 
sf = shapefile.Reader("./shapefiles/world_countries_boundary_file_world_2002") 
recs = sf.records() 
shapes = sf.shapes() 
Nshp = len(shapes) 
cns  = [] 
for nshp in xrange(Nshp): 
    cns.append(recs[nshp][1]) 
cns = array(cns) 
cm = get_cmap('Dark2') 
cccol = cm(1.*arange(Nshp)/Nshp) 
# -- plot -- 
fig  = plt.figure() 
ax  = fig.add_subplot(111) 
for nshp in xrange(Nshp): 
    ptchs = [] 
    pts  = array(shapes[nshp].points) 
    prt  = shapes[nshp].parts 
    par  = list(prt) + [pts.shape[0]] 
    for pij in xrange(len(prt)): 
    ptchs.append(Polygon(pts[par[pij]:par[pij+1]])) 
    ax.add_collection(PatchCollection(ptchs,facecolor=cccol[nshp,:],edgecolor='k', linewidths=.1)) 
ax.set_xlim(-180,+180) 
ax.set_ylim(-90,90) 
fig.savefig('test.png') 

Sau đó, nó sẽ giống như thế này: enter image description here

+0

Rất đẹp. Nếu tôi chạy kịch bản này tôi thấy rằng nó bỏ qua các vòng bên trong trong đa giác, không phải luôn luôn là một vấn đề, nhưng một cái gì đó để nhận thức được. –

+0

Nơi nào bạn xác định hàm 'get_cmap'? –

+0

Bạn có thể sửa đổi câu trả lời dựa trên câu trả lời của @Elendil về 'cm = matplotlib.cm.get_cmap ('Dark2')' –

3

Dưới đây là một đoạn mã tôi đã sử dụng để vẽ shapefile đa giác. Nó sử dụng GDAL/OGR để đọc shapefile và lô bánh rán đúng đa giác hình dạng:

from osgeo import ogr 
import numpy as np 
import matplotlib.path as mpath 
import matplotlib.patches as mpatches 
import matplotlib.pyplot as plt 

# Extract first layer of features from shapefile using OGR 
ds = ogr.Open('world_countries_boundary_file_world_2002.shp') 
nlay = ds.GetLayerCount() 
lyr = ds.GetLayer(0) 

# Get extent and calculate buffer size 
ext = lyr.GetExtent() 
xoff = (ext[1]-ext[0])/50 
yoff = (ext[3]-ext[2])/50 

# Prepare figure 
fig = plt.figure() 
ax = fig.add_subplot(111) 
ax.set_xlim(ext[0]-xoff,ext[1]+xoff) 
ax.set_ylim(ext[2]-yoff,ext[3]+yoff) 

paths = [] 
lyr.ResetReading() 

# Read all features in layer and store as paths 
for feat in lyr: 
    geom = feat.geometry() 
    codes = [] 
    all_x = [] 
    all_y = [] 
    for i in range(geom.GetGeometryCount()): 
     # Read ring geometry and create path 
     r = geom.GetGeometryRef(i) 
     x = [r.GetX(j) for j in range(r.GetPointCount())] 
     y = [r.GetY(j) for j in range(r.GetPointCount())] 
     # skip boundary between individual rings 
     codes += [mpath.Path.MOVETO] + \ 
        (len(x)-1)*[mpath.Path.LINETO] 
     all_x += x 
     all_y += y 
    path = mpath.Path(np.column_stack((all_x,all_y)), codes) 
    paths.append(path) 

# Add paths as patches to axes 
for path in paths: 
    patch = mpatches.PathPatch(path, \ 
      facecolor='blue', edgecolor='black') 
    ax.add_patch(patch) 

ax.set_aspect(1.0) 
plt.show() 
1
from fiona import collection 
import matplotlib.pyplot as plt 
from descartes import PolygonPatch 
from matplotlib.collections import PatchCollection 
from itertools import imap 
from matplotlib.cm import get_cmap 

cm = get_cmap('Dark2') 

figure, axes = plt.subplots(1) 

source_path = "./shapefiles/world_countries_boundary_file_world_2002" 

with collection(source_path, 'r') as source: 
    patches = imap(PolygonPatch, (record['geometry'] for record in source) 

axes.add_collection(PatchCollection (patches, cmap=cm, linewidths=0.1)) 

axes.set_xlim(-180,+180) 
axes.set_ylim(-90,90) 

plt.show() 

Lưu ý điều này giả định đa giác, MultiPolygons có thể xử lý một cách tương tự với

map(PolygonPatch, MultiPolygon(record['geometry'])) 
1

Về @hannesk ' Câu trả lời, bạn nên thêm các mục nhập sau: from numpy import arrayimport matplotlib và thay thế dòng cm = get_cmap('Dark2') theo cm = matplotlib.cm.get_cmap('Dark2')

(Tôi không nổi tiếng để thêm nhận xét vào bài đăng được chú ý.)

+0

Tôi hiểu tại sao bạn đã làm điều đó nhưng không nên làm điều đó. –

+0

Khi bạn có vẻ hiểu, điều này sẽ thích hợp hơn làm nhận xét. Để phê bình hoặc yêu cầu làm rõ từ tác giả, để lại nhận xét bên dưới bài đăng của họ - bạn luôn có thể nhận xét về bài đăng của riêng bạn và sau khi bạn có đủ [danh tiếng] (http://stackoverflow.com/help/whats-reputation), bạn sẽ có thể [nhận xét về bài đăng bất kỳ] (http://stackoverflow.com/help/privileges/comment). – Makyen

+0

Lần sau tôi để lỗi trong bài đăng mà không nói bất cứ điều gì anh chàng khác sẽ sửa nó trong tương lai nếu họ tìm thấy nó! Tôi chỉ muốn giúp đỡ như tôi có thể nhưng tôi không thể làm điều đó như Stack mong đợi bằng cách viết bình luận ... Và nhờ lá phiếu của bạn tôi sẽ không bao giờ có quyền đăng bình luận :-) – Elendil

Các vấn đề liên quan