使用LineCollection绘制shapefile会显示所有边缘，但会部分填充它们 [英] Plotting shapefile using LineCollection shows all edges, but partially fills them

问题描述

在过去的几天里，我一直在尝试仅在针对我的国家的地图中插入气象站数据.我这样做如下:

For the past few days I have been trying to get weather station data interpolated in a map for my country only. I do this as follows:

• 加载数据，我使用插值创建一个网格
• 基于此网格，我绘制了contour和contourf图像
• 然后我在地图顶部绘制德国，比利时和法国的shapefile，以覆盖不相关的contour/contourf元素.我将本教程用于那.
• 最后，我还使用了一个海洋形状文件(IHO Sea Areas； www.marineregions.org/downloads.php#iho)进行绘图，以覆盖北海.使用QGIS，我编辑了这个海洋shapefile并删除了除北海以外的所有内容-给定时间限制:)

• Loading the data, I create a grid using interpolation
• Based on this grid I draw a contour and contourf image
• I then draw shapefiles for Germany, Belgium and France on top of the map in order to cover the irrelevant contour/contourf elements. I used this tutorial for that.
• Finally, I use an oceans shapefile (IHO Sea Areas; www.marineregions.org/downloads.php#iho) to plot as well in order to cover the North Sea. Using QGIS, I edited this oceans shapefile and removed everything except for the North Sea - given time constraints :)

您会说一切都很好-但由于某些原因，该国的部分地区以及这些岛屿被解释为是水.我猜这是因为这些都是截然不同的部分，都没有连接到主要土地(由于水域/河流).

You would say that everything goes fine -- but for some reason parts of the country as well as the islands are interpreted as being water. I guess this is because these are distinct parts, all not connected to the main land (due to waters/rivers).

奇怪的是，绘制了边缘，但未填充边缘.

Strangely enough, the edges are drawn, but they are not filled.

我已经尝试和搜索了很多，但是对如何解决这个问题一无所知.我猜它在LineCollection中的某个位置，因为在QGIS中shapefile是正确的(即，当单击岛等时，它识别出 no shape ，这是正确的，因为它仅在单击时才识别出形状大海). 我衷心希望您能帮助我指出哪里出了问题以及我该如何解决！

I've tried and searched a lot but have no clue about how to fix this. I guess it is somewhere in the LineCollection because in QGIS the shapefile is correct (i.e. it recognizes no shape when clicking the islands etc., which is correct since it should only recognize a shape when clicking on the sea). I sincerely hope that you could help me point out where I'm wrong and how I can fix this!

非常感谢！

这是我得到的地图:

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);

推荐答案

您的问题是LineCollection并未执行您认为的操作.您想要的是一个外部填充形状，带有打孔"，具有您的LineCollection中其他线条的形状(即北海中的岛屿).但是，LineCollection填充列表中的每个线段，所以填充的形状彼此重叠.

Your problem is that LineCollection does not do what you think it does. What you want is an outer filled shape with 'holes punched in' that have the shapes of the other lines in your LineCollection (i.e. the islands in the north sea). LineCollection, however, fills every line segment in the list, so the filled shapes just overlay each other.

受这篇文章的启发，我写了一个答案，似乎可以使用Patches解决您的问题.但是，我不能完全确定解决方案的鲁棒性:根据链接的(未答复)帖子，外形的顶点应按顺时针方向排序，而打孔"孔的顶点应按逆时针方向排序. (我也检查过，它似乎是正确的；

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 shows the same behaviour). As the shapefiles are binary, it is hard to verify the ordering of the vertices, but the result appears to be correct. In the example below I assume that in each shapefile the longest line segment is the outline of the country (or north sea), while shorter segments are islands or some such. I thus first order the segments of each shapefile by length and create a Path and a PathPatch thereafter. I took the freedom to use a different colour for each shapefile in order to make sure that everything works as it should.

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)

结果如下:

希望这会有所帮助.

这篇关于使用LineCollection绘制shapefile会显示所有边缘，但会部分填充它们的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！