3D CartoPy simile a Matplotlib-Basemap

Question

Sono nuovo di Python con una domanda su Cartopy che può essere utilizzato in una trama 3D. Di seguito è riportato un esempio utilizzando matplotlibBasemap.

import matplotlib.pyplot as plt 
from mpl_toolkits.mplot3d import Axes3D 
from mpl_toolkits.basemap import Basemap 

m = Basemap(projection='merc', 
      llcrnrlat=52.0,urcrnrlat=58.0, 
      llcrnrlon=19.0,urcrnrlon=40.0, 
      rsphere=6371200.,resolution='h',area_thresh=10) 

fig = plt.figure() 
ax = Axes3D(fig) 
ax.add_collection3d(m.drawcoastlines(linewidth=0.25)) 
ax.add_collection3d(m.drawcountries(linewidth=0.35)) 
ax.add_collection3d(m.drawrivers(color='blue')) 

ax.set_xlabel('X') 
ax.set_ylabel('Y') 
ax.set_zlabel('Height') 

fig.show() 

Questo crea una mappa all'interno di un asse 3D in modo che sia possibile tracciare gli oggetti sulla superficie. Ma con Cartopy restituisce un matplotlib.axes.GeoAxesSubplot. Non è chiaro come prendere questo e aggiungere a una figura/asse 3D come sopra con matplotlib-basemap.

Quindi, qualcuno può dare qualche suggerimento su come fare una trama 3D simile con Cartopy?

Answer 1

La mappa di base mpl3d è un hack abbastanza curato, ma non è stato progettato per funzionare nel modo descritto. Di conseguenza, al momento non è possibile utilizzare la stessa tecnica per molto altro che semplici litorali. Ad esempio, i continenti pieni semplicemente non funzionano AFAICT.

Detto questo, un analogo hack è disponibile quando si utilizza il cartopy. Dal momento che possiamo accedere genericamente a informazioni sullo shapefile, questa soluzione dovrebbe funzionare per qualsiasi shapefile di polilinea come le linee costiere.

Il primo passo è procurarsi shapefile, e le rispettive geometrie:

feature = cartopy.feature.NaturalEarthFeature('physical', 'coastline', '110m') 
geoms = feature.geometries() 

Successivamente, possiamo convertire questi per la proiezione desiderata:

target_projection = ccrs.PlateCarree() 
geoms = [target_projection.project_geometry(geom, feature.crs) 
     for geom in geoms] 

Trattandosi di geometrie tornite , quindi vogliamo convertirli in percorsi matplotlib con:

from cartopy.mpl.patch import geos_to_path 
import itertools 

paths = list(itertools.chain.from_iterable(geos_to_path(geom) 
              for geom in geoms)) 

Con pa Dovremmo essere in grado di creare semplicemente un PathCollection in matplotlib e aggiungerlo agli assi, ma purtroppo Axes3D non sembra affrontare le istanze di PathCollection, quindi dobbiamo aggirarlo creando un LineCollection (come fa la mappa di base) . Purtroppo LineCollections non prendono percorsi, ma segmenti, che siamo in grado di calcolare con:

segments = [] 
for path in paths: 
    vertices = [vertex for vertex, _ in path.iter_segments()] 
    vertices = np.asarray(vertices) 
    segments.append(vertices) 

Tirando questo tutti insieme, si finisce con un risultato simile alla trama mappa di base che il codice produce:

import itertools 

from mpl_toolkits.mplot3d import Axes3D 
import matplotlib.pyplot as plt 
from matplotlib.collections import LineCollection 
import numpy as np 

import cartopy.feature 
from cartopy.mpl.patch import geos_to_path 
import cartopy.crs as ccrs 


fig = plt.figure() 
ax = Axes3D(fig, xlim=[-180, 180], ylim=[-90, 90]) 
ax.set_zlim(bottom=0) 


target_projection = ccrs.PlateCarree() 

feature = cartopy.feature.NaturalEarthFeature('physical', 'coastline', '110m') 
geoms = feature.geometries() 

geoms = [target_projection.project_geometry(geom, feature.crs) 
     for geom in geoms] 

paths = list(itertools.chain.from_iterable(geos_to_path(geom) for geom in geoms)) 

# At this point, we start working around mpl3d's slightly broken interfaces. 
# So we produce a LineCollection rather than a PathCollection. 
segments = [] 
for path in paths: 
    vertices = [vertex for vertex, _ in path.iter_segments()] 
    vertices = np.asarray(vertices) 
    segments.append(vertices) 

lc = LineCollection(segments, color='black') 

ax.add_collection3d(lc) 

ax.set_xlabel('X') 
ax.set_ylabel('Y') 
ax.set_zlabel('Height') 

plt.show() 

in cima a questo, mpl3d sembra gestire PolyCollection bene, che sarebbe il percorso che avrebbe indagato per geometrie pieni, come ad esempio il contorno terreno (in contrapposizione alla linea di costa, che è strettamente un contorno).

Il passo importante è quello di convertire i percorsi di poligoni, e utilizzare questi in un oggetto PolyCollection:

concat = lambda iterable: list(itertools.chain.from_iterable(iterable)) 

polys = concat(path.to_polygons() for path in paths) 
lc = PolyCollection(polys, edgecolor='black', 
        facecolor='green', closed=False) 

il codice completo per questo caso sarebbe un aspetto simile:

import itertools 

from mpl_toolkits.mplot3d import Axes3D 
import matplotlib.pyplot as plt 
from matplotlib.collections import LineCollection, PolyCollection 
import numpy as np 

import cartopy.feature 
from cartopy.mpl.patch import geos_to_path 
import cartopy.crs as ccrs 


fig = plt.figure() 
ax = Axes3D(fig, xlim=[-180, 180], ylim=[-90, 90]) 
ax.set_zlim(bottom=0) 


concat = lambda iterable: list(itertools.chain.from_iterable(iterable)) 

target_projection = ccrs.PlateCarree() 

feature = cartopy.feature.NaturalEarthFeature('physical', 'land', '110m') 
geoms = feature.geometries() 

geoms = [target_projection.project_geometry(geom, feature.crs) 
     for geom in geoms] 

paths = concat(geos_to_path(geom) for geom in geoms) 

polys = concat(path.to_polygons() for path in paths) 

lc = PolyCollection(polys, edgecolor='black', 
        facecolor='green', closed=False) 

ax.add_collection3d(lc) 

ax.set_xlabel('X') 
ax.set_ylabel('Y') 
ax.set_zlabel('Height') 

plt.show() 

a cedere :