3D Cartopy 中的轮廓
mtb*_*-za 3 matplotlib cartopy
我正在寻找在 3D 图上绘制（可变）数量的填充轮廓的帮助。问题是这些点需要正确地进行地理参考。我已经使用 Cartopy 处理了 2D 案例，但不能简单地使用mpl_toolkits.mplot3d，因为只能将一个投影传递到figure()方法中。
这个问题很有用，但主要集中在绘制 shapefile，而我拥有所有点和每个点的值以用于轮廓绘制。
这个问题看起来也很有希望，但不涉及 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()cartopy 可以在 3D 中处理的东西制作 x 和 y 值呢？
注意事项：
- 每当我与主要维护者（Ben Root，GitHub 上的 @weathergod）交谈时，matplotlib 中的 3d 内容经常被称为 2.5d。这应该表明它在 3d 中真正渲染的能力存在一些问题，而且 matplotlib 似乎不太可能解决其中的一些问题（例如具有非常量 z 顺序的艺术家）。当渲染工作时，它非常棒。当它没有时，没有太多可以做的。
- Cartopy 和 Basemap 都有一些技巧，可以让您在 matplotlib 中使用 3d 模式进行可视化。他们真的是黑客 - YMMV，我想这不太可能进入核心 Basemap 或 Cartopy。
顺便说一下，我从Cartopy + Matplotlib (contourf) -你从那里引用和构建的地图覆盖数据中得到了我的答案。
由于您想在轮廓之上构建，我采用了具有两个 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 Robinson 图：
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
| 归档时间： |
|
| 查看次数： |
2315 次 |
| 最近记录： |