2016-02-22 89 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