python-2.7 - Matplotlib-Basemap に似た 3D CartoPy

Question

Cartopy を 3D プロットで使用できるかどうかについて質問がある Python は初めてです。以下は を使用した例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()

これにより、3D 軸内にマップが作成され、サーフェス上にオブジェクトをプロットできるようになります。しかし、Cartopy ではmatplotlib.axes.GeoAxesSubplot. これを取得して、上記の のように 3D フィギュア/軸に追加する方法が明確ではありませんmatplotlib-basemap。

それで、Cartopy で同様の 3D プロットを行う方法について誰かが何か指針を与えることができますか?

Answer 1

ベースマップ mpl3d は非常に巧妙なハックですが、説明されている方法で機能するようには設計されていません。その結果、現在、単純な海岸線以外には同じ手法を使用することはできません。たとえば、塗りつぶされた大陸は AFAICT では機能しません。

とはいえ、cartopy を使用する場合にも同様のハックが利用できます。シェープファイル情報に一般的にアクセスできるため、このソリューションは海岸線などのポリライン シェープファイルで機能するはずです。

最初のステップは、シェープファイルとそれぞれのジオメトリを取得することです。

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

次に、これらを目的の投影に変換できます。

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

これらは形の整ったジオメトリであるため、次を使用して matplotlib パスに変換します。

from cartopy.mpl.patch import geos_to_path
import itertools

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

パスを使用すると、matplotlib で PathCollection を作成し、それを軸に追加するだけでよいはずですが、残念ながら、Axes3D は PathCollection インスタンスに対応していないようです。 ）。悲しいことに、LineCollections はパスではなく、次のように計算できるセグメントを取ります。

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

これをすべてまとめると、コードが生成するベースマップ プロットと同様の結果になります。

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()

これに加えて、mpl3d は PolyCollection を適切に処理しているように見えます。これは、土地の輪郭 (厳密には輪郭である海岸線とは対照的に) などの塗りつぶされたジオメトリを調査するルートになります。

重要なステップは、パスをポリゴンに変換し、これらを 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)

この場合の完全なコードは次のようになります。

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()

得た：