Cómo utilizar la ruta de acceso recortada para el polígono de mapa base

Quiero usar imshow (por ejemplo) para mostrar algunos datos dentro de los límites de un país (para los fines de ejemplo, elegí los EE. UU.) El siguiente ejemplo ilustra lo que quiero:

import numpy as np import matplotlib.pyplot as plt from matplotlib.patches import RegularPolygon data = np.arange(100).reshape(10, 10) fig = plt.figure() ax = fig.add_subplot(111) im = ax.imshow(data) poly = RegularPolygon([ 0.5, 0.5], 6, 0.4, fc='none', ec='k', transform=ax.transAxes) im.set_clip_path(poly) ax.add_patch(poly) ax.axis('off') plt.show() 

El resultado es:

introduzca la descripción de la imagen aquí

Ahora quiero hacer esto, pero en lugar de un simple polígono, quiero usar la forma compleja de los EE. UU. He creado algunos datos de ejemplo contenidos en la matriz de “Z” como se puede ver en el código a continuación. Es esta información la que quiero mostrar, usando un mapa de colores pero solo dentro de los límites de los Estados Unidos continentales.

Hasta ahora he intentado lo siguiente. Desde aquí obtengo un archivo de forma contenido en “nationp010g.shp.tar.gz” y uso el módulo de mapa base en Python para trazar los EE. UU. Tenga en cuenta que este es el único método que he encontrado que me permite obtener un polígono del área que necesito. Si hay métodos alternativos también me interesaría. Luego creo un polígono llamado “mainpoly” que es casi el polígono que quiero coloreado en azul:

    introduzca la descripción de la imagen aquí

    Observe que solo un cuerpo ha sido coloreado, todos los otros polígonos separados permanecen blancos:

    introduzca la descripción de la imagen aquí

    Entonces, el área de color azul es casi lo que quiero, tenga en cuenta que hay límites no deseados cerca de Canadá porque la frontera en realidad atraviesa algunos lagos, pero ese es un problema menor. El problema real es, ¿por qué no se muestran mis datos de imshow dentro de los Estados Unidos? Al comparar mis códigos de ejemplo primero y segundo, no puedo ver por qué no obtengo una imagen recortada en mi segundo ejemplo, como lo hago en el primero. Cualquier ayuda sería apreciada en entender lo que me estoy perdiendo.

     import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.basemap import Basemap as Basemap from matplotlib.patches import Polygon # Lambert Conformal map of lower 48 states. m = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49, projection='lcc',lat_1=33,lat_2=45,lon_0=-95) shp_info = m.readshapefile('nationp010g/nationp010g', 'borders', drawbounds=True) # draw country boundaries. for nshape,seg in enumerate(m.borders): if nshape == 1873: #This nshape denotes the large continental body of the USA, which we want mainseg = seg mainpoly = Polygon(mainseg,facecolor='blue',edgecolor='k') nx, ny = 10, 10 lons, lats = m.makegrid(nx, ny) # get lat/lons of ny by nx evenly space grid. x, y = m(lons, lats) # compute map proj coordinates. Z = np.zeros((nx,ny)) Z[:] = np.NAN for i in np.arange(len(x)): for j in np.arange(len(y)): Z[i,j] = x[0,i] ax = plt.gca() im = ax.imshow(Z, cmap = plt.get_cmap('coolwarm') ) im.set_clip_path(mainpoly) ax.add_patch(mainpoly) plt.show() 

    Actualizar

    Me doy cuenta de que la línea

     ax.add_patch(mainpoly) 

    Ni siquiera agrega la forma de polígono a un gráfico. ¿No lo estoy usando correctamente? Que yo sepa, mainpoly se calculó correctamente utilizando el método Polygon (). Verifiqué que las entradas de coordenadas son sensibles:

     plt.plot(mainseg[:,0], mainseg[:,1] ,'.') 

    lo que da

    introduzca la descripción de la imagen aquí

    También he considerado acerca de este problema durante mucho tiempo.
    Y encontré que el lenguaje NCL tiene la función de enmascarar los datos fuera de algún borde.
    Aquí está el ejemplo:

    http://sofes.miximages.com/python/bdb1a6c007b82645.png

    La ttwig de contourf solo se muestra dentro de la frontera china. Haga clic aquí para el código.

    Sé que Python tiene un paquete llamado PyNCL que admite todo el código NCL en el marco de Python.
    Pero realmente quiero trazar este tipo de figura utilizando el mapa base. Si lo has descubierto, por favor publica en internet. Lo aprenderé a la primera.

    ¡Gracias!

    Añadir 2016-01-16

    En cierto modo, lo he descubierto.
    Esta es mi idea y código, y está inspirada en esta pregunta que formulé hoy.

    Mi metodo
    1. Convierta el shapefile del área interesante (como US) en shapely.polygon.
    2. Pruebe cada punto de valor dentro / fuera del polígono.
    3. Si el punto de valor está fuera del área de estudio, ocúltelo como np.nan

    Introducción * el polígono xxx era una ciudad en China en formato de archivo shape ESRI. * Fiona, paquete bien formado fueron utilizados aquí.

     # generate the shapely.polygon shape = fiona.open("xxx.shp") pol = shape.next() geom = shape(pol['geometry']) poly_data = pol["geometry"]["coordinates"][0] poly = Polygon(poly_data) 

    Se muestra como:

    http://i4.tietuku.com/2012307faec02634.png

     ### test the value point ### generate the grid network which represented by the grid midpoints. lon_med = np.linspace((xi[0:2].mean()),(xi[-2:].mean()),len(x_grid)) lat_med = np.linspace((yi[0:2].mean()),(yi[-2:].mean()),len(y_grid)) value_test_mean = dsu.mean(axis = 0) value_mask = np.zeros(len(lon_med)*len(lat_med)).reshape(len(lat_med),len(lon_med)) for i in range(0,len(lat_med),1): for j in range(0,len(lon_med),1): points = np.array([lon_med[j],lat_med[i]]) mask = np.array([poly.contains(Point(points[0], points[1]))]) if mask == False: value_mask[i,j] = np.nan if mask == True: value_mask[i,j] = value_test_mean[i,j] # Mask the np.nan value Z_mask = np.ma.masked_where(np.isnan(so2_mask),so2_mask) # plot! fig=plt.figure(figsize=(6,4)) ax=plt.subplot() map = Basemap(llcrnrlon=x_map1,llcrnrlat=y_map1,urcrnrlon=x_map2,urcrnrlat=y_map2) map.drawparallels(np.arange(y_map1+0.1035,y_map2,0.2),labels= [1,0,0,1],size=14,linewidth=0,color= '#FFFFFF') lon_grid = np.linspace(x_map1,x_map2,len(x_grid)) lat_grid = np.linspace(y_map1,y_map2,len(y_grid)) xx,yy = np.meshgrid(lon_grid,lat_grid) pcol =plt.pcolor(xx,yy,Z_mask,cmap = plt.cm.Spectral_r ,alpha =0.75,zorder =2) 

    resultado

    http://i4.tietuku.com/c6620c5b6730a5f0.png

    http://i4.tietuku.com/a22ad484fee627b9.png

    resultado original

    http://i4.tietuku.com/011584fbc36222c9.png