John
John

Reputation: 401

How to convert contourf values and plot to surface of 3d sphere

I am trying to use matplotlib contourf to generate a contour plot of temperature data. I would then like to map this data onto a 3d sphere. I am using vpython to render the image. The issue I have is that the polygons are not generated across the surface of the sphere, and also there are many gaps in the data. Can someone please explain how this can be achieved ?

Please note I do not wish to solve this by generating a texture from the contour and then wrapping the sphere with the texture. I wish to translate the polygon paths generated by contourf, and directly translate those paths onto the sphere.

from vpython import vector

import numpy as np
from matplotlib import pyplot as plt

from vpython import triangle, vertex, vec


from scipy.interpolate import griddata


def spherical_to_cartesian(lat, lon, radius=3):
    lons = np.radians(lon)
    lats = np.radians(lat)


    x = radius * np.cos(lats) * np.cos(lons)
    y = radius * np.cos(lats) * np.sin(lons)
    z = radius * np.sin(lats)
    return np.array([x, y, z])


shape = (721, 1440)


lats = np.linspace(-90, 90, shape[0])
lons = np.linspace(-180, 180, shape[1])


min_temp = -30
max_temp = 50
temps = np.random.uniform(min_temp, max_temp, size=shape)



lons_grid, lats_grid = np.meshgrid(lons, lats)

new_lons = np.linspace(lons.min(), lons.max(), 72)  # 72 points in longitude
new_lats = np.linspace(lats.min(), lats.max(), 36)  # 36 points in latitude
new_lons_grid, new_lats_grid = np.meshgrid(new_lons, new_lats)

radius = 3

lats_flat = lats_grid.flatten()
lons_flat = lons_grid.flatten()
temps_flat = temps.flatten()

coarse_temps = griddata(
    (lons_flat, lats_flat),
    temps_flat,
    (new_lons_grid, new_lats_grid),
    method='linear'  # Use 'linear', 'nearest', or 'cubic' as needed
)


norm = plt.Normalize(coarse_temps.min(), vmax=coarse_temps.max())

cmap = plt.get_cmap("inferno", 100)

levels = 100

levels = np.linspace(coarse_temps.min(), coarse_temps.max(), levels + 1)

bucketed_data = np.digitize(coarse_temps, levels) - 1

fig, ax = plt.subplots()
contours = ax.contourf(new_lons_grid, new_lats_grid, bucketed_data, levels=np.arange(len(levels)))
plt.close(fig)



def create_polygon(region, temp_color):
    if len(region) < 3:
        # Can't form a polygon with fewer than 3 points
        return

    # Calculate the centroid of the region
    centroid = vec(sum(p[0] for p in region) / len(region),
                   sum(p[1] for p in region) / len(region),
                   sum(p[2] for p in region) / len(region))

    # Tessellate the region into triangles
    for i in range(len(region) - 1):
        v0 = vec(*region[i])
        v1 = vec(*region[i + 1])
        triangle(
            v0=vertex(pos=v0, color=temp_color),
            v1=vertex(pos=v1, color=temp_color),
            v2=vertex(pos=centroid, color=temp_color)
        )


    v0 = vec(*region[-1])
    v1 = vec(*region[0])
    triangle(
        v0=vertex(pos=v0, color=temp_color),
        v1=vertex(pos=v1, color=temp_color),
        v2=vertex(pos=centroid, color=temp_color)
    )


def split_contours(segs, kinds=None):
    if kinds is None:
        return segs  # nothing to be done
    new_segs = []
    for i, seg in enumerate(segs):
        segkinds = kinds[i]
        boundaries = [0] + list(np.nonzero(segkinds == 79)[0])
        for b in range(len(boundaries) - 1):
            new_segs.append(seg[boundaries[b] + (1 if b > 0 else 0):boundaries[b + 1]])
    return new_segs



allsegs = contours.allsegs
allkinds = contours.allkinds
colors = cmap(norm(coarse_temps))
rgb_colors = [
    tuple(int(c * 255) for c in color[:3])
    for color in colors.reshape(-1, colors.shape[-1])
]
for clev in range(len(contours.allsegs)):
    kinds = None if allkinds is None else allkinds[clev]
    segs = split_contours(allsegs[clev], kinds)
    rgb = cmap(clev)[:3]
    coords = [spherical_to_cartesian(lat, lon, radius=3) for seg in segs for lon, lat in seg]
    temp_color = vector(*rgb)
    create_polygon(coords, temp_color)



import time
while True:
    time.sleep(0.03)

This is how the sphere is rendered. I would like the contours to be rendered on the surface of the sphere, as if the 2d contours were projected onto the 3d surface.

enter image description here

Upvotes: 2

Views: 190

Answers (0)

Related Questions