如何使用 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
线条,如上图所示更平滑。总的来说,我想得到这样的轮廓图像:-
我想到的可能的解决方案是
- 找出绘制等高线和mask/omit那些线条所需的点数,如果它们的数量很少。
或
- 找到等高线的面积(因为我只想省略带圆圈的等高线)并且 omit/mask 那些更小。
或
- 通过将距离增加到 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
数组
如果您需要上述代码中使用的 x
和 y
数据,我可以上传以供您查看。
平滑线
您编写代码以完美地移除较小的圆圈。然而,我在原来的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
个对象的列表。 现在如果你制作一个Figure
和Axes
实例,你可以使用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
在一行中直接生成 level
s,这几乎可以得到相同的结果。可以在两行中实现完全相同的结果(生成 arange
,然后选择介于 p1
和 p2
之间的那些)。此外,由于您在 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
进行下采样。
我正在尝试绘制 contour
压力水平线。我正在使用包含更高分辨率数据(范围从 3 公里到 27 公里)的 netCDF 文件。由于更高分辨率的数据集,我得到了很多不需要绘制的压力值(相反,我不介意省略某些无关紧要的值的等高线)。我已经根据 link http://matplotlib.org/basemap/users/examples.html 中给出的示例编写了一些绘图脚本。
绘制后的图像看起来像这样
从图像中我圈出了小的不需要绘制的等高线。此外,我想绘制所有 contour
线条,如上图所示更平滑。总的来说,我想得到这样的轮廓图像:-
我想到的可能的解决方案是
- 找出绘制等高线和mask/omit那些线条所需的点数,如果它们的数量很少。
或
- 找到等高线的面积(因为我只想省略带圆圈的等高线)并且 omit/mask 那些更小。
或
- 通过将距离增加到 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
数组
如果您需要上述代码中使用的 x
和 y
数据,我可以上传以供您查看。
平滑线
您编写代码以完美地移除较小的圆圈。然而,我在原来的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
个对象的列表。 现在如果你制作一个Figure
和Axes
实例,你可以使用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
在一行中直接生成 level
s,这几乎可以得到相同的结果。可以在两行中实现完全相同的结果(生成 arange
,然后选择介于 p1
和 p2
之间的那些)。此外,由于您在 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
进行下采样。