如何使用matplotlib删除/省略较小的轮廓线

Question

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

绘图后图像看起来像这样

从图像中我已经包围了小而不需要绘制的轮廓.此外,我想contour如上图所示绘制所有线条更平滑.总的来说我想得到这样的轮廓图像: -

我想到的可能解决方案是

1. 找出绘制轮廓和掩模所需的点数/如果它们的数量很少,则忽略这些线.

要么

2. 找到轮廓的区域(因为我只想省略圆圈轮廓)并省略/掩盖那些较小的轮廓.

要么

3. 通过将距离增加到50 km - 100 km来降低分辨率(仅限轮廓).

我能够使用SO线程Python成功获得点数:从matplotlib.pyplot.contour()中查找轮廓线

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

任何实现上述建议解决方案的解决方案都非常受欢迎.

编辑:-

@ Andras Deak我print 'diameter is ', diameter在线上方使用了del(level.get_paths()[kp])一行来检查代码是否过滤掉了所需的直径.这是我设置时的过滤消息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()

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

您可以看到线条尚未移除.以下是mslp我们尝试使用的数组链接http://www.mediafire.com/download/7vi0mxqoe0y6pm9/mslp.txt

如果您需要x以及y上述代码中使用的数据,我可以上传以供您查看.

流畅的线条

您可以编码以删除较小的圆圈.然而,我在原帖(平滑线)中提出的另一个问题似乎不起作用.我已经使用您的代码对数组进行切片以获得最小值并对其进行轮廓分析.我使用以下代码来减少数组大小: -

slice = 15

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

结果如下.

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

重新定期的netcdf数据

但是那里提供的解决方案都没有起作用.上面类似于我的问题没有适当的解决方案.如果这个问题得到解决,那么代码就是完美而完整的.

Answer 1

大概的概念

你的问题似乎有两个非常不同的一半:一个是关于省略小轮廓,另一个关于平滑轮廓线.后者更简单,因为除了降低你的contour()通话分辨率之外我无法想到其他任何事情,就像你说的那样.

至于去除一些轮廓线,这里是一个基于直接单独去除轮廓线的解决方案.您必须遍历collections返回的对象contour(),并为每个元素检查每个元素Path,并删除不需要的元素.重绘画figure布将摆脱不必要的线条:

# 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:

总而言之,它似乎很合理.

你的实际案例

由于您已添加了实际数据,因此这是您的案例的应用程序.请注意,您可以level使用单行直接生成s np,这几乎可以为您提供相同的结果.完全相同的可在2条线来实现(产生arange,然后选择那些落入之间p1和p2).此外,由于您正在设置levels调用contour,我相信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'更快,但不太流畅.

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

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