使用 LineCollection 绘制 shapefile 显示所有边缘，但部分填充它们
Plotting shapefile using LineCollection shows all edges, but partially fills them
在过去的几天里，我一直在尝试将气象站数据插值到我所在国家/地区的地图中。我这样做如下：
- 加载数据，我使用插值法创建了一个网格
- 基于这个网格我画了一个
contour 和 contourf 图像
- 然后我在地图顶部绘制德国、比利时和法国的 shapefile，以覆盖不相关的
contour/contourf 元素。为此，我使用了 this tutorial。
- 最后，我还使用海洋 shapefile (IHO Sea Areas; www.marineregions.org/downloads.php#iho) 绘制以覆盖北海。使用 QGIS，我编辑了这个海洋 shapefile 并删除了除北海以外的所有内容 - 给定时间限制:)
你会说一切顺利 -- 但由于某种原因，这个国家的部分地区以及岛屿被解释为水域。我想这是因为这些是不同的部分，都没有连接到主要土地（由于 waters/rivers）。
奇怪的是，画了边，但没有填充。
我已经尝试并搜索了很多，但不知道如何解决这个问题。我猜它在 LineCollection 中的某个地方，因为在 QGIS 中，shapefile 是正确的（即它在单击岛屿时识别 no shape 等，这是正确的，因为它应该只识别单击海时的形状）。
我真诚地希望你能帮我指出我错在哪里以及我该如何解决这个问题！
非常感谢！
这是我得到的地图：
https://i.imgur.com/GHISN7n.png
我的代码如下（你可能会看到我对这种编程很陌生，我昨天开始:)):
import numpy as np
import matplotlib
matplotlib.use('Agg')
from scipy.interpolate import griddata
from mpl_toolkits.basemap import Basemap, maskoceans
import matplotlib.pyplot as plt
from numpy.random import seed
import shapefile as shp
from matplotlib.collections import LineCollection
from matplotlib import cm
# Set figure size
plt.figure(figsize=(15,15), dpi=80)
# Define map bounds
xMin, xMax = 2.5, 8.0
yMin, yMax = 50.6, 53.8
# Create map
m = Basemap(projection='merc',llcrnrlon=xMin,llcrnrlat=yMin,urcrnrlon=xMax,urcrnrlat=yMax,resolution='h')
m.drawmapboundary(fill_color='#d4dadc',linewidth=0.25)
# m.drawcoastlines(linewidth=0.5,color='#333333')
# Load data
y = [54.325666666667,52.36,53.269444444444,55.399166666667,54.116666666667,53.614444444444,53.491666666667,53.824130555556,52.918055555556,54.03694,52.139722,52.926865008825,54.853888888889,52.317222,53.240026656696,52.642696895243,53.391265948394,52.505333893732,52.098821802977,52.896643913235,52.457270486008,53.223000488316,52.701902388132,52.0548617826,53.411581103636,52.434561756559,52.749056395511,53.123676213651,52.067534268959,53.194409573306,52.27314817052,51.441334059998,51.224757511326,51.990941918858,51.447744494043,51.960667359998,51.969031121385,51.564889021961,51.857593837453,51.449772459909,51.658528382201,51.196699902606,50.905256257898,51.497306260089,yMin,yMin,yMax,yMax]
x = [2.93575,3.3416666666667,3.6277777777778,3.8102777777778,4.0122222222222,4.9602777777778,5.9416666666667,2.9452777777778,4.1502777777778,6.04167,4.436389,4.7811453228565,4.6961111111111,4.789722,4.9207907082729,4.9787572406902,5.3458010937365,4.6029300588208,5.1797058644882,5.383478899702,5.5196324030324,5.7515738887123,5.8874461671401,5.8723225499118,6.1990994508938,6.2589770334531,6.5729701105864,6.5848470019087,6.6567253619722,7.1493220605216,6.8908745111116,3.5958241584686,3.8609657214986,4.121849767852,4.342014,4.4469005114756,4.9259216999194,4.9352386335384,5.1453989235756,5.3770039280214,5.7065946674719,5.7625447234516,5.7617834850481,6.1961067840608,xMin,xMax,xMin,xMax]
z = [4.8,5.2,5.8,5.4,5,5.3,5.4,4.6,5.8,6.3,4.8,5.4,5.3,4.6,5.4,4.4,4.1,5.5,4.5,4.2,3.9,3.7,4.2,3.2,4,3.8,2.7,2.3,3.4,2.5,3.7,5.2,2.9,5.1,3.8,4.4,4.2,3.9,3.8,3.2,2.6,2.8,2.4,3.1]
avg = np.average(z)
z.extend([avg,avg,avg,avg])
x,y = m(x,y)
# target grid to interpolate to
xis = np.arange(min(x),max(x),2000)
yis = np.arange(min(y),max(y),2000)
xi,yi = np.meshgrid(xis,yis)
# interpolate
zi = griddata((x,y),z,(xi,yi),method='cubic')
# Decide on proper values for colour bar (todo)
vrange = max(z)-min(z)
mult = 2
vmin = min(z)-(mult*vrange)
vmax = max(z)+(mult*vrange)
# Draw contours
cs = m.contour(xi, yi, zi, 5, linewidths=0.25, colors='k')
cs = m.contourf(xi, yi, zi, 5,vmax=vmax,vmin=vmin,cmap=plt.get_cmap('jet'))
# Plot seas from shapefile
sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/northsea')
shapes = sf.shapes()
print shapes[0].parts
records = sf.records()
ax = plt.subplot(111)
for record, shape in zip(records,shapes):
lons,lats = zip(*shape.points)
data = np.array(m(lons, lats)).T
print len(shape.parts)
if len(shape.parts) == 1:
segs = [data,]
else:
segs = []
for i in range(1,len(shape.parts)):
index = shape.parts[i-1]
index2 = shape.parts[i]
segs.append(data[index:index2])
segs.append(data[index2:])
lines = LineCollection(segs,antialiaseds=(1,),zorder=3)
lines.set_facecolors('#abc0d3')
lines.set_edgecolors('red')
lines.set_linewidth(0.5)
ax.add_collection(lines)
# Plot France from shapefile
sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/FRA_adm0')
shapes = sf.shapes()
records = sf.records()
ax = plt.subplot(111)
for record, shape in zip(records,shapes):
lons,lats = zip(*shape.points)
data = np.array(m(lons, lats)).T
if len(shape.parts) == 1:
segs = [data,]
else:
segs = []
for i in range(1,len(shape.parts)):
index = shape.parts[i-1]
index2 = shape.parts[i]
segs.append(data[index:index2])
segs.append(data[index2:])
lines = LineCollection(segs,antialiaseds=(1,))
lines.set_facecolors('#fafaf8')
lines.set_edgecolors('#333333')
lines.set_linewidth(0.5)
ax.add_collection(lines)
# Plot Belgium from shapefile
sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/BEL_adm0')
shapes = sf.shapes()
records = sf.records()
ax = plt.subplot(111)
for record, shape in zip(records,shapes):
lons,lats = zip(*shape.points)
data = np.array(m(lons, lats)).T
if len(shape.parts) == 1:
segs = [data,]
else:
segs = []
for i in range(1,len(shape.parts)):
index = shape.parts[i-1]
index2 = shape.parts[i]
segs.append(data[index:index2])
segs.append(data[index2:])
lines = LineCollection(segs,antialiaseds=(1,))
lines.set_facecolors('#fafaf8')
lines.set_edgecolors('#333333')
lines.set_linewidth(0.5)
ax.add_collection(lines)
# Plot Germany from shapefile
sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/DEU_adm0')
shapes = sf.shapes()
records = sf.records()
ax = plt.subplot(111)
for record, shape in zip(records,shapes):
lons,lats = zip(*shape.points)
data = np.array(m(lons, lats)).T
if len(shape.parts) == 1:
segs = [data,]
else:
segs = []
for i in range(1,len(shape.parts)):
index = shape.parts[i-1]
index2 = shape.parts[i]
segs.append(data[index:index2])
segs.append(data[index2:])
lines = LineCollection(segs,antialiaseds=(1,))
lines.set_facecolors('#fafaf8')
lines.set_edgecolors('#333333')
lines.set_linewidth(0.5)
ax.add_collection(lines)
# Finish plot
plt.axis('off')
plt.savefig("test2.png",bbox_inches='tight',pad_inches=0);
你还没有绘制荷兰。
# Create map
...
# Draw contours
...
# Plot seas from shapefile
sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/northsea')
...
# Plot France from shapefile
sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/FRA_adm0')
...
# Plot Belgium from shapefile
sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/BEL_adm0')
...
# Plot Germany from shapefile
sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/DEU_adm0')
...
# Finish plot
所以荷兰的部分地区填充了您绘制的数据，海岸线在那里（通过绘制海洋？），但该国其他地区不存在，因为您没有绘制它。
我会避开这种从每个 shapefile 单独绘制的方法，并尝试只为您的海岸线等使用底图。如果您必须从 shapefile 中执行此操作，我鼓励您定义一个函数来执行类似于：
def drawContentsOfShapeFile(mymap, path_to_shapefile):
<whatever>
#and use it:
m = Basemap(...)
myshapefiles = []
myshapefiles.append("/DIR_TO_SHP/shapefiles/FRA_adm0")
myshapefiles.append("/DIR_TO_SHP/shapefiles/norway")
for sf in myshapefiles: drawContentOfShapefile(m, sf)
你的问题是 LineCollection 没有按照你的想法去做。你想要的是一个带有 'holes punched in' 的外部填充形状，它具有 LineCollection 中其他线条的形状（即北海中的岛屿）。 LineCollection，然而，填充列表中的每个线段，所以填充的形状只是相互重叠。
Inspired by this post, I wrote an answer that seems to solve your problem using Patches. However, I'm not entirely sure how robust the solution is: according to the linked (unanswered) post, the vertices of the outer shape should be ordered clockwise while the vertices of the 'punched in' holes should be ordered counter-clock wise (I checked that too and it appears to be correct; this matplotlib example 显示相同的行为）。由于 shapefile 是二进制的，因此很难验证顶点的顺序，但结果似乎是正确的。在下面的示例中，我假设在每个 shapefile 中最长的线段是国家（或北海）的轮廓，而较短的线段是岛屿或类似的。因此，我首先按长度对每个 shapefile 的段进行排序，然后创建一个 Path 和一个 PathPatch。我自由地为每个 shapefile 使用不同的颜色，以确保一切正常。
import numpy as np
import matplotlib
matplotlib.use('Agg')
from scipy.interpolate import griddata
from mpl_toolkits.basemap import Basemap, maskoceans
import matplotlib.pyplot as plt
from numpy.random import seed
import shapefile as shp
from matplotlib.collections import LineCollection
from matplotlib.patches import Path, PathPatch
from matplotlib import cm
# Set figure size
fig, ax = plt.subplots(figsize=(15,15), dpi = 80)
# Define map bounds
xMin, xMax = 2.5, 8.0
yMin, yMax = 50.6, 53.8
shapefiles = [
'shapefiles/BEL_adm0',
'shapefiles/FRA_adm0',
'shapefiles/DEU_adm0',
'shapefiles/northsea',
]
colors = ['red', 'green', 'yellow', 'blue']
y = [54.325666666667,52.36,53.269444444444,55.399166666667,54.116666666667,53.614444444444,53.491666666667,53.824130555556,52.918055555556,54.03694,52.139722,52.926865008825,54.853888888889,52.317222,53.240026656696,52.642696895243,53.391265948394,52.505333893732,52.098821802977,52.896643913235,52.457270486008,53.223000488316,52.701902388132,52.0548617826,53.411581103636,52.434561756559,52.749056395511,53.123676213651,52.067534268959,53.194409573306,52.27314817052,51.441334059998,51.224757511326,51.990941918858,51.447744494043,51.960667359998,51.969031121385,51.564889021961,51.857593837453,51.449772459909,51.658528382201,51.196699902606,50.905256257898,51.497306260089,yMin,yMin,yMax,yMax]
x = [2.93575,3.3416666666667,3.6277777777778,3.8102777777778,4.0122222222222,4.9602777777778,5.9416666666667,2.9452777777778,4.1502777777778,6.04167,4.436389,4.7811453228565,4.6961111111111,4.789722,4.9207907082729,4.9787572406902,5.3458010937365,4.6029300588208,5.1797058644882,5.383478899702,5.5196324030324,5.7515738887123,5.8874461671401,5.8723225499118,6.1990994508938,6.2589770334531,6.5729701105864,6.5848470019087,6.6567253619722,7.1493220605216,6.8908745111116,3.5958241584686,3.8609657214986,4.121849767852,4.342014,4.4469005114756,4.9259216999194,4.9352386335384,5.1453989235756,5.3770039280214,5.7065946674719,5.7625447234516,5.7617834850481,6.1961067840608,xMin,xMax,xMin,xMax]
z = [4.8,5.2,5.8,5.4,5,5.3,5.4,4.6,5.8,6.3,4.8,5.4,5.3,4.6,5.4,4.4,4.1,5.5,4.5,4.2,3.9,3.7,4.2,3.2,4,3.8,2.7,2.3,3.4,2.5,3.7,5.2,2.9,5.1,3.8,4.4,4.2,3.9,3.8,3.2,2.6,2.8,2.4,3.1]
avg = np.average(z)
z.extend([avg,avg,avg,avg])
# Create map
m = Basemap(
ax = ax,
projection='merc',
llcrnrlon=xMin,
llcrnrlat=yMin,
urcrnrlon=xMax,
urcrnrlat=yMax,
resolution='h'
)
x,y = m(x,y)
m.drawmapboundary(fill_color='#d4dadc',linewidth=0.25)
# target grid to interpolate to
xis = np.arange(min(x),max(x),2000)
yis = np.arange(min(y),max(y),2000)
xi,yi = np.meshgrid(xis,yis)
# interpolate
zi = griddata((x,y),z,(xi,yi),method='cubic')
# Decide on proper values for colour bar (todo)
vrange = max(z)-min(z)
mult = 2
vmin = min(z)-(mult*vrange)
vmax = max(z)+(mult*vrange)
# Draw contours
cs = m.contour(xi, yi, zi, 5, linewidths=0.25, colors='k')
cs = m.contourf(xi, yi, zi, 5,vmax=vmax,vmin=vmin,cmap=plt.get_cmap('jet'))
for sf_name,color in zip(shapefiles, colors):
print(sf_name)
# Load data
#drawing shapes:
sf = shp.Reader(sf_name)
shapes = sf.shapes()
##print shapes[0].parts
records = sf.records()
##ax = plt.subplot(111)
for record, shape in zip(records,shapes):
lons,lats = zip(*shape.points)
data = np.array(m(lons, lats)).T
if len(shape.parts) == 1:
segs = [data,]
else:
segs = []
for i in range(1,len(shape.parts)):
index = shape.parts[i-1]
index2 = shape.parts[i]
seg = data[index:index2]
segs.append(seg)
segs.append(data[index2:])
##assuming that the longest segment is the enclosing
##line and ordering the segments by length:
lens=np.array([len(s) for s in segs])
order = lens.argsort()[::-1]
segs = [segs[i] for i in order]
##producing the outlines:
lines = LineCollection(segs,antialiaseds=(1,),zorder=4)
##note: leaving the facecolors out:
##lines.set_facecolors('#abc0d3')
lines.set_edgecolors('red')
lines.set_linewidth(0.5)
ax.add_collection(lines)
##producing a path from the line segments:
segs_lin = [v for s in segs for v in s]
codes = [
[Path.MOVETO]+
[Path.LINETO for p in s[1:]]
for s in segs]
codes_lin = [c for s in codes for c in s]
path = Path(segs_lin, codes_lin)
##patch = PathPatch(path, facecolor="#abc0d3", lw=0, zorder = 3)
patch = PathPatch(path, facecolor=color, lw=0, zorder = 3)
ax.add_patch(patch)
plt.axis('off')
fig.savefig("shapefiles.png",bbox_inches='tight',pad_inches=0)
结果如下所示：
希望对您有所帮助。
在过去的几天里，我一直在尝试将气象站数据插值到我所在国家/地区的地图中。我这样做如下：
- 加载数据，我使用插值法创建了一个网格
- 基于这个网格我画了一个
contour和contourf图像 - 然后我在地图顶部绘制德国、比利时和法国的 shapefile，以覆盖不相关的
contour/contourf元素。为此，我使用了 this tutorial。 - 最后，我还使用海洋 shapefile (IHO Sea Areas; www.marineregions.org/downloads.php#iho) 绘制以覆盖北海。使用 QGIS，我编辑了这个海洋 shapefile 并删除了除北海以外的所有内容 - 给定时间限制:)
你会说一切顺利 -- 但由于某种原因，这个国家的部分地区以及岛屿被解释为水域。我想这是因为这些是不同的部分，都没有连接到主要土地（由于 waters/rivers）。
奇怪的是，画了边，但没有填充。
我已经尝试并搜索了很多，但不知道如何解决这个问题。我猜它在 LineCollection 中的某个地方，因为在 QGIS 中，shapefile 是正确的（即它在单击岛屿时识别 no shape 等，这是正确的，因为它应该只识别单击海时的形状）。
我真诚地希望你能帮我指出我错在哪里以及我该如何解决这个问题！
非常感谢！
这是我得到的地图：
https://i.imgur.com/GHISN7n.png
我的代码如下（你可能会看到我对这种编程很陌生，我昨天开始:)):
import numpy as np
import matplotlib
matplotlib.use('Agg')
from scipy.interpolate import griddata
from mpl_toolkits.basemap import Basemap, maskoceans
import matplotlib.pyplot as plt
from numpy.random import seed
import shapefile as shp
from matplotlib.collections import LineCollection
from matplotlib import cm
# Set figure size
plt.figure(figsize=(15,15), dpi=80)
# Define map bounds
xMin, xMax = 2.5, 8.0
yMin, yMax = 50.6, 53.8
# Create map
m = Basemap(projection='merc',llcrnrlon=xMin,llcrnrlat=yMin,urcrnrlon=xMax,urcrnrlat=yMax,resolution='h')
m.drawmapboundary(fill_color='#d4dadc',linewidth=0.25)
# m.drawcoastlines(linewidth=0.5,color='#333333')
# Load data
y = [54.325666666667,52.36,53.269444444444,55.399166666667,54.116666666667,53.614444444444,53.491666666667,53.824130555556,52.918055555556,54.03694,52.139722,52.926865008825,54.853888888889,52.317222,53.240026656696,52.642696895243,53.391265948394,52.505333893732,52.098821802977,52.896643913235,52.457270486008,53.223000488316,52.701902388132,52.0548617826,53.411581103636,52.434561756559,52.749056395511,53.123676213651,52.067534268959,53.194409573306,52.27314817052,51.441334059998,51.224757511326,51.990941918858,51.447744494043,51.960667359998,51.969031121385,51.564889021961,51.857593837453,51.449772459909,51.658528382201,51.196699902606,50.905256257898,51.497306260089,yMin,yMin,yMax,yMax]
x = [2.93575,3.3416666666667,3.6277777777778,3.8102777777778,4.0122222222222,4.9602777777778,5.9416666666667,2.9452777777778,4.1502777777778,6.04167,4.436389,4.7811453228565,4.6961111111111,4.789722,4.9207907082729,4.9787572406902,5.3458010937365,4.6029300588208,5.1797058644882,5.383478899702,5.5196324030324,5.7515738887123,5.8874461671401,5.8723225499118,6.1990994508938,6.2589770334531,6.5729701105864,6.5848470019087,6.6567253619722,7.1493220605216,6.8908745111116,3.5958241584686,3.8609657214986,4.121849767852,4.342014,4.4469005114756,4.9259216999194,4.9352386335384,5.1453989235756,5.3770039280214,5.7065946674719,5.7625447234516,5.7617834850481,6.1961067840608,xMin,xMax,xMin,xMax]
z = [4.8,5.2,5.8,5.4,5,5.3,5.4,4.6,5.8,6.3,4.8,5.4,5.3,4.6,5.4,4.4,4.1,5.5,4.5,4.2,3.9,3.7,4.2,3.2,4,3.8,2.7,2.3,3.4,2.5,3.7,5.2,2.9,5.1,3.8,4.4,4.2,3.9,3.8,3.2,2.6,2.8,2.4,3.1]
avg = np.average(z)
z.extend([avg,avg,avg,avg])
x,y = m(x,y)
# target grid to interpolate to
xis = np.arange(min(x),max(x),2000)
yis = np.arange(min(y),max(y),2000)
xi,yi = np.meshgrid(xis,yis)
# interpolate
zi = griddata((x,y),z,(xi,yi),method='cubic')
# Decide on proper values for colour bar (todo)
vrange = max(z)-min(z)
mult = 2
vmin = min(z)-(mult*vrange)
vmax = max(z)+(mult*vrange)
# Draw contours
cs = m.contour(xi, yi, zi, 5, linewidths=0.25, colors='k')
cs = m.contourf(xi, yi, zi, 5,vmax=vmax,vmin=vmin,cmap=plt.get_cmap('jet'))
# Plot seas from shapefile
sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/northsea')
shapes = sf.shapes()
print shapes[0].parts
records = sf.records()
ax = plt.subplot(111)
for record, shape in zip(records,shapes):
lons,lats = zip(*shape.points)
data = np.array(m(lons, lats)).T
print len(shape.parts)
if len(shape.parts) == 1:
segs = [data,]
else:
segs = []
for i in range(1,len(shape.parts)):
index = shape.parts[i-1]
index2 = shape.parts[i]
segs.append(data[index:index2])
segs.append(data[index2:])
lines = LineCollection(segs,antialiaseds=(1,),zorder=3)
lines.set_facecolors('#abc0d3')
lines.set_edgecolors('red')
lines.set_linewidth(0.5)
ax.add_collection(lines)
# Plot France from shapefile
sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/FRA_adm0')
shapes = sf.shapes()
records = sf.records()
ax = plt.subplot(111)
for record, shape in zip(records,shapes):
lons,lats = zip(*shape.points)
data = np.array(m(lons, lats)).T
if len(shape.parts) == 1:
segs = [data,]
else:
segs = []
for i in range(1,len(shape.parts)):
index = shape.parts[i-1]
index2 = shape.parts[i]
segs.append(data[index:index2])
segs.append(data[index2:])
lines = LineCollection(segs,antialiaseds=(1,))
lines.set_facecolors('#fafaf8')
lines.set_edgecolors('#333333')
lines.set_linewidth(0.5)
ax.add_collection(lines)
# Plot Belgium from shapefile
sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/BEL_adm0')
shapes = sf.shapes()
records = sf.records()
ax = plt.subplot(111)
for record, shape in zip(records,shapes):
lons,lats = zip(*shape.points)
data = np.array(m(lons, lats)).T
if len(shape.parts) == 1:
segs = [data,]
else:
segs = []
for i in range(1,len(shape.parts)):
index = shape.parts[i-1]
index2 = shape.parts[i]
segs.append(data[index:index2])
segs.append(data[index2:])
lines = LineCollection(segs,antialiaseds=(1,))
lines.set_facecolors('#fafaf8')
lines.set_edgecolors('#333333')
lines.set_linewidth(0.5)
ax.add_collection(lines)
# Plot Germany from shapefile
sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/DEU_adm0')
shapes = sf.shapes()
records = sf.records()
ax = plt.subplot(111)
for record, shape in zip(records,shapes):
lons,lats = zip(*shape.points)
data = np.array(m(lons, lats)).T
if len(shape.parts) == 1:
segs = [data,]
else:
segs = []
for i in range(1,len(shape.parts)):
index = shape.parts[i-1]
index2 = shape.parts[i]
segs.append(data[index:index2])
segs.append(data[index2:])
lines = LineCollection(segs,antialiaseds=(1,))
lines.set_facecolors('#fafaf8')
lines.set_edgecolors('#333333')
lines.set_linewidth(0.5)
ax.add_collection(lines)
# Finish plot
plt.axis('off')
plt.savefig("test2.png",bbox_inches='tight',pad_inches=0);
你还没有绘制荷兰。
# Create map
...
# Draw contours
...
# Plot seas from shapefile
sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/northsea')
...
# Plot France from shapefile
sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/FRA_adm0')
...
# Plot Belgium from shapefile
sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/BEL_adm0')
...
# Plot Germany from shapefile
sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/DEU_adm0')
...
# Finish plot
所以荷兰的部分地区填充了您绘制的数据，海岸线在那里（通过绘制海洋？），但该国其他地区不存在，因为您没有绘制它。
我会避开这种从每个 shapefile 单独绘制的方法，并尝试只为您的海岸线等使用底图。如果您必须从 shapefile 中执行此操作，我鼓励您定义一个函数来执行类似于：
def drawContentsOfShapeFile(mymap, path_to_shapefile):
<whatever>
#and use it:
m = Basemap(...)
myshapefiles = []
myshapefiles.append("/DIR_TO_SHP/shapefiles/FRA_adm0")
myshapefiles.append("/DIR_TO_SHP/shapefiles/norway")
for sf in myshapefiles: drawContentOfShapefile(m, sf)
你的问题是 LineCollection 没有按照你的想法去做。你想要的是一个带有 'holes punched in' 的外部填充形状，它具有 LineCollection 中其他线条的形状（即北海中的岛屿）。 LineCollection，然而，填充列表中的每个线段，所以填充的形状只是相互重叠。
Inspired by this post, I wrote an answer that seems to solve your problem using Patches. However, I'm not entirely sure how robust the solution is: according to the linked (unanswered) post, the vertices of the outer shape should be ordered clockwise while the vertices of the 'punched in' holes should be ordered counter-clock wise (I checked that too and it appears to be correct; this matplotlib example 显示相同的行为）。由于 shapefile 是二进制的，因此很难验证顶点的顺序，但结果似乎是正确的。在下面的示例中，我假设在每个 shapefile 中最长的线段是国家（或北海）的轮廓，而较短的线段是岛屿或类似的。因此，我首先按长度对每个 shapefile 的段进行排序，然后创建一个 Path 和一个 PathPatch。我自由地为每个 shapefile 使用不同的颜色，以确保一切正常。
import numpy as np
import matplotlib
matplotlib.use('Agg')
from scipy.interpolate import griddata
from mpl_toolkits.basemap import Basemap, maskoceans
import matplotlib.pyplot as plt
from numpy.random import seed
import shapefile as shp
from matplotlib.collections import LineCollection
from matplotlib.patches import Path, PathPatch
from matplotlib import cm
# Set figure size
fig, ax = plt.subplots(figsize=(15,15), dpi = 80)
# Define map bounds
xMin, xMax = 2.5, 8.0
yMin, yMax = 50.6, 53.8
shapefiles = [
'shapefiles/BEL_adm0',
'shapefiles/FRA_adm0',
'shapefiles/DEU_adm0',
'shapefiles/northsea',
]
colors = ['red', 'green', 'yellow', 'blue']
y = [54.325666666667,52.36,53.269444444444,55.399166666667,54.116666666667,53.614444444444,53.491666666667,53.824130555556,52.918055555556,54.03694,52.139722,52.926865008825,54.853888888889,52.317222,53.240026656696,52.642696895243,53.391265948394,52.505333893732,52.098821802977,52.896643913235,52.457270486008,53.223000488316,52.701902388132,52.0548617826,53.411581103636,52.434561756559,52.749056395511,53.123676213651,52.067534268959,53.194409573306,52.27314817052,51.441334059998,51.224757511326,51.990941918858,51.447744494043,51.960667359998,51.969031121385,51.564889021961,51.857593837453,51.449772459909,51.658528382201,51.196699902606,50.905256257898,51.497306260089,yMin,yMin,yMax,yMax]
x = [2.93575,3.3416666666667,3.6277777777778,3.8102777777778,4.0122222222222,4.9602777777778,5.9416666666667,2.9452777777778,4.1502777777778,6.04167,4.436389,4.7811453228565,4.6961111111111,4.789722,4.9207907082729,4.9787572406902,5.3458010937365,4.6029300588208,5.1797058644882,5.383478899702,5.5196324030324,5.7515738887123,5.8874461671401,5.8723225499118,6.1990994508938,6.2589770334531,6.5729701105864,6.5848470019087,6.6567253619722,7.1493220605216,6.8908745111116,3.5958241584686,3.8609657214986,4.121849767852,4.342014,4.4469005114756,4.9259216999194,4.9352386335384,5.1453989235756,5.3770039280214,5.7065946674719,5.7625447234516,5.7617834850481,6.1961067840608,xMin,xMax,xMin,xMax]
z = [4.8,5.2,5.8,5.4,5,5.3,5.4,4.6,5.8,6.3,4.8,5.4,5.3,4.6,5.4,4.4,4.1,5.5,4.5,4.2,3.9,3.7,4.2,3.2,4,3.8,2.7,2.3,3.4,2.5,3.7,5.2,2.9,5.1,3.8,4.4,4.2,3.9,3.8,3.2,2.6,2.8,2.4,3.1]
avg = np.average(z)
z.extend([avg,avg,avg,avg])
# Create map
m = Basemap(
ax = ax,
projection='merc',
llcrnrlon=xMin,
llcrnrlat=yMin,
urcrnrlon=xMax,
urcrnrlat=yMax,
resolution='h'
)
x,y = m(x,y)
m.drawmapboundary(fill_color='#d4dadc',linewidth=0.25)
# target grid to interpolate to
xis = np.arange(min(x),max(x),2000)
yis = np.arange(min(y),max(y),2000)
xi,yi = np.meshgrid(xis,yis)
# interpolate
zi = griddata((x,y),z,(xi,yi),method='cubic')
# Decide on proper values for colour bar (todo)
vrange = max(z)-min(z)
mult = 2
vmin = min(z)-(mult*vrange)
vmax = max(z)+(mult*vrange)
# Draw contours
cs = m.contour(xi, yi, zi, 5, linewidths=0.25, colors='k')
cs = m.contourf(xi, yi, zi, 5,vmax=vmax,vmin=vmin,cmap=plt.get_cmap('jet'))
for sf_name,color in zip(shapefiles, colors):
print(sf_name)
# Load data
#drawing shapes:
sf = shp.Reader(sf_name)
shapes = sf.shapes()
##print shapes[0].parts
records = sf.records()
##ax = plt.subplot(111)
for record, shape in zip(records,shapes):
lons,lats = zip(*shape.points)
data = np.array(m(lons, lats)).T
if len(shape.parts) == 1:
segs = [data,]
else:
segs = []
for i in range(1,len(shape.parts)):
index = shape.parts[i-1]
index2 = shape.parts[i]
seg = data[index:index2]
segs.append(seg)
segs.append(data[index2:])
##assuming that the longest segment is the enclosing
##line and ordering the segments by length:
lens=np.array([len(s) for s in segs])
order = lens.argsort()[::-1]
segs = [segs[i] for i in order]
##producing the outlines:
lines = LineCollection(segs,antialiaseds=(1,),zorder=4)
##note: leaving the facecolors out:
##lines.set_facecolors('#abc0d3')
lines.set_edgecolors('red')
lines.set_linewidth(0.5)
ax.add_collection(lines)
##producing a path from the line segments:
segs_lin = [v for s in segs for v in s]
codes = [
[Path.MOVETO]+
[Path.LINETO for p in s[1:]]
for s in segs]
codes_lin = [c for s in codes for c in s]
path = Path(segs_lin, codes_lin)
##patch = PathPatch(path, facecolor="#abc0d3", lw=0, zorder = 3)
patch = PathPatch(path, facecolor=color, lw=0, zorder = 3)
ax.add_patch(patch)
plt.axis('off')
fig.savefig("shapefiles.png",bbox_inches='tight',pad_inches=0)
结果如下所示：
希望对您有所帮助。