使用Python绘制球体上的温度分布

Question

我遇到了以下问题：

我在一个球体上有N个点，这些点的坐标保存在一个数组x中，数组的形状是(N,3)。这个数组包含了这些点的笛卡尔坐标。此外，在每个点上，我还有一个指定的温度，这个温度存储在一个数组T中，数组的形状是(N,)。

有没有什么简单的方法可以把这个温度分布用不同的颜色映射到平面上呢？

如果这样能简化任务，坐标也可以用极坐标(\theta,\phi)来表示。

Answer 1

一种实现这个效果的方法是通过颜色映射来设置 facecolors，也就是把你的热度数据映射到颜色上。

下面是一个例子：

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

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

u = np.linspace(0, 2 * np.pi, 80)
v = np.linspace(0, np.pi, 80)

# create the sphere surface
x=10 * np.outer(np.cos(u), np.sin(v))
y=10 * np.outer(np.sin(u), np.sin(v))
z=10 * np.outer(np.ones(np.size(u)), np.cos(v))

# simulate heat pattern (striped)
myheatmap = np.abs(np.sin(y))

ax.plot_surface(x, y, z, cstride=1, rstride=1, facecolors=cm.hot(myheatmap))

plt.show()

在这个例子中，我的“热图”只是沿着y轴的条纹，我是用函数 np.abs(np.sin(y)) 制作的，但任何从0到1的数值都可以使用（当然，这些数值的形状需要和 x 等的形状匹配）。

Answer 2

要绘制你的数据，你可以使用Basemap这个工具。唯一的问题是，contour和contourf这两个方法需要有网格化的数据。下面是一个在球面上使用简单（而且速度慢）的IDW样插值的例子。欢迎大家提出意见。

import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

def cart2sph(x, y, z):
    dxy = np.sqrt(x**2 + y**2)
    r = np.sqrt(dxy**2 + z**2)
    theta = np.arctan2(y, x)
    phi = np.arctan2(z, dxy)
    theta, phi = np.rad2deg([theta, phi])
    return theta % 360, phi, r

def sph2cart(theta, phi, r=1):
    theta, phi = np.deg2rad([theta, phi])
    z = r * np.sin(phi)
    rcosphi = r * np.cos(phi)
    x = rcosphi * np.cos(theta)
    y = rcosphi * np.sin(theta)
    return x, y, z

# random data
pts = 1 - 2 * np.random.rand(500, 3)
l = np.sqrt(np.sum(pts**2, axis=1))
pts = pts / l[:, np.newaxis]
T = 150 * np.random.rand(500)

# naive IDW-like interpolation on regular grid
theta, phi, r = cart2sph(*pts.T)
nrows, ncols = (90,180)
lon, lat = np.meshgrid(np.linspace(0,360,ncols), np.linspace(-90,90,nrows))
xg,yg,zg = sph2cart(lon,lat)
Ti = np.zeros_like(lon)
for r in range(nrows):
    for c in range(ncols):
        v = np.array([xg[r,c], yg[r,c], zg[r,c]])
        angs = np.arccos(np.dot(pts, v))
        idx = np.where(angs == 0)[0]
        if idx.any():
            Ti[r,c] = T[idx[0]]
        else:
            idw = 1 / angs**2 / sum(1 / angs**2)
            Ti[r,c] = np.sum(T * idw)

# set up map projection
map = Basemap(projection='ortho', lat_0=45, lon_0=15)
# draw lat/lon grid lines every 30 degrees.
map.drawmeridians(np.arange(0, 360, 30))
map.drawparallels(np.arange(-90, 90, 30))
# compute native map projection coordinates of lat/lon grid.
x, y = map(lon, lat)
# contour data over the map.
cs = map.contourf(x, y, Ti, 15)
plt.title('Contours of T')
plt.show()