El trazado de shapefile usando LineCollection muestra todos los bordes, pero los llena parcialmente

Durante los últimos días he estado tratando de interpolar los datos de la estación meteorológica en un mapa de mi país solamente. Hago esto de la siguiente manera:

  • Cargando los datos, creo una grilla utilizando interpolación.
  • Basándome en esta cuadrícula dibujo un contour y contourf imagen de contour contourf
  • Luego dibujo shapefiles para Alemania, Bélgica y Francia en la parte superior del mapa para cubrir los elementos de contour / contour irrelevantes. Utilicé este tutorial para eso.
  • Finalmente, uso un shapefile de océanos (Áreas Marinas de la OHI; www.marineregions.org/downloads.php#iho) para trazar también para cubrir el Mar del Norte. Usando QGIS, edité este shapefile de océanos y eliminé todo, excepto las restricciones de tiempo del Mar del Norte 🙂

Diría que todo va bien, pero, por alguna razón, partes del país y de las islas se interpretan como agua. Supongo que esto se debe a que estas son partes distintas, todas no conectadas a la tierra principal (debido a las aguas / ríos).

Curiosamente, los bordes están dibujados, pero no están rellenos.

Lo he intentado y buscado mucho, pero no tengo idea de cómo solucionarlo. Supongo que está en algún lugar de LineCollection porque en QGIS el shapefile es correcto (es decir, no reconoce ninguna forma al hacer clic en las islas, etc., lo cual es correcto ya que solo debería reconocer una forma al hacer clic en el mar). Espero sinceramente que puedas ayudarme a señalar dónde estoy equivocado y cómo puedo solucionarlo.

¡Muchas gracias!

Este es el mapa que estoy obteniendo:

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

Mi código es el siguiente ( y probablemente vea que soy muy nuevo en este tipo de progtwigción, comencé ayer 🙂 ):

 import numpy as np import matplotlib matplotlib.use('Agg') from scipy.interpolate import griddata from mpl_toolkits.basemap import Basemap, maskoceans import matplotlib.pyplot as plt from numpy.random import seed import shapefile as shp from matplotlib.collections import LineCollection from matplotlib import cm # Set figure size plt.figure(figsize=(15,15), dpi=80) # Define map bounds xMin, xMax = 2.5, 8.0 yMin, yMax = 50.6, 53.8 # Create map m = Basemap(projection='merc',llcrnrlon=xMin,llcrnrlat=yMin,urcrnrlon=xMax,urcrnrlat=yMax,resolution='h') m.drawmapboundary(fill_color='#d4dadc',linewidth=0.25) # m.drawcoastlines(linewidth=0.5,color='#333333') # Load data y = [54.325666666667,52.36,53.269444444444,55.399166666667,54.116666666667,53.614444444444,53.491666666667,53.824130555556,52.918055555556,54.03694,52.139722,52.926865008825,54.853888888889,52.317222,53.240026656696,52.642696895243,53.391265948394,52.505333893732,52.098821802977,52.896643913235,52.457270486008,53.223000488316,52.701902388132,52.0548617826,53.411581103636,52.434561756559,52.749056395511,53.123676213651,52.067534268959,53.194409573306,52.27314817052,51.441334059998,51.224757511326,51.990941918858,51.447744494043,51.960667359998,51.969031121385,51.564889021961,51.857593837453,51.449772459909,51.658528382201,51.196699902606,50.905256257898,51.497306260089,yMin,yMin,yMax,yMax] x = [2.93575,3.3416666666667,3.6277777777778,3.8102777777778,4.0122222222222,4.9602777777778,5.9416666666667,2.9452777777778,4.1502777777778,6.04167,4.436389,4.7811453228565,4.6961111111111,4.789722,4.9207907082729,4.9787572406902,5.3458010937365,4.6029300588208,5.1797058644882,5.383478899702,5.5196324030324,5.7515738887123,5.8874461671401,5.8723225499118,6.1990994508938,6.2589770334531,6.5729701105864,6.5848470019087,6.6567253619722,7.1493220605216,6.8908745111116,3.5958241584686,3.8609657214986,4.121849767852,4.342014,4.4469005114756,4.9259216999194,4.9352386335384,5.1453989235756,5.3770039280214,5.7065946674719,5.7625447234516,5.7617834850481,6.1961067840608,xMin,xMax,xMin,xMax] z = [4.8,5.2,5.8,5.4,5,5.3,5.4,4.6,5.8,6.3,4.8,5.4,5.3,4.6,5.4,4.4,4.1,5.5,4.5,4.2,3.9,3.7,4.2,3.2,4,3.8,2.7,2.3,3.4,2.5,3.7,5.2,2.9,5.1,3.8,4.4,4.2,3.9,3.8,3.2,2.6,2.8,2.4,3.1] avg = np.average(z) z.extend([avg,avg,avg,avg]) x,y = m(x,y) # target grid to interpolate to xis = np.arange(min(x),max(x),2000) yis = np.arange(min(y),max(y),2000) xi,yi = np.meshgrid(xis,yis) # interpolate zi = griddata((x,y),z,(xi,yi),method='cubic') # Decide on proper values for colour bar (todo) vrange = max(z)-min(z) mult = 2 vmin = min(z)-(mult*vrange) vmax = max(z)+(mult*vrange) # Draw contours cs = m.contour(xi, yi, zi, 5, linewidths=0.25, colors='k') cs = m.contourf(xi, yi, zi, 5,vmax=vmax,vmin=vmin,cmap=plt.get_cmap('jet')) # Plot seas from shapefile sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/northsea') shapes = sf.shapes() print shapes[0].parts records = sf.records() ax = plt.subplot(111) for record, shape in zip(records,shapes): lons,lats = zip(*shape.points) data = np.array(m(lons, lats)).T print len(shape.parts) if len(shape.parts) == 1: segs = [data,] else: segs = [] for i in range(1,len(shape.parts)): index = shape.parts[i-1] index2 = shape.parts[i] segs.append(data[index:index2]) segs.append(data[index2:]) lines = LineCollection(segs,antialiaseds=(1,),zorder=3) lines.set_facecolors('#abc0d3') lines.set_edgecolors('red') lines.set_linewidth(0.5) ax.add_collection(lines) # Plot France from shapefile sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/FRA_adm0') shapes = sf.shapes() records = sf.records() ax = plt.subplot(111) for record, shape in zip(records,shapes): lons,lats = zip(*shape.points) data = np.array(m(lons, lats)).T if len(shape.parts) == 1: segs = [data,] else: segs = [] for i in range(1,len(shape.parts)): index = shape.parts[i-1] index2 = shape.parts[i] segs.append(data[index:index2]) segs.append(data[index2:]) lines = LineCollection(segs,antialiaseds=(1,)) lines.set_facecolors('#fafaf8') lines.set_edgecolors('#333333') lines.set_linewidth(0.5) ax.add_collection(lines) # Plot Belgium from shapefile sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/BEL_adm0') shapes = sf.shapes() records = sf.records() ax = plt.subplot(111) for record, shape in zip(records,shapes): lons,lats = zip(*shape.points) data = np.array(m(lons, lats)).T if len(shape.parts) == 1: segs = [data,] else: segs = [] for i in range(1,len(shape.parts)): index = shape.parts[i-1] index2 = shape.parts[i] segs.append(data[index:index2]) segs.append(data[index2:]) lines = LineCollection(segs,antialiaseds=(1,)) lines.set_facecolors('#fafaf8') lines.set_edgecolors('#333333') lines.set_linewidth(0.5) ax.add_collection(lines) # Plot Germany from shapefile sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/DEU_adm0') shapes = sf.shapes() records = sf.records() ax = plt.subplot(111) for record, shape in zip(records,shapes): lons,lats = zip(*shape.points) data = np.array(m(lons, lats)).T if len(shape.parts) == 1: segs = [data,] else: segs = [] for i in range(1,len(shape.parts)): index = shape.parts[i-1] index2 = shape.parts[i] segs.append(data[index:index2]) segs.append(data[index2:]) lines = LineCollection(segs,antialiaseds=(1,)) lines.set_facecolors('#fafaf8') lines.set_edgecolors('#333333') lines.set_linewidth(0.5) ax.add_collection(lines) # Finish plot plt.axis('off') plt.savefig("test2.png",bbox_inches='tight',pad_inches=0); 

Su problema es que LineCollection no hace lo que usted cree que hace. Lo que desea es una forma rellena exterior con ‘agujeros perforados’ que tengan las formas de las otras líneas en su LineCollection (es decir, las islas en el mar del norte). LineCollection , sin embargo, llena cada segmento de línea en la lista, por lo que las formas rellenas simplemente se superponen entre sí.

Inspirado por esta publicación , escribí una respuesta que parece resolver su problema con los Patches . Sin embargo, no estoy completamente seguro de cuán robusta es la solución: según el poste vinculado (sin respuesta), los vértices de la forma exterior deben ordenarse en el sentido de las agujas del reloj, mientras que los vértices de los orificios “perforados” deben ordenarse en sentido contrario a las agujas del reloj (También verifiqué eso y parece ser correcto; este ejemplo de matplotlib muestra el mismo comportamiento). Como los shapefiles son binarios, es difícil verificar el orden de los vértices, pero el resultado parece ser correcto. En el siguiente ejemplo, asumo que en cada shapefile el segmento de línea más largo es el perfil del país (o mar del norte), mientras que los segmentos más cortos son islas o algo similar. Por lo tanto, primero ordeno los segmentos de cada shapefile por longitud y luego creo un Path y un PathPatch . Me tomé la libertad de usar un color diferente para cada shapefile para asegurarme de que todo funcione como debería.

 import numpy as np import matplotlib matplotlib.use('Agg') from scipy.interpolate import griddata from mpl_toolkits.basemap import Basemap, maskoceans import matplotlib.pyplot as plt from numpy.random import seed import shapefile as shp from matplotlib.collections import LineCollection from matplotlib.patches import Path, PathPatch from matplotlib import cm # Set figure size fig, ax = plt.subplots(figsize=(15,15), dpi = 80) # Define map bounds xMin, xMax = 2.5, 8.0 yMin, yMax = 50.6, 53.8 shapefiles = [ 'shapefiles/BEL_adm0', 'shapefiles/FRA_adm0', 'shapefiles/DEU_adm0', 'shapefiles/northsea', ] colors = ['red', 'green', 'yellow', 'blue'] y = [54.325666666667,52.36,53.269444444444,55.399166666667,54.116666666667,53.614444444444,53.491666666667,53.824130555556,52.918055555556,54.03694,52.139722,52.926865008825,54.853888888889,52.317222,53.240026656696,52.642696895243,53.391265948394,52.505333893732,52.098821802977,52.896643913235,52.457270486008,53.223000488316,52.701902388132,52.0548617826,53.411581103636,52.434561756559,52.749056395511,53.123676213651,52.067534268959,53.194409573306,52.27314817052,51.441334059998,51.224757511326,51.990941918858,51.447744494043,51.960667359998,51.969031121385,51.564889021961,51.857593837453,51.449772459909,51.658528382201,51.196699902606,50.905256257898,51.497306260089,yMin,yMin,yMax,yMax] x = [2.93575,3.3416666666667,3.6277777777778,3.8102777777778,4.0122222222222,4.9602777777778,5.9416666666667,2.9452777777778,4.1502777777778,6.04167,4.436389,4.7811453228565,4.6961111111111,4.789722,4.9207907082729,4.9787572406902,5.3458010937365,4.6029300588208,5.1797058644882,5.383478899702,5.5196324030324,5.7515738887123,5.8874461671401,5.8723225499118,6.1990994508938,6.2589770334531,6.5729701105864,6.5848470019087,6.6567253619722,7.1493220605216,6.8908745111116,3.5958241584686,3.8609657214986,4.121849767852,4.342014,4.4469005114756,4.9259216999194,4.9352386335384,5.1453989235756,5.3770039280214,5.7065946674719,5.7625447234516,5.7617834850481,6.1961067840608,xMin,xMax,xMin,xMax] z = [4.8,5.2,5.8,5.4,5,5.3,5.4,4.6,5.8,6.3,4.8,5.4,5.3,4.6,5.4,4.4,4.1,5.5,4.5,4.2,3.9,3.7,4.2,3.2,4,3.8,2.7,2.3,3.4,2.5,3.7,5.2,2.9,5.1,3.8,4.4,4.2,3.9,3.8,3.2,2.6,2.8,2.4,3.1] avg = np.average(z) z.extend([avg,avg,avg,avg]) # Create map m = Basemap( ax = ax, projection='merc', llcrnrlon=xMin, llcrnrlat=yMin, urcrnrlon=xMax, urcrnrlat=yMax, resolution='h' ) x,y = m(x,y) m.drawmapboundary(fill_color='#d4dadc',linewidth=0.25) # target grid to interpolate to xis = np.arange(min(x),max(x),2000) yis = np.arange(min(y),max(y),2000) xi,yi = np.meshgrid(xis,yis) # interpolate zi = griddata((x,y),z,(xi,yi),method='cubic') # Decide on proper values for colour bar (todo) vrange = max(z)-min(z) mult = 2 vmin = min(z)-(mult*vrange) vmax = max(z)+(mult*vrange) # Draw contours cs = m.contour(xi, yi, zi, 5, linewidths=0.25, colors='k') cs = m.contourf(xi, yi, zi, 5,vmax=vmax,vmin=vmin,cmap=plt.get_cmap('jet')) for sf_name,color in zip(shapefiles, colors): print(sf_name) # Load data #drawing shapes: sf = shp.Reader(sf_name) shapes = sf.shapes() ##print shapes[0].parts records = sf.records() ##ax = plt.subplot(111) for record, shape in zip(records,shapes): lons,lats = zip(*shape.points) data = np.array(m(lons, lats)).T if len(shape.parts) == 1: segs = [data,] else: segs = [] for i in range(1,len(shape.parts)): index = shape.parts[i-1] index2 = shape.parts[i] seg = data[index:index2] segs.append(seg) segs.append(data[index2:]) ##assuming that the longest segment is the enclosing ##line and ordering the segments by length: lens=np.array([len(s) for s in segs]) order = lens.argsort()[::-1] segs = [segs[i] for i in order] ##producing the outlines: lines = LineCollection(segs,antialiaseds=(1,),zorder=4) ##note: leaving the facecolors out: ##lines.set_facecolors('#abc0d3') lines.set_edgecolors('red') lines.set_linewidth(0.5) ax.add_collection(lines) ##producing a path from the line segments: segs_lin = [v for s in segs for v in s] codes = [ [Path.MOVETO]+ [Path.LINETO for p in s[1:]] for s in segs] codes_lin = [c for s in codes for c in s] path = Path(segs_lin, codes_lin) ##patch = PathPatch(path, facecolor="#abc0d3", lw=0, zorder = 3) patch = PathPatch(path, facecolor=color, lw=0, zorder = 3) ax.add_patch(patch) plt.axis('off') fig.savefig("shapefiles.png",bbox_inches='tight',pad_inches=0) 

El resultado se ve así:

resultado del código anterior

Espero que esto ayude.

No has trazado los Países Bajos.

 # Create map ... # Draw contours ... # Plot seas from shapefile sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/northsea') ... # Plot France from shapefile sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/FRA_adm0') ... # Plot Belgium from shapefile sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/BEL_adm0') ... # Plot Germany from shapefile sf = shp.Reader(r'/DIR_TO_SHP/shapefiles/DEU_adm0') ... # Finish plot