【问题标题】:contour lines from the edge of a map don't show up on basemap地图边缘的等高线未显示在底图上
【发布时间】:2017-09-10 10:17:29
【问题描述】:

我在底图投影上绘制了几条等高线,如下图所示:

有 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)

【问题讨论】:

    标签: matplotlib contour matplotlib-basemap


    【解决方案1】:

    问题是网格没有覆盖完整的地图。网格网格在您要绘制高斯轮廓线的位置根本没有任何点。

    以下是重现此行为的示例,其中 x 方向中的网格网格从 -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()
    

    问题的代码中出现了类似的问题。
    在 Y 方向上,您使用完整的贴图 Y = np.linspace(lllat, urlat, numrows),在 X 方向上,您将网格限制为从 mlon.min() 开始,

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

    解决方案当然不是在波特兰开始网格,而是在海洋的某个地方,即在所示地图的边缘。

    【讨论】:

      猜你喜欢
      • 2019-03-22
      • 2023-02-04
      • 1970-01-01
      • 2021-11-30
      • 1970-01-01
      • 2015-01-08
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多