Python shapefile 和 matplotlib：绘制 shapefile 坐标的多边形集合
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/15968762/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
shapefile and matplotlib: plot polygon collection of shapefile coordinates
提问by hannesk
I'm trying to plot filled polygons of countries on the world map with matplotlib in python.
我正在尝试使用 python 中的 matplotlib 在世界地图上绘制国家的填充多边形。
I've got a shapefile with country boundary coordinates of every country. Now, I want to convert these coordinates (for each country) into a polygon with matplotlib. Without using Basemap. Unfortunately, the parts are crossing or overlapping. Is there a workarund, maybe using the distance from point to point.. or reordering them ?
我有一个包含每个国家/地区的国家边界坐标的 shapefile。现在，我想使用 matplotlib 将这些坐标（对于每个国家/地区）转换为多边形。不使用底图。不幸的是，这些部分交叉或重叠。是否有解决方法，也许使用点到点的距离......或重新排序它们？
回答by hannesk
Ha! I found out, how.. I completely neglected, the sf.shapes[i].parts information! Then it comes down to:
哈！我发现了，我完全忽略了 sf.shapes[i].parts 信息！然后归结为：
# -- 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')
Then it will look like this:
然后它看起来像这样：
回答by Ondro
Here is another piece of code I used to plot polygon shapefiles. It uses GDAL/OGR to read shapefile and plots correctly donut shape polygons:
这是我用来绘制多边形 shapefile 的另一段代码。它使用 GDAL/OGR 读取 shapefile 并正确绘制甜甜圈形状的多边形：
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()
回答by JMJR
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()
Note this assumes polygons, MultiPolygons can be handles in a similar manner with
请注意，这是假设多边形，MultiPolygons 可以以类似的方式处理
map(PolygonPatch, MultiPolygon(record['geometry']))
回答by Elendil
Regarding to @hannesk's answer, you should add the following imports: from numpy import arrayand import matplotliband replace the line cm = get_cmap('Dark2')by cm = matplotlib.cm.get_cmap('Dark2')
关于到@hannesk的回答，您应该添加以下的进口：from numpy import array和import matplotlib并更换线cm = get_cmap('Dark2')通过cm = matplotlib.cm.get_cmap('Dark2')
(I'm not so famous to add a comment to the noticed post.)
（我不是很出名，无法在注意到的帖子中添加评论。）
