3D Cartopy中的轮廓线 [英] contourf in 3D Cartopy

问题描述

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

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

请注意，尽管输入坐标是经/纬度，但3d轴的坐标是墨卡托坐标系的坐标-这告诉我们，在我们要进行Cartopy的变换方面，我们处在正确的轨道上为我们.

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

综合起来:

将 cartopy.crs 导入为 ccrs导入 cartopy.feature从cartopy.mpl.patch导入geos_to_path导入迭代工具导入matplotlib.pyplot作为plt从 mpl_toolkits.mplot3d.art3d 导入 Poly3DCollection从matplotlib.collections导入PolyCollection将numpy导入为npdef 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, 纽约 = 256, 512X = np.linspace(-180，10，nx)Y = np.linspace(-90, 90, ny)Z = f(np.linspace(-3, 3, nx), np.linspace(-3, 3, ny))无花果= plt.figure()ax3d = fig.add_axes([0, 0, 1, 1], 投影='3d')#制作一个可用于在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)对于 zlev，zip 中的集合(cs.levels，cs.collections):路径 = collection.get_paths()#找出matplotlib变换以从X，Y坐标处获取我们# 到投影坐标.trans_to_proj = collection.get_transform()-proj_ax.transData路径= [trans_to_proj.transform_path(path)用于路径中的路径]verts3d = [np.concatenate([path.vertices，np.tile(zlev，[path.vertices.shape [0]，1])]，轴=1)对于路径中的路径]代码= [路径中路径的path.codes]pc = Poly3DCollection([])pc.set_verts_and_codes(verts3d，代码)#手动复制轮廓(如颜色)中的所有参数.# 理想情况下，我们会使用 update_from，但这也会复制诸如#进行转换，并弄乱3d图.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())#现在添加海岸线.concat = lambda iterable: list(itertools.chain.from_iterable(iterable))target_projection = proj_ax.projection功能= cartopy.feature.NaturalEarthFeature('physical'，'land'，'110m')geoms = feature.geometries()# 使用便捷(私有)方法将范围作为匀称几何体获取.边界= proj_ax._get_extent_geom()# 将 PlateCarree 中的几何图形转换为所需的投影.geoms = [target_projection.project_geometry(geom，feature.crs)对于geom中的geom]#根据地图的范围裁剪几何形状(因为mpl3d不能为我们做这件事)geoms = [boundary.intersection(geom) for geom in geoms]# 将几何图形转换为路径，以便我们可以在 matplotlib 中使用它们.路径 = concat(geos_to_path(geom) for geom in geoms)polys = concat(path.to_polygons() 用于路径中的路径)lc = PolyCollection(polys，edgecolor ='black'，facecolor='绿色'，关闭=真)ax3d.add_collection3d(lc，zs = ax3d.get_zlim()[0])plt.close(proj_ax.figure)plt.show()

对此进行一点处理，并将一些概念抽象为函数，这将非常有用:

将 cartopy.crs 导入为 ccrs导入 cartopy.feature从 cartopy.mpl.patch 导入 geos_to_path导入迭代工具导入matplotlib.pyplot作为plt导入 mpl_toolkits.mplot3d从matplotlib.collections导入PolyCollection，LineCollection将numpy导入为npdef add_contourf3d(ax, contour_set):proj_ax = contour_set.collections[0].axes对于zlev，以zip格式收集(contour_set.levels，contour_set.collections):路径 = collection.get_paths()#找出matplotlib转换以将我们从X，Y中引出#坐标为投影坐标.trans_to_proj = collection.get_transform() - proj_ax.transData路径= [trans_to_proj.transform_path(path)用于路径中的路径]verts = [path.vertices路径中的路径]代码= [路径中路径的path.codes]pc = PolyCollection([])pc.set_verts_and_codes(verts, 代码)# 手动复制轮廓中的所有参数(如颜色).# 理想情况下，我们会使用 update_from，但这也会复制诸如# 变换，弄乱了 3d 图.pc.set_facecolor(collection.get_facecolor())pc.set_edgecolor(collection.get_edgecolor())pc.set_alpha(collection.get_alpha())ax3d.add_collection3d(pc，zs = zlev)# 根据轴的限制更新 3d 轴的限制#产生轮廓.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):"""将给定特征添加到给定轴."""concat = lambda iterable: list(itertools.chain.from_iterable(iterable))target_projection = ax.projectiongeoms = list(feature.geometries())如果target_projection！= feature.crs:#将几何图形从要素的CRS转换为#所需的投影.geoms = [target_projection.project_geometry(geom, feature.crs)对于 geoms 中的 geom]如果clip_geom:# 根据地图的范围裁剪几何图形(因为 mpl3d#无法为我们做到)geoms = [geom中的geom的geom.intersection(clip_geom)]# 将几何图形转换为路径，以便我们可以在 matplotlib 中使用它们.path = concat(geom中geom的geos_to_path(geom))# 错误:mpl3d 无法处理 edgecolor='face'kwargs = feature.kwargs如果 kwargs.get('edgecolor') == 'face':kwargs['edgecolor'] = kwargs['facecolor']polys = concat(path.to_polygons(closed_only=False) 用于路径中的路径)如果 kwargs.get('facecolor', 'none') == 'none':lc = LineCollection(polys, **kwargs)别的:lc = PolyCollection(polys，closed = False，** kwargs)ax3d.add_collection3d(lc，zs = zs)

我曾经制作过以下有趣的3D Robinson图:

  def f(x，y):x, y = np.meshgrid(x, y)返回 (1 - x/2 + x**5 + y**3 + x*y**2) * np.exp(-x**2 -y**2)nx, 纽约 = 256, 512X = np.linspace(-180，10，nx)Y = np.linspace(-89，89，ny)Z = f(np.linspace(-3，3，nx)，np.linspace(-3，3，ny))无花果= plt.figure()ax3d = fig.add_axes([0, 0, 1, 1], 投影='3d')#制作一个可用于在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.projectionadd_contourf3d(ax3d，cs)# 使用便捷(私有)方法将范围作为匀称几何体获取.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)# 放置投影的轮廓(neatline).轮廓= cartopy.feature.ShapelyFeature([proj_ax.projection.boundary]，proj_ax.projection，edgecolor='black', facecolor='none')add_feature3d(ax3d，轮廓，clip_geom，zs=zbase)#关闭中间图(2d)plt.close(proj_ax.figure)plt.show()

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

HTH

I am looking for help in plotting a (variable) number of filled contours onto a 3D plot. The rub is that the points need to be correctly geo-referenced. I have got the 2D case working, using Cartopy, but one can not simply use mpl_toolkits.mplot3d, since one can only pass one projection into the figure() method.

This question was useful, but is focused mostly on plotting a shapefile, while I have all the points and the values at each point for use in the contouring.

This question also looked promising, but does not deal with a 3D axis.

I have a method working using straight mpl_toolkits.mplot3d, but it is distorting the data, since it is in the wrong CRS. I would use Basemap, but it does not handle UTM projections very well for some reason.

It looks something like this though (the plot ends up much less spotted, the data forms linear features, but this should serve to give an idea of how it works):

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

So how can I make the x and y values for the contourf() something that cartopy can handle in 3D?

解决方案

Caveats:

1. The 3d stuff in matplotlib is frequently referred to as 2.5d whenever I speak to the primary maintainer (Ben Root, @weathergod on GitHub). This should give you an indication that there are a few issues with its ability to truly render in 3d, and it seems unlikely that matplotlib will ever be able to address some of these issues (like artists having a non-constant z order). When the rendering works, it is pretty awesome. When it doesn't, there isn't a lot that can be done about it.
2. Cartopy and Basemap both have hacks that allow you to visualise with the 3d mode in matplotlib. They really are hacks - YMMV, and I imagine this isn't something that is likely to go into core Basemap or Cartopy.

With that out of the way, I took my answer from Cartopy + Matplotlib (contourf) - Map Overriding data that you referenced and built-up from there.

Since you want to build on top of contours, I took the approach of having two Axes instances (and two figures). The first is the primitive 2d (cartopy) GeoAxes, the second is the non-cartopy 3D axes. Right before I do a plt.show (or savefig), I simply close the 2d GeoAxes (with plt.close(ax)).

Next, I use the fact that the return value from a plt.contourf is a collection of artists, from which we can take the coordinates and properties (including color) of the contours.

Using the 2d coordinates that are generated by the contourf in the 2d GeoAxes and contour collection, I insert the z dimension (the contour level) into the coordinates and construct a Poly3DCollection.

This turns out as something like:

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

Of course, there is a bunch of factorisation we can do here, as well as bringing in the georeferenced component you were referring to (like having coastlines etc.).

Notice that despite the input coordinates being lat/longs, the coordinates of the 3d axes are those of a Mercator coordinate system - this tells us that we are on the right track with regards to the transforms that we are getting cartopy to do for us.

Next, I take the code from the answer you referenced to include land polygons. The matplotlib 3d axes currently has no ability to clip polygons that fall outside of the x/y limits, so I needed to do that manually.

Bringing it all together:

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

Rounding this off a bit, and abstracting a few of the concepts to functions makes this pretty useful:

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)

Which I used to produce the following fun 3D Robinson plot:

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

Answering this question was a lot of fun, and reminded me of some of the matplotlib & cartopy transform internals. There is no doubt that it has the power to produce some useful visualisations, but I personally wouldn't be using it in production due to the issues inherent with matplotlib's 3d (2.5d) implementation.

HTH

这篇关于3D Cartopy中的轮廓线的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！