3D Cartopy 中的轮廓

Question

我正在寻求将（可变）数量的填充轮廓绘制到 3D 图上的帮助。问题是这些点需要正确的地理参考。我已经使用 Cartopy 完成了 2D 案例，但不能简单地使用 mpl_toolkits.mplot3d，因为只能将一个投影传递给 figure() 方法。

This question 很有用，但主要专注于绘制 shapefile，而我拥有所有点和每个点的值以用于轮廓绘制。

This question 看起来也很有前途，但不处理 3D 轴。

我有一个使用直接mpl_toolkits.mplot3d 的方法，但它会扭曲数据，因为它位于错误的 CRS 中。我会使用Basemap，但由于某种原因它不能很好地处理 UTM 投影。

虽然看起来像这样（情节最终不那么明显，数据形成线性特征，但这应该有助于了解它是如何工作的）：

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

the_data = {'grdx': range(0, 100),
            'grdy': range(0, 100),
            'grdz': [[np.random.rand(100) for ii in range(100)]
                       for jj in range(100)]}
data_heights = range(0, 300, 50)

fig = plt.figure(figsize=(17, 17))
ax = fig.add_subplot(111, projection='3d')
x = the_data['grdx']
y = the_data['grdy']
ii = 0
for height in data_heights:
    print(height)
    z = the_data['grdz'][ii]
    shape = np.shape(z)
    print(shape)
    flat = np.ravel(z)
    flat[np.isclose(flat, 0.5, 0.2)] = height
    flat[~(flat == height)] = np.nan
    z = np.reshape(flat, shape)
    print(z)
    ax.contourf(y, x, z, alpha=.35)
    ii += 1
plt.show()

那么，我怎样才能使 contourf() 的 x 和 y 值成为 cartopy 可以在 3D 中处理的东西？

Answer 1

注意事项：

1. 每当我与主要维护者（Ben Root，GitHub 上的@weathergod）交谈时，matplotlib 中的 3d 内容通常被称为 2.5d。这应该让您知道它在 3d 中真正渲染的能力存在一些问题，而且 matplotlib 似乎不太可能解决其中一些问题（例如具有非常量 z 顺序的艺术家）。当渲染工作时，它非常棒。如果不这样做，则无能为力。
2. Cartopy 和 Basemap 都有一些技巧，可让您在 matplotlib 中使用 3d 模式进行可视化。他们真的是黑客 - YMMV，我想这不太可能进入核心 Basemap 或 Cartopy。

除此之外，我从Cartopy + Matplotlib (contourf) - Map Overriding data 那里得到了我的答案，你引用并从那里建立起来。

由于您想在轮廓之上构建，我采用了具有两个 Axes 实例（和两个图形）的方法。第一个是原始 2d（cartopy）GeoAxes，第二个是非 cartopy 3D 轴。在我执行 plt.show（或 savefig）之前，我只需关闭 2d GeoAxes（使用 plt.close(ax)）。

接下来，我利用 plt.contourf 的返回值是艺术家的集合这一事实，我们可以从中获取轮廓的坐标和属性（包括颜色）。

使用由 2d GeoAxes 和轮廓集合中的 contourf 生成的 2d 坐标，我将 z 维度（轮廓级别）插入到坐标中并构造一个 Poly3DCollection。

结果是这样的：

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import numpy as np


def f(x,y):
    x, y = np.meshgrid(x, y)
    return (1 - x / 2 + x**5 + y**3 + x*y**2) * np.exp(-x**2 -y**2)

nx, ny = 256, 512
X = np.linspace(-180, 10, nx)
Y = np.linspace(-90, 90, ny)
Z = f(np.linspace(-3, 3, nx), np.linspace(-3, 3, ny))


fig = plt.figure()
ax3d = fig.add_axes([0, 0, 1, 1], projection='3d')

# Make an axes that we can use for mapping the data in 2d.
proj_ax = plt.figure().add_axes([0, 0, 1, 1], projection=ccrs.Mercator())
cs = proj_ax.contourf(X, Y, Z, transform=ccrs.PlateCarree(), alpha=0.4)


for zlev, collection in zip(cs.levels, cs.collections):
    paths = collection.get_paths()
    # Figure out the matplotlib transform to take us from the X, Y coordinates
    # to the projection coordinates.
    trans_to_proj = collection.get_transform() - proj_ax.transData

    paths = [trans_to_proj.transform_path(path) for path in paths]
    verts3d = [np.concatenate([path.vertices,
                               np.tile(zlev, [path.vertices.shape[0], 1])],
                              axis=1)
               for path in paths]
    codes = [path.codes for path in paths]
    pc = Poly3DCollection([])
    pc.set_verts_and_codes(verts3d, codes)

    # Copy all of the parameters from the contour (like colors) manually.
    # Ideally we would use update_from, but that also copies things like
    # the transform, and messes up the 3d plot.
    pc.set_facecolor(collection.get_facecolor())
    pc.set_edgecolor(collection.get_edgecolor())
    pc.set_alpha(collection.get_alpha())

    ax3d.add_collection3d(pc)

proj_ax.autoscale_view()

ax3d.set_xlim(*proj_ax.get_xlim())
ax3d.set_ylim(*proj_ax.get_ylim())
ax3d.set_zlim(Z.min(), Z.max())


plt.close(proj_ax.figure)
plt.show()

当然，我们可以在这里进行大量分解，以及引入您所指的地理参考组件（例如拥有海岸线等）。

请注意，尽管输入坐标是纬度/经度，但 3d 轴的坐标是墨卡托坐标系的坐标 - 这告诉我们，关于我们要进行 cartopy 的变换，我们处于正确的轨道上对我们来说。

接下来，我从您引用的答案中提取包含土地多边形的代码。 matplotlib 3d 轴目前无法裁剪超出 x/y 限制的多边形，因此我需要手动进行。

把它们放在一起：

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

import itertools
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from matplotlib.collections import PolyCollection
import numpy as np


def f(x,y):
    x, y = np.meshgrid(x, y)
    return (1 - x / 2 + x**5 + y**3 + x*y**2) * np.exp(-x**2 -y**2)

nx, ny = 256, 512
X = np.linspace(-180, 10, nx)
Y = np.linspace(-90, 90, ny)
Z = f(np.linspace(-3, 3, nx), np.linspace(-3, 3, ny))


fig = plt.figure()
ax3d = fig.add_axes([0, 0, 1, 1], projection='3d')

# Make an axes that we can use for mapping the data in 2d.
proj_ax = plt.figure().add_axes([0, 0, 1, 1], projection=ccrs.Mercator())
cs = proj_ax.contourf(X, Y, Z, transform=ccrs.PlateCarree(), alpha=0.4)


for zlev, collection in zip(cs.levels, cs.collections):
    paths = collection.get_paths()
    # Figure out the matplotlib transform to take us from the X, Y coordinates
    # to the projection coordinates.
    trans_to_proj = collection.get_transform() - proj_ax.transData

    paths = [trans_to_proj.transform_path(path) for path in paths]
    verts3d = [np.concatenate([path.vertices,
                               np.tile(zlev, [path.vertices.shape[0], 1])],
                              axis=1)
               for path in paths]
    codes = [path.codes for path in paths]
    pc = Poly3DCollection([])
    pc.set_verts_and_codes(verts3d, codes)

    # Copy all of the parameters from the contour (like colors) manually.
    # Ideally we would use update_from, but that also copies things like
    # the transform, and messes up the 3d plot.
    pc.set_facecolor(collection.get_facecolor())
    pc.set_edgecolor(collection.get_edgecolor())
    pc.set_alpha(collection.get_alpha())

    ax3d.add_collection3d(pc)

proj_ax.autoscale_view()

ax3d.set_xlim(*proj_ax.get_xlim())
ax3d.set_ylim(*proj_ax.get_ylim())
ax3d.set_zlim(Z.min(), Z.max())


# Now add coastlines.
concat = lambda iterable: list(itertools.chain.from_iterable(iterable))

target_projection = proj_ax.projection

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

# Use the convenience (private) method to get the extent as a shapely geometry.
boundary = proj_ax._get_extent_geom()

# Transform the geometries from PlateCarree into the desired projection.
geoms = [target_projection.project_geometry(geom, feature.crs)
         for geom in geoms]
# Clip the geometries based on the extent of the map (because mpl3d can't do it for us)
geoms = [boundary.intersection(geom) for geom in geoms]

# Convert the geometries to paths so we can use them in matplotlib.
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=True)
ax3d.add_collection3d(lc, zs=ax3d.get_zlim()[0])

plt.close(proj_ax.figure)
plt.show() 

稍微完善一下，并将一些概念抽象为函数，这非常有用：

import cartopy.crs as ccrs
import cartopy.feature
from cartopy.mpl.patch import geos_to_path
import itertools
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d
from matplotlib.collections import PolyCollection, LineCollection
import numpy as np


def add_contourf3d(ax, contour_set):
    proj_ax = contour_set.collections[0].axes
    for zlev, collection in zip(contour_set.levels, contour_set.collections):
        paths = collection.get_paths()
        # Figure out the matplotlib transform to take us from the X, Y
        # coordinates to the projection coordinates.
        trans_to_proj = collection.get_transform() - proj_ax.transData

        paths = [trans_to_proj.transform_path(path) for path in paths]
        verts = [path.vertices for path in paths]
        codes = [path.codes for path in paths]
        pc = PolyCollection([])
        pc.set_verts_and_codes(verts, codes)

        # Copy all of the parameters from the contour (like colors) manually.
        # Ideally we would use update_from, but that also copies things like
        # the transform, and messes up the 3d plot.
        pc.set_facecolor(collection.get_facecolor())
        pc.set_edgecolor(collection.get_edgecolor())
        pc.set_alpha(collection.get_alpha())

        ax3d.add_collection3d(pc, zs=zlev)

    # Update the limit of the 3d axes based on the limit of the axes that
    # produced the contour.
    proj_ax.autoscale_view()

    ax3d.set_xlim(*proj_ax.get_xlim())
    ax3d.set_ylim(*proj_ax.get_ylim())
    ax3d.set_zlim(Z.min(), Z.max())

def add_feature3d(ax, feature, clip_geom=None, zs=None):
    """
    Add the given feature to the given axes.
    """
    concat = lambda iterable: list(itertools.chain.from_iterable(iterable))

    target_projection = ax.projection
    geoms = list(feature.geometries())

    if target_projection != feature.crs:
        # Transform the geometries from the feature's CRS into the
        # desired projection.
        geoms = [target_projection.project_geometry(geom, feature.crs)
                 for geom in geoms]

    if clip_geom:
        # Clip the geometries based on the extent of the map (because mpl3d
        # can't do it for us)
        geoms = [geom.intersection(clip_geom) for geom in geoms]

    # Convert the geometries to paths so we can use them in matplotlib.
    paths = concat(geos_to_path(geom) for geom in geoms)

    # Bug: mpl3d can't handle edgecolor='face'
    kwargs = feature.kwargs
    if kwargs.get('edgecolor') == 'face':
        kwargs['edgecolor'] = kwargs['facecolor']

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

    if kwargs.get('facecolor', 'none') == 'none':
        lc = LineCollection(polys, **kwargs)
    else:
        lc = PolyCollection(polys, closed=False, **kwargs)
    ax3d.add_collection3d(lc, zs=zs)

我用来制作以下有趣的 3D 罗宾逊情节：

def f(x, y):
    x, y = np.meshgrid(x, y)
    return (1 - x / 2 + x**5 + y**3 + x*y**2) * np.exp(-x**2 -y**2)


nx, ny = 256, 512
X = np.linspace(-180, 10, nx)
Y = np.linspace(-89, 89, ny)
Z = f(np.linspace(-3, 3, nx), np.linspace(-3, 3, ny))


fig = plt.figure()
ax3d = fig.add_axes([0, 0, 1, 1], projection='3d')

# Make an axes that we can use for mapping the data in 2d.
proj_ax = plt.figure().add_axes([0, 0, 1, 1], projection=ccrs.Robinson())
cs = proj_ax.contourf(X, Y, Z, transform=ccrs.PlateCarree(), alpha=1)

ax3d.projection = proj_ax.projection
add_contourf3d(ax3d, cs)

# Use the convenience (private) method to get the extent as a shapely geometry.
clip_geom = proj_ax._get_extent_geom().buffer(0)


zbase = ax3d.get_zlim()[0]
add_feature3d(ax3d, cartopy.feature.OCEAN, clip_geom, zs=zbase)
add_feature3d(ax3d, cartopy.feature.LAND, clip_geom, zs=zbase)
add_feature3d(ax3d, cartopy.feature.COASTLINE, clip_geom, zs=zbase)

# Put the outline (neatline) of the projection on.
outline = cartopy.feature.ShapelyFeature(
    [proj_ax.projection.boundary], proj_ax.projection,
    edgecolor='black', facecolor='none')
add_feature3d(ax3d, outline, clip_geom, zs=zbase)


# Close the intermediate (2d) figure
plt.close(proj_ax.figure)
plt.show()

回答这个问题很有趣，让我想起了一些 matplotlib 和 cartopy 转换的内部结构。毫无疑问，它有能力产生一些有用的可视化，但由于 matplotlib 的 3d (2.5d) 实现固有的问题，我个人不会在生产中使用它。

HTH

Answer 2

在我的环境中，错误 无法执行“GEOSIntersection_r”。可能的原因是几何 无效 只需删除导致错误的那些即可解决

geoms2 = []
for i in range(len(geoms)) :
    if geoms[i].is_valid :
        geoms2.append(geoms[i])
geoms = geoms2 

在路口之前。到目前为止，结果对我来说看起来不错。