2016-02-22 61 views
1

我需要屏蔽掉海洋/水面上的数据,以便只能看到陆地上的数据。以下是我正在使用的matplotlib脚本的示例。该数据通过griddata()进行插值。掩蔽海洋上的插值数据

的Python脚本:

def mapformat(): 

    m = Basemap(llcrnrlon=-85,llcrnrlat=36.5,urcrnrlon=-64,urcrnrlat=47.5, 
     projection='lcc',lat_1=33,lat_2=45,lon_0=-70,lon_1=-60, resolution='h') 
    # resolution c, l, i, h, f in that order 

    m.drawmapboundary(fill_color='aqua', zorder=-1) 
    m.fillcontinents(color='green', lake_color='aqua', zorder=0) 

    m.drawcounties(color='0.1', linewidth=0.05, antialiased=True) 
    m.drawcoastlines(color='0.0', linewidth=0.25, antialiased=True) 
    m.drawcountries(color='0.0', linewidth=0.5, antialiased=True) 
    m.drawstates(color='0.0', linewidth=0.25, antialiased=True) 
    #m.drawparallels(np.arange(35.,45.,5), labels=[1,0,0,1], dashes=[1,1], linewidth=0.25, color='0.5') 
    #m.drawmeridians(np.arange(0., 360., 5.), labels=[1,0,0,1], dashes=[1,1], linewidth=0.25, color='0.5') 

    return m 

data = np.loadtxt('/home/.../.../.../maxs', delimiter=',', skiprows=1) 

m = mapformat() 

dx = 0.25 

grid_x, grid_y = np.mgrid[-85:-60:dx, 34:50:dx] #Northeast 

temp = data[:,0] 

#print temp 

figure = plt.gcf() # get current figure 
figure.set_size_inches(8, 4.5) 

grid_z = griddata((data[:,2],data[:,1]), data[:,0], (grid_x,grid_y), method='cubic') 

x,y = m(data[:,2], data[:,1]) # flip lat/lon 

grid_x,grid_y = m(grid_x,grid_y) 

#m.plot(x,y, 'ko', markersize=2) 

for i in range(len(temp)): 
    fmt=r"%.f" % (temp[i]) 
    #plt.text(x[i], y[i], fmt, va="center", ha="center", fontsize='12') 
    plt.annotate(fmt,xy=(x[i], y[i]), xytext=None, va="center", ha="center", fontsize="3") 

    clevs1 =[-30,-29,-28,-27,-26,-25,-24,-23,-22,-21, -20, -19,-18,-17,-16,-15,-14,-13,-12,-11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18, 
19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57, 
58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97, 
98,99,100,101,102,103,104,105,106,107,108,109,110] 

custom = mcolors.LinearSegmentedColormap.from_list('own2', [ltsilver,medblue,darkpurple,brightpink,medred,white, 
medblue,cyan,medgreen,yellow,orange,red,darkred,darkpeach,ltsilver]) 

m.contourf(grid_x,grid_y,grid_z,clevs1,cmap=custom) 

地块样本:

ne.png

我已阅读,屏蔽掉可以通过在底图is_land可以实现,但我不知道它将与插值数据一起工作。另外,海洋/水域上没有数据点。

回答

0

我记得我试图做类似的事情。我会尽我所能,但这是来自挖掘旧代码。

的关键步骤似乎是:

  1. 获得了良好的陆/海边界
  2. 做一个多边形对象并将
  3. 添加多边形到轴
  4. 展现您的数据,使用关键字clip_pathclip_on

的特定代码段我是

USPatch = mpl.patch.Polygon(USXY, facecolor='none', edgecolor='none') 
ax.add_patch(USPatch) 
im = bm.pcolormesh(lonE, latE, data, 
        ax=ax, clip_path=USPatch, clip_on=True, zorder=1) 
im.set_clip_path(USPatch) 

其中USXYnumpy数组美国边境和bm是我的底图对象。我想我从bm.drawcoastlines().get_segments()得到了USXY,但我不知道(或不记得了)如何从线段去一个良好的,完整的边界。

不幸的是,我一直无法从单独的旧代码片断做出最小的工作示例。尝试解析this on Basemap and clipping