2017-04-14 190 views
2

我在底圖投影上繪製了幾條輪廓線,如下圖所示:enter image description here地圖邊緣的輪廓線不顯示在底圖上

有3個輪廓沒有完全繪製(在俄勒岡州,華盛頓州和加利福尼亞州),似乎這條線已經在相同的緯度上全部削減了3條。我不知道如何解決這個問題。 我加了插補點數,沒有幫助。改變了ll和ur點,包括更多的區域沒有幫助。

的代碼如下(不可再生,但可能會幫助):

def visualise_bigaus(mus, sigmas, corxys , output_type='pdf', **kwargs): 





    lllat = 24.396308 
    lllon = -124.848974 
    urlat = 49.384358 
    urlon = -66.885444 

    fig = plt.figure(figsize=(4, 2.5)) 
    ax = fig.add_subplot(111, axisbg='w', frame_on=False) 
    m = Basemap(llcrnrlat=lllat, 
    urcrnrlat=urlat, 
    llcrnrlon=lllon, 
    urcrnrlon=urlon, 
    resolution='i', projection='cyl') 

    m.drawmapboundary(fill_color = 'white') 
    #m.drawcoastlines(linewidth=0.2) 
    m.drawcountries(linewidth=0.2) 
    m.drawstates(linewidth=0.2, color='lightgray') 
    #m.fillcontinents(color='white', lake_color='#0000ff', zorder=2) 
    #m.drawrivers(color='#0000ff') 
    m.drawlsmask(land_color='gray',ocean_color="#b0c4de", lakes=True) 
    lllon, lllat = m(lllon, lllat) 
    urlon, urlat = m(urlon, urlat) 
    mlon, mlat = m(*(mus[:,1], mus[:,0])) 
    numcols, numrows = 1000, 1000 
    X = np.linspace(mlon.min(), urlon, numcols) 
    Y = np.linspace(lllat, urlat, numrows) 

    X, Y = np.meshgrid(X, Y) 
    m.scatter(mlon, mlat, s=0.2, c='red') 
    shp_info = m.readshapefile('./data/us_states_st99/st99_d00','states',drawbounds=True, zorder=0) 
    printed_names = [] 
    ax = plt.gca() 
    ax.xaxis.set_visible(False) 
    ax.yaxis.set_visible(False) 
    for spine in ax.spines.itervalues(): 
     spine.set_visible(False) 


    for k in xrange(mus.shape[0]): 
     #here x is longitude and y is latitude 
     #apply softplus to sigmas (to make them positive) 
     sigmax=np.log(1 + np.exp(sigmas[k][1])) 
     sigmay=np.log(1 + np.exp(sigmas[k][0])) 
     mux=mlon[k] 
     muy=mlat[k] 
     corxy = corxys[k] 
     #apply the soft sign 
     corxy = corxy/(1 + np.abs(corxy)) 
     #now given corxy find sigmaxy 
     sigmaxy = corxy * sigmax * sigmay 

     #corxy = 1.0/(1 + np.abs(sigmaxy)) 
     Z = mlab.bivariate_normal(X, Y, sigmax=sigmax, sigmay=sigmay, mux=mux, muy=muy, sigmaxy=sigmaxy) 

     #Z = maskoceans(X, Y, Z) 


     con = m.contour(X, Y, Z, levels=[0.02], linewidths=0.5, colors='darkorange', antialiased=True) 
     ''' 
     num_levels = len(con.collections) 
     if num_levels > 1: 
      for i in range(0, num_levels): 
       if i != (num_levels-1): 
        con.collections[i].set_visible(False) 
     ''' 
     contour_labels = False 
     if contour_labels: 
      plt.clabel(con, [con.levels[-1]], inline=True, fontsize=10) 

    ''' 
    world_shp_info = m.readshapefile('./data/CNTR_2014_10M_SH/Data/CNTR_RG_10M_2014','world',drawbounds=False, zorder=100) 
    for shapedict,state in zip(m.world_info, m.world): 
     if shapedict['CNTR_ID'] not in ['CA', 'MX']: continue 
     poly = MplPolygon(state,facecolor='gray',edgecolor='gray') 
     ax.add_patch(poly) 
    '''     
    if iter: 
     iter = str(iter).zfill(3) 
    else: 
     iter = '' 
    plt.tight_layout() 
    plt.savefig('./maps/video/gaus_' + iter + '.' + output_type, frameon=False, dpi=200) 

回答

2

的問題不是覆蓋整個地圖的meshgrid。網格在您想要繪製高斯輪廓線的位置根本沒有任何點。

重現此行爲的示例如下,其中x directio中的網格網格起始於-1,使得低於該網格點的點不會繪製。

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

fig, ax=plt.subplots() 
ax.plot([-2,2],[-2,-2], alpha=0) 
X,Y = np.meshgrid(np.linspace(-1,2),np.linspace(-2,2)) 
Z = mlab.bivariate_normal(X, Y, sigmax=1., sigmay=1., mux=0.1, muy=0.1, sigmaxy=0) 
con = ax.contour(X, Y, Z, levels=[Z.max()/3, Z.max()/2., Z.max()*0.8],colors='darkorange') 
plt.show() 

enter image description here

類似的問題發生在從所述問題的代碼。
而在Y方向,您可以使用完整的地圖,Y = np.linspace(lllat, urlat, numrows),在X方向你限制網在mlon.min()開始,

X = np.linspace(mlon.min(), urlon, numcols)

的解決辦法當然也可以不啓動在波特蘭的網格,但地方在海洋中,即在所示地圖的邊緣。