Plotting the temperature distribution on a sphere with python

3
votes

I have the following problem:

a have N points on a sphere specified by a array x, with x.shape=(N,3). This array contains their cartesian coordinates. Furthermore, at each point, I have a specified temperature. This quantity is saved in an array T, with T.shape=(N,).

Is there any straight forward way to map this temperature distribution into the plane using different colors?

If it simplifies the task, the position can also be given in polar coordinates (\theta,\phi).

2
pythonmatplotlibplot

2 Answers

10
votes

To plot your data, you can use Basemap. The only problem is, that both contour and contourf routines needs gridded data. Here is example with naive (and slow) IDW-like interpolation on sphere. Any comments are welcome.

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

enter image description here

5
votes

One way to do this is to set facecolors by mapping your heat data through the colormap.

Here's an example:

enter image description here

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

Here, my "heatmap" is just stripes along the y-axis, which I made using the function np.abs(np.sin(y)), but anything that goes form 0 to 1 will work (and, of course, it needs to match the shapes on x, etc.