2017-04-14 4 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です。 meshgridには、ガウスの等高線を描画したい位置にポイントがありません。

この動作を再現する例は、x directioのmeshgridが-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)

ソリューションはもちろん、ポートランドにメッシュを開始しないことだろうが、どこかで海、すなわち、地図の端にある。

関連する問題