shapefile e matplotlib: raccolta poligonale trama delle coordinate dello shapefile

Question

Sto provando a tracciare poligoni pieni di paesi sulla mappa del mondo con matplotlib in python.

Ho uno shapefile con le coordinate del confine di ogni paese. Ora, voglio convertire queste coordinate (per ogni paese) in un poligono con matplotlib. Senza usare Basemap. Sfortunatamente, le parti si incrociano o si sovrappongono. C'è un workarund, magari usando la distanza da un punto all'altro ... o riordinandoli?

Answer 1

Ha! Ho scoperto, come ... ho completamente trascurato, le informazioni di sf.shapes [i] .parts! Poi si tratta di:

# -- import -- 
import shapefile 
import matplotlib.pyplot as plt 
import matplotlib.patches as patches 
from matplotlib.patches import Polygon 
from matplotlib.collections import PatchCollection 
# -- input -- 
sf = shapefile.Reader("./shapefiles/world_countries_boundary_file_world_2002") 
recs = sf.records() 
shapes = sf.shapes() 
Nshp = len(shapes) 
cns  = [] 
for nshp in xrange(Nshp): 
    cns.append(recs[nshp][1]) 
cns = array(cns) 
cm = get_cmap('Dark2') 
cccol = cm(1.*arange(Nshp)/Nshp) 
# -- plot -- 
fig  = plt.figure() 
ax  = fig.add_subplot(111) 
for nshp in xrange(Nshp): 
    ptchs = [] 
    pts  = array(shapes[nshp].points) 
    prt  = shapes[nshp].parts 
    par  = list(prt) + [pts.shape[0]] 
    for pij in xrange(len(prt)): 
    ptchs.append(Polygon(pts[par[pij]:par[pij+1]])) 
    ax.add_collection(PatchCollection(ptchs,facecolor=cccol[nshp,:],edgecolor='k', linewidths=.1)) 
ax.set_xlim(-180,+180) 
ax.set_ylim(-90,90) 
fig.savefig('test.png') 

allora sarà simile a questa:

Answer 2

Ecco un altro pezzo di codice che ho usato per tracciare shapefile poligonali. Esso utilizza GDAL/OGR per leggere shapefile e trame ciambella correttamente poligoni di forma:

from osgeo import ogr 
import numpy as np 
import matplotlib.path as mpath 
import matplotlib.patches as mpatches 
import matplotlib.pyplot as plt 

# Extract first layer of features from shapefile using OGR 
ds = ogr.Open('world_countries_boundary_file_world_2002.shp') 
nlay = ds.GetLayerCount() 
lyr = ds.GetLayer(0) 

# Get extent and calculate buffer size 
ext = lyr.GetExtent() 
xoff = (ext[1]-ext[0])/50 
yoff = (ext[3]-ext[2])/50 

# Prepare figure 
fig = plt.figure() 
ax = fig.add_subplot(111) 
ax.set_xlim(ext[0]-xoff,ext[1]+xoff) 
ax.set_ylim(ext[2]-yoff,ext[3]+yoff) 

paths = [] 
lyr.ResetReading() 

# Read all features in layer and store as paths 
for feat in lyr: 
    geom = feat.geometry() 
    codes = [] 
    all_x = [] 
    all_y = [] 
    for i in range(geom.GetGeometryCount()): 
     # Read ring geometry and create path 
     r = geom.GetGeometryRef(i) 
     x = [r.GetX(j) for j in range(r.GetPointCount())] 
     y = [r.GetY(j) for j in range(r.GetPointCount())] 
     # skip boundary between individual rings 
     codes += [mpath.Path.MOVETO] + \ 
        (len(x)-1)*[mpath.Path.LINETO] 
     all_x += x 
     all_y += y 
    path = mpath.Path(np.column_stack((all_x,all_y)), codes) 
    paths.append(path) 

# Add paths as patches to axes 
for path in paths: 
    patch = mpatches.PathPatch(path, \ 
      facecolor='blue', edgecolor='black') 
    ax.add_patch(patch) 

ax.set_aspect(1.0) 
plt.show() 

Answer 3

from fiona import collection 
import matplotlib.pyplot as plt 
from descartes import PolygonPatch 
from matplotlib.collections import PatchCollection 
from itertools import imap 
from matplotlib.cm import get_cmap 

cm = get_cmap('Dark2') 

figure, axes = plt.subplots(1) 

source_path = "./shapefiles/world_countries_boundary_file_world_2002" 

with collection(source_path, 'r') as source: 
    patches = imap(PolygonPatch, (record['geometry'] for record in source) 

axes.add_collection(PatchCollection (patches, cmap=cm, linewidths=0.1)) 

axes.set_xlim(-180,+180) 
axes.set_ylim(-90,90) 

plt.show() 

Nota questo presuppone poligoni, multipoligoni possono essere le maniglie in modo simile con

map(PolygonPatch, MultiPolygon(record['geometry'])) 

Answer 4

Riguardo al @hannesk ' s risposta, è necessario aggiungere le seguenti importazioni: from numpy import array e import matplotlib e sostituire la riga cm = get_cmap('Dark2') per cm = matplotlib.cm.get_cmap('Dark2')

(Non sono così famoso per aggiungere un commento al post notato.)