如何使用 matplotlib remove/omit 更小的轮廓线

How to remove/omit smaller contour lines using matplotlib

我正在尝试绘制 contour 压力水平线。我正在使用包含更高分辨率数据(范围从 3 公里到 27 公里)的 netCDF 文件。由于更高分辨率的数据集,我得到了很多不需要绘制的压力值(相反,我不介意省略某些无关紧要的值的等高线)。我已经根据 link http://matplotlib.org/basemap/users/examples.html 中给出的示例编写了一些绘图脚本。

绘制后的图像看起来像这样

从图像中我圈出了小的不需要绘制的等高线。此外,我想绘制所有 contour 线条,如上图所示更平滑。总的来说,我想得到这样的轮廓图像:-

我想到的可能的解决方案是

  1. 找出绘制等高线和mask/omit那些线条所需的点数,如果它们的数量很少。

  1. 找到等高线的面积(因为我只想省略带圆圈的等高线)并且 omit/mask 那些更小。

  1. 通过将距离增加到 50 公里 - 100 公里来降低分辨率(仅等高线)。

我能够使用SO线程成功获得积分Python: find contour lines from matplotlib.pyplot.contour()

但我无法使用这些要点实施上述任何建议的解决方案。

非常感谢任何实现上述建议解决方案的解决方案。

编辑:-

@安德拉斯迪克 我在 del(level.get_paths()[kp]) 行上方使用 print 'diameter is ', diameter 行来检查代码是否过滤掉了所需的直径。这是我设置 if diameter < 15000::

时过滤的消息
diameter is  9099.66295612
diameter is  13264.7838257
diameter is  445.574234531
diameter is  1618.74618114
diameter is  1512.58974168

但是生成的图像没有任何效果。所有看起来都与上面的照片相同。我很确定我已经保存了这个图形(在绘制风刺之后)。

关于降低分辨率的解决方案,plt.contour(x[::2,::2],y[::2,::2],mslp[::2,::2])有效。我必须应用一些过滤器来使曲线平滑。

用于删除行的完整工作示例代码:-

这是供您查看的示例代码

#!/usr/bin/env python
from netCDF4 import Dataset
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
import numpy as np
import scipy.ndimage
from mpl_toolkits.basemap import interp
from mpl_toolkits.basemap import Basemap


# Set default map
west_lon = 68
east_lon = 93
south_lat = 7
north_lat = 23

nc = Dataset('ncfile.nc')
# Get this variable for later calucation
temps = nc.variables['T2']
time = 0  # We will take only first interval for this example
# Draw basemap
m = Basemap(projection='merc', llcrnrlat=south_lat, urcrnrlat=north_lat,
                llcrnrlon=west_lon, urcrnrlon=east_lon, resolution='l')
m.drawcoastlines()
m.drawcountries(linewidth=1.0)
# This sets the standard grid point structure at full resolution
x, y = m(nc.variables['XLONG'][0], nc.variables['XLAT'][0])

# Set figure margins
width = 10
height = 8

plt.figure(figsize=(width, height))
plt.rc("figure.subplot", left=.001)
plt.rc("figure.subplot", right=.999)
plt.rc("figure.subplot", bottom=.001)
plt.rc("figure.subplot", top=.999)

plt.figure(figsize=(width, height), frameon=False)

# Convert Surface Pressure to Mean Sea Level Pressure
stemps = temps[time] + 6.5 * nc.variables['HGT'][time] / 1000.
mslp = nc.variables['PSFC'][time] * np.exp(9.81 / (287.0 * stemps) * nc.variables['HGT'][time]) * 0.01 + (
    6.7 * nc.variables['HGT'][time] / 1000)

# Contour only at 2 hpa interval
level = []
for i in range(mslp.min(), mslp.max(), 1):
    if i % 2 == 0:
        if i >= 1006 and i <= 1018:
            level.append(i)

# Save mslp values to upload to SO thread
# np.savetxt('mslp.txt', mslp, fmt='%.14f', delimiter=',')

P = plt.contour(x, y, mslp, V=2, colors='b', linewidths=2, levels=level)


# Solution suggested by Andras Deak
for level in P.collections:
    for kp,path in enumerate(level.get_paths()):
        # include test for "smallness" of your choice here:
        # I'm using a simple estimation for the diameter based on the
        #    x and y diameter...
        verts = path.vertices # (N,2)-shape array of contour line coordinates
        diameter = np.max(verts.max(axis=0) - verts.min(axis=0))

        if diameter < 15000: # threshold to be refined for your actual dimensions!
            #print 'diameter is ', diameter
            del(level.get_paths()[kp])  # no remove() for Path objects:(
            #level.remove() # This does not work. produces ValueError: list.remove(x): x not in list

plt.gcf().canvas.draw()


plt.savefig('dummy', bbox_inches='tight')
plt.close()

绘图保存后我得到相同的图像

可以看到线条还没有去掉。这是我们尝试使用 http://www.mediafire.com/download/7vi0mxqoe0y6pm9/mslp.txt

的 link 到 mslp 数组

如果您需要上述代码中使用的 xy 数据,我可以上传以供您查看。

平滑线

您编写代码以完美地移除较小的圆圈。然而,我在原来的post(平滑线)中提出的另一个问题似乎不起作用。我已经使用您的代码对数组进行切片以获得最小值并绘制轮廓。我使用以下代码来减小数组大小:-

slice = 15

CS = plt.contour(x[::slice,::slice],y[::slice,::slice],mslp[::slice,::slice], colors='b', linewidths=1, levels=levels)

结果如下。

搜索了几个小时后,我发现这个 SO 线程有类似的问题:-

Regridding regular netcdf data

但是none 那里提供的解决方案 works.The 与我上面类似的问题没有正确的解决方案。如果这个问题得到解决,那么代码就是完美和完整的。

这是一个非常糟糕的解决方案,但这是我提出的唯一解决方案。使用您链接到的 this solution 中的 get_contour_verts 函数,可能与 matplotlib._cntr 模块一起使用,以便最初不会绘制任何内容。这为您提供了等高线、截面、顶点等的列表。然后您必须浏览该列表和 pop 您不想要的等高线。例如,您可以通过计算最小直径来做到这一点;如果点之间的最大距离小于某个截止值,则将其丢弃。

这会给您留下 LineCollection 个对象的列表。 现在如果你制作一个FigureAxes实例,你可以使用Axes.add_collection来添加列表中的所有LineCollection .

我很快就检查过了,但它似乎有效。如果有机会,我会带着一个最低限度的工作示例回来。希望对您有所帮助!


编辑:这是基本思想的 MWE。我不熟悉 plt._cntr.Cntr,所以我最终使用 plt.contour 来获取初始轮廓对象。结果,您最终得到了两个数字;你只需要关闭第一个。您可以将 checkDiameter 替换为任何有效的函数。我 认为 你可以把线段变成 Polygon 并计算面积,但你必须自己弄清楚。如果您 运行 遇到此代码的问题,请告诉我,但它至少对我有用。

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

def checkDiameter(seg, tol=.3):
    # Function for screening line segments. NB: Not actually a proper diameter.
    diam = (seg[:,0].max() - seg[:,0].min(),
            seg[:,1].max() - seg[:,1].min())
    return not (diam[0] < tol or diam[1] < tol)

# Create testing data
x = np.linspace(-1,1, 21)
xx, yy = np.meshgrid(x,x)
z = np.exp(-(xx**2 + .5*yy**2))

# Original plot with plt.contour
fig0, ax0 = plt.subplots()
# Make sure this contour object actually has a tiny contour to remove
cntrObj = ax0.contour(xx,yy,z, levels=[.2,.4,.6,.8,.9,.95,.99,.999])

# Primary loop: Copy contours into a new LineCollection
lineNew = list()
for lineOriginal in cntrObj.collections:
    # Get properties of the original LineCollection
    segments = lineOriginal.get_segments()
    propDict = lineOriginal.properties()
    propDict = {key: value for (key,value) in propDict.items()
        if key in ['linewidth','color','linestyle']}  # Whatever parameters you want to carry over
    # Filter out the lines with small diameters
    segments = [seg for seg in segments if checkDiameter(seg)]
    # Create new LineCollection out of the OK segments
    if len(segments) > 0:
        lineNew.append(mpl.collections.LineCollection(segments, **propDict))

# Make new plot with only these line collections; display results
fig1, ax1 = plt.subplots()
ax1.set_xlim(ax0.get_xlim())
ax1.set_ylim(ax0.get_ylim())
for line in lineNew:
    ax1.add_collection(line)
plt.show()

仅供参考:带有propDict的位只是为了自动从原始图中引入一些线属性。但是,您不能一次使用整个词典。首先,它包含旧绘图的线段,但您可以将它们换成新的。但其次,它似乎包含许多相互冲突的参数:多个线宽、面色等。{key for key in propDict if I want key} 解决方法是我绕过它的方法,但我相信其他人可以做到更干净

总体思路

您的问题似乎有两个截然不同的部分:一个是关于省略小轮廓,另一个是关于平滑轮廓线。后者更简单,因为除了降低 contour() 调用的分辨率外,我真的想不出其他任何东西,就像你说的那样。

至于去除几条轮廓线,这里有一个基于直接逐条去除轮廓线的解决方案。您必须遍历 contour() 返回的对象的 collections,并为每个元素检查每个 Path,并删除不需要的元素。重新绘制 figure 的 canvas 将去掉不必要的行:

# dummy example based on matplotlib.pyplot.clabel example:
import matplotlib
import numpy as np
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt

delta = 0.025
x = np.arange(-3.0, 3.0, delta)
y = np.arange(-2.0, 2.0, delta)
X, Y = np.meshgrid(x, y)
Z1 = mlab.bivariate_normal(X, Y, 1.0, 1.0, 0.0, 0.0)
Z2 = mlab.bivariate_normal(X, Y, 1.5, 0.5, 1, 1)
# difference of Gaussians
Z = 10.0 * (Z2 - Z1)


plt.figure()
CS = plt.contour(X, Y, Z)

for level in CS.collections:
    for kp,path in reversed(list(enumerate(level.get_paths()))):
        # go in reversed order due to deletions!

        # include test for "smallness" of your choice here:
        # I'm using a simple estimation for the diameter based on the
        #    x and y diameter...
        verts = path.vertices # (N,2)-shape array of contour line coordinates
        diameter = np.max(verts.max(axis=0) - verts.min(axis=0))

        if diameter<1: # threshold to be refined for your actual dimensions!
            del(level.get_paths()[kp])  # no remove() for Path objects:(

# this might be necessary on interactive sessions: redraw figure
plt.gcf().canvas.draw()

这是直径阈值为 1 的原始版本(左)和删除版本(右)(注意顶部 0 级别的小部分):

请注意,顶部的小线已被删除,而中间的巨大青色线则没有,即使两者对应于相同的 collections 元素,即相同的轮廓级别。如果我们不想允许这样做,我们可以调用 CS.collections[k].remove(),这可能是做同样事情的更安全的方法(但它不允许我们区分对应于相同的等高线水平)。

为了表明摆弄截止直径能按预期工作,这是阈值 2 的结果:

总而言之,这似乎很合理。


您的实际案例

由于您已经添加了实际数据,下面是您案例的应用程序。请注意,您可以使用 np 在一行中直接生成 levels,这几乎可以得到相同的结果。可以在两行中实现完全相同的结果(生成 arange,然后选择介于 p1p2 之间的那些)。此外,由于您在 contour 的调用中设置了 levels,我相信函数调用的 V=2 部分没有效果。

import numpy as np
import matplotlib.pyplot as plt

# insert actual data here...
Z = np.loadtxt('mslp.txt',delimiter=',')
X,Y = np.meshgrid(np.linspace(0,300000,Z.shape[1]),np.linspace(0,200000,Z.shape[0]))
p1,p2 = 1006,1018

# this is almost the same as the original, although it will produce
# [p1, p1+2, ...] instead of `[Z.min()+n, Z.min()+n+2, ...]`
levels = np.arange(np.maximum(Z.min(),p1),np.minimum(Z.max(),p2),2)


#control
plt.figure()
CS = plt.contour(X, Y, Z, colors='b', linewidths=2, levels=levels)


#modified
plt.figure()
CS = plt.contour(X, Y, Z, colors='b', linewidths=2, levels=levels)

for level in CS.collections:
    for kp,path in reversed(list(enumerate(level.get_paths()))):
        # go in reversed order due to deletions!

        # include test for "smallness" of your choice here:
        # I'm using a simple estimation for the diameter based on the
        #    x and y diameter...
        verts = path.vertices # (N,2)-shape array of contour line coordinates
        diameter = np.max(verts.max(axis=0) - verts.min(axis=0))

        if diameter<15000: # threshold to be refined for your actual dimensions!
            del(level.get_paths()[kp])  # no remove() for Path objects:(

# this might be necessary on interactive sessions: redraw figure
plt.gcf().canvas.draw()
plt.show()

结果,原始(左)与新(右):


通过重采样平滑

我也决定解决平滑问题。我所能想到的就是对原始数据进行下采样,然后使用 griddata (插值)再次上采样。下采样部分也可以通过插值来完成,尽管输入数据中的小规模变化可能会使这个问题不适定。所以这是粗略的版本:

import scipy.interpolate as interp   #the new one

# assume you have X,Y,Z,levels defined as before

# start resampling stuff
dN = 10 # use every dN'th element of the gridded input data
my_slice = [slice(None,None,dN),slice(None,None,dN)]

# downsampled data
X2,Y2,Z2 = X[my_slice],Y[my_slice],Z[my_slice]
# same as X2 = X[::dN,::dN] etc.

# upsampling with griddata over original mesh
Zsmooth = interp.griddata(np.array([X2.ravel(),Y2.ravel()]).T,Z2.ravel(),(X,Y),method='cubic')

# plot
plt.figure()
CS = plt.contour(X, Y, Zsmooth, colors='b', linewidths=2, levels=levels)

您可以随意使用用于插值的网格,在这种情况下,我只使用手头的原始网格。您还可以尝试不同类型的插值:默认 'linear' 会更快,但不太流畅。

下采样(左)和上采样(右)后的结果:

当然你应该在这个重采样业务之后仍然应用小线去除算法,并记住这会严重扭曲你的输入数据(因为如果它没有扭曲,那么它就不会平滑).另外,请注意,由于在下采样步骤中使用了粗略的方法,我们在所考虑区域的 top/right 边缘附近引入了一些缺失值。如果这是一个问题,您应该考虑根据我之前提到的 griddata 进行下采样。