2017-04-14 4 views
2

Ich zeichne mehrere Konturlinien über eine Grundkartenprojektion, wie in der folgenden Abbildung gezeigt: enter image description here.Konturlinien vom Rand einer Karte werden nicht in der Grundkarte angezeigt

Es gibt 3 Konturen, die nicht vollständig gezeichnet sind (in Oregon, Washington und Kalifornien) und scheint, als gäbe es diese Linie, die alle 3 von ihnen in der gleichen Breite geschnitten hat. Ich bin mir nicht sicher, wie ich dieses Problem lösen soll. Ich habe die Anzahl der Stützpunkte hinzugefügt, hat nicht geholfen. änderte die ll und ur Punkte, um mehr Fläche zu enthalten, hat nicht geholfen.

Der Code ist unten (nicht reproduzierbar aber könnte helfen):

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) 

Antwort

2

Das Problem der meshgrid ist die komplette Karte nicht abdeckt. Das Meshgrid hat einfach keine Punkte an den Positionen, an denen Sie die Gaußsche Konturlinie zeichnen möchten.

Ein Beispiel, um dieses Verhalten zu reproduzieren, ist das folgende, wo das Meshgrid in x directio bei -1 beginnt, so dass Punkte niedriger als das nicht gezeichnet werden.

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

Ein ähnliches Problem tritt von der Frage im Code.
Während in Y-Richtung, verwenden Sie die komplette Karte, Y = np.linspace(lllat, urlat, numrows), in X-Richtung Sie das Netz beschränken auf mlon.min() zu starten,

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

Die Lösung wäre natürlich nicht das Netz in Portland zu starten, aber irgendwo im Ozean, also am Rand der gezeigten Karte.

Verwandte Themen