Smooth 2D Plot of Hemispherical Coordinates using matplotlib

Question

I have some n points on a hemisphere (theta in range (0, 90) and phi in range (0, 180)). I want to have a 2D plot of the heatmap since 3D plots have occlusion. In addition, since the n points are located at spaced interval, a smooth plot (say, a Gaussian smoothing) will probably look better.

My attempt: I found polar plot of matplotlib here which looks something like what I want although (a) the grid coordinates are mislabeled and (b) it is not smoothened for spaced points.

Edit: My minimum working example

import numpy as np
import matplotlib.pyplot as plt

def to_degrees(x):
    return x*np.pi/180.0

def get_projection(phi, lmda, phi_0=0.0, lmda_0=to_degrees(90.0)):
    # Credits : https://en.wikipedia.org/wiki/Orthographic_map_projection
    x = np.cos(phi)*np.sin(lmda - lmda_0)
    y = np.cos(phi_0)*np.sin(phi) - np.sin(phi_0)*np.cos(phi)*np.cos(lmda-lmda_0)
    return [x, y]

# Adding latitudes and longitudes to give the appearance of a sphere
latitudes = [60, 30, 0, -30, -60] #elevations
longitudes = [0, 30, 60, 90, 120, 150, 180] #azimuths

plt.gca().set_aspect('equal', adjustable='box')
for longitude in longitudes:
    prev_point = get_projection(to_degrees(-90.), to_degrees(0))
    for latitude in range(-90, 90):
        curr_point = get_projection(to_degrees(latitude), to_degrees(longitude))
        plt.plot([prev_point[0], curr_point[0]], [prev_point[1], curr_point[1]], 'k', alpha=0.3)
        prev_point = curr_point

for latitude in latitudes:
    prev_point = get_projection(to_degrees(latitude), to_degrees(0))
    for longitude in range(0, 180):
        curr_point = get_projection(to_degrees(latitude), to_degrees(longitude))
        plt.plot([prev_point[0], curr_point[0]], [prev_point[1], curr_point[1]], 'k', alpha=0.3)
        prev_point = curr_point

views = [[-60, 0], [60, 0]] # and similar points of the format [azimuth, elevation]
frequency = [0.5, 0.3] # and similar numbers in range [0,1] for heatmap

for view_idx in range(len(views)):
    loc = get_projection(to_degrees(views[view_idx][0]), to_degrees(views[view_idx][1]))
    plt.scatter(loc[0], loc[1], s=300, c=np.array(plt.cm.jet(frequency[view_idx])).reshape(1, -1))

plt.show()

to get this

Since I have 11-12 such points spread all over the hemisphere, I want to make the heatmap smooth as well.

Answer 1

Based on this post,

you could create a mesh, calculate the colours with a function then use imshow with interpolation

I wrote a function create_gaussian_mesh to solve the problem

def create_gaussian_mesh(views,cmap_names,t_x,t_y,radii,ax):
    """
    views: the points
    cmap_names: for heatmap
    t_x: the number of grids in x direction 
    t_y: the number of grids in y direction 
    radii: the radii of the Gaussians to plot
    ax: the canvas 
    """
    def gaussian(view,radius):
        
        # initialize a patch and grids  
        patch = np.empty((t_x,t_y))
        patch[:,:] = np.nan
        x = np.linspace(-1,1,t_x)
        y = np.linspace(-1,1,t_y)
        x_grid,y_grid = np.meshgrid(x, y)
    
        loc = get_projection(to_degrees(view[0]),to_degrees(view[1]))
        # threshold controls the size of the gaussian 
        circle_mask = (x_grid-loc[0])**2 + (y_grid-loc[1])**2 < radius
        gaussian_value = np.exp((x_grid-loc[0])**2+(y_grid-loc[1])**2)
        patch[circle_mask] = gaussian_value[circle_mask]
    
        return patch

    # modify the patch
    for view,cmap_name,radius in zip(views,cmap_names,radii): 
        patch = gaussian(view,radius)
        extent = -1,1,-1,1
        cmap = plt.get_cmap(cmap_name)
        ax.imshow(patch,cmap=cmap,alpha=.6,interpolation='bilinear',extent=extent)

Using this function, you can plot your data

import numpy as np
import matplotlib.pyplot as plt

to_degrees = lambda x: np.deg2rad(x)

def get_projection(phi, lmda, phi_0=0.0, lmda_0=to_degrees(90.0)):
    # Credits : https://en.wikipedia.org/wiki/Orthographic_map_projection
    x = np.cos(phi)*np.sin(lmda - lmda_0)
    y = np.cos(phi_0)*np.sin(phi) - np.sin(phi_0)*np.cos(phi)*np.cos(lmda-lmda_0)
    return x, y

# Adding latitudes and longitudes to give the appearance of a sphere
latitudes = [60, 30, 0, -30, -60] #elevations
longitudes = [0, 30, 60, 90, 120, 150, 180] #azimuths

fig,ax = plt.subplots(figsize=(8,8))
ax.set_aspect('equal', adjustable='box')
ax.set_xlim(-1,1)
ax.set_ylim(-1,1)
# set the right ticks 
x_ticks = ['$%i\degree$' % ind for ind in np.linspace(0,180,10).astype(int)]
y_ticks = ['$%i\degree$' % ind for ind in np.linspace(-90,90,10).astype(int)]
ax.set_xticks(np.linspace(-1,1,10));ax.set_xticklabels(x_ticks)
ax.set_yticks(np.linspace(-1,1,10));ax.set_yticklabels(y_ticks)

for longitude in longitudes:
    prev_point = get_projection(to_degrees(-90.), to_degrees(0))
    for latitude in range(-90, 90):
        curr_point = get_projection(to_degrees(latitude), to_degrees(longitude))
        ax.plot([prev_point[0], curr_point[0]], [prev_point[1], curr_point[1]], 'k', alpha=0.3)
        prev_point = curr_point
for latitude in latitudes:
    prev_point = get_projection(to_degrees(latitude), to_degrees(0))
    for longitude in range(0, 180):
        curr_point = get_projection(to_degrees(latitude), to_degrees(longitude))
        ax.plot([prev_point[0], curr_point[0]], [prev_point[1], curr_point[1]], 'k', alpha=0.3)
        prev_point = curr_point

views = [[-60, 0], [60, 0], [20,10],[30,45]] # and similar points of the format [azimuth, elevation]
# instead of frequencies, you need a list of names of cmaps 
cmap_names= ['gray','hot','cool','Greys']
    
# the radius of the gaussians to plot 
radii = np.linspace(0.08,0.1,len(views))
create_gaussian_mesh(views,cmap_names,t_x=300,t_y=300,radii=radii,ax=ax)

in which the ticks of both x and y axis are properly set. The output figure of the code is