3D CartoPy类似于Matplotlib-Basemap

6za6bjd0  于 2023-05-01  发布在  其他
关注(0)|答案(1)|浏览(113)

我是Python的新手,有一个关于Cartopy能够在3D绘图中使用的问题。下面是一个使用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()

这将在三维轴内创建贴图,以便可以在曲面上打印对象。但是Cartopy返回一个matplotlib.axes.GeoAxesSubplot。不清楚如何将其与matplotlib-basemap一起添加到3D图形/轴上。
那么,有人可以给予任何指针,如何做一个类似的3D绘图与Cartopy?

wgmfuz8q

wgmfuz8q1#

底图mpl3d是一个非常巧妙的技巧,但它并没有按照所描述的方式来设计。因此,除了简单的海岸线之外,目前还不能将相同的技术用于其他地方。例如,填充的大陆只是不工作AFAICT。
也就是说,在使用cartopy时,类似的黑客也是可用的。由于我们可以一般地访问shapefile信息,因此此解决方案应该适用于任何折线shapefile,例如海岸线。
第一步是获取shapefile和相应的几何体:

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]

由于这些都是shapely的几何体,所以我们希望将它们转换为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示例,因此我们需要通过构建LineCollection来解决这个问题(就像底图一样)。遗憾的是,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()

使屈服:

相关问题