Interpolating spherical grid to regular grid?

Question

I am trying to properly interpolate values from a circular grid onto a regular grid with Python3. The data points are sparse compared to my grid goal of 400x400. My goal is to be able to take these values and display them accurately onto an image of the earth. My input data is in the form of [x, y, value].

The following is an image of my data.

I have tried using scipy griddata and several different interpolation methods in numpy but none of them produce accurate results. I believe a potential way to get accurate results is to do spherical interpolation to create a high res spherical grid, then use griddata to map it to a rectangular grid, but I have no idea as to using spherical interpolation for this. Following is a couple images, ignore the orientation of the photos as they are from different times.

Using numpy interp2d, I get this:

What I would like to get is something similar to this, where it is very smooth as it should be:

Here is code to reproduce the problem. Only numpy, matplotlib, and scipy are required. The get_rotation_array() function with no args, creates a pretty close example of what the sample data could be, for anyone testing.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from scipy import interpolate

# GLOBALS
EARTH_RADIUS = 6370997.0
SOLAR_GRID_RES_KM = 750000
EARTH_GRID_RES_KM = 5*100000
CUT_OFF_VAL = 1000000

# Earth Patches
earth_circle1 = plt.Circle((EARTH_RADIUS, EARTH_RADIUS), EARTH_RADIUS, edgecolor='black', fill=False, linewidth=1)
earth_circle2 = plt.Circle((EARTH_RADIUS, EARTH_RADIUS), EARTH_RADIUS, edgecolor='black', fill=False, linewidth=1)

# This function is messy but it roughly simulates
# what kind of data I am expecting
def get_rotation_array(steps=20, per_line=9):
    x_vals = []
    y_vals = []
    z_vals = []
    r = EARTH_RADIUS - 10000
    for el in range(1, per_line):
        for t in np.linspace(0, 2*np.pi, num=steps):
            x = (el/float(per_line - 1))*r*np.cos(t) + EARTH_RADIUS
            y = (el/float(per_line - 1))*r*np.sin(t) + EARTH_RADIUS
            z = el - 2*(el/float(per_line - 1))*np.abs((1.5*np.pi) - t)
            if y < (EARTH_RADIUS + CUT_OFF_VAL):
                x_vals.append(x)
                y_vals.append(y)
                z_vals.append(z)

    x_vals.append(EARTH_RADIUS)
    y_vals.append(EARTH_RADIUS)
    z_vals.append(1)

    return np.array(x_vals), np.array(y_vals), np.array(z_vals)

# Get "Sample" Data
x, y, z = get_rotation_array()

# Create Sublots
fig, ax = plt.subplots(1, 2)

# Get Values for raw plot
cmap = cm.get_cmap("jet", 1000)
alpha = np.interp(z, [z.min(), z.max()], [0, 1])
colour = cmap(alpha)

# Plot Raw Plot
ax[0].set_title("Sample Data")
ax[0].scatter(x, y, c=colour)
ax[0].add_patch(earth_circle1)
ax[0].set_xlim([0,EARTH_RADIUS*2])
ax[0].set_ylim([0,EARTH_RADIUS*2])

# Use griddata interpolation
x_solar_interp = np.arange(0, EARTH_RADIUS*2, EARTH_GRID_RES_KM)
y_solar_interp = np.arange(0, EARTH_RADIUS + CUT_OFF_VAL, EARTH_GRID_RES_KM)
xx_interp, yy_interp = np.meshgrid(x_solar_interp, y_solar_interp)

z_interp = interpolate.griddata((x, y), z, (xx_interp, yy_interp), method='linear')

# Plot the Colormesh
plt.pcolormesh(xx_interp, yy_interp, z_interp, cmap=cmap, shading='flat')

# Plot Interpolated Data
ax[1].set_title("Interpolated")
ax[1].add_patch(earth_circle2)
ax[1].set_xlim([0,EARTH_RADIUS*2])
ax[1].set_ylim([0,EARTH_RADIUS*2])

# Show the plots
plt.show()

The interpolation breaks down because the data isn't dependent on the x,y values and is dependent on the angle from the center of the earth. So after all that it comes down to, how to do proper spherical interpolation in Python3 with data like this? Sorry if I have missed anything this is my first time posting on StackOverflow!

Answer 1

There are different way to do this. I think the main point is the distinction between unstructured data (i.e. only coordinates of points are given, not the mesh) and structured data (i.e. points are on a grid). In your case, the data is initially structured (points obtained using meshgrid), but the structure is lost by using the loop to compute z.

To plot a surface using unstructured data (and for interpolation), a mesh have to be first computed (using Delaunay triangulation).

The function plt.tripcolor from matplotlib do this for you directly: The shading option can be set to 'gouraud' to obtain a smooth rendering. I set it to 'flat' to see the triangles obtained from the meshing.

plt.figure(figsize=(8,8))
ax = plt.subplot(aspect='equal')
cmap = cm.get_cmap('jet')


plt.tripcolor(x, y, z, cmap=cmap, shading='flat'); # use shading='gouraud' to smooth
ax.plot(x, y, '.', color='red', label='data points');

earth_circle = plt.Circle((EARTH_RADIUS, EARTH_RADIUS), EARTH_RADIUS,
                           edgecolor='black', fill=False, linewidth=1);

ax.add_artist(earth_circle);
ax.set_xlabel('x (m)'); ax.set_ylabel('y (m)');
cbar = plt.colorbar();
cbar.set_label('z')
ax.legend();

If the data are still wanted on a cartesian grid, they can be interpolated using griddata. The interpolation is based on a similar Delaunay triangulation. Then, the function pcolormesh can be used to plot the surface:

# Get Values for griddata plot
# Use griddata interpolation
EARTH_GRID_RES_KM = 5*100000 # changed! to emphasis what is really plotted
x_solar_interp = np.arange(0, EARTH_RADIUS + CUT_OFF_VAL, EARTH_GRID_RES_KM)
y_solar_interp = np.arange(0, EARTH_RADIUS*2, EARTH_GRID_RES_KM)
xx_interp, yy_interp = np.meshgrid(x_solar_interp, y_solar_interp)

z_interp = interpolate.griddata((x, y), z, (xx_interp, yy_interp),
                                method='linear', fill_value=np.nan)

# Graph
plt.figure(figsize=(8,8))
ax = plt.subplot(aspect='equal')

cmap = cm.get_cmap('jet')
cmap.set_bad(color='white')

plt.pcolormesh(xx_interp, yy_interp, z_interp, cmap=cmap,
               shading='flat'); # try shading='gouraud'
# note about pcolormesh dealing with NaN: https://stackoverflow.com/a/33667896/8069403

earth_circle = plt.Circle((EARTH_RADIUS, EARTH_RADIUS), EARTH_RADIUS,
                           edgecolor='black', fill=False, linewidth=1);
ax.add_artist(earth_circle);

ax.plot(xx_interp.flatten(), yy_interp.flatten(), '.',
        color='black', label='data points');

ax.set_xlabel('x (m)'); ax.set_ylabel('y (m)');
cbar = plt.colorbar(cmap=cmap);
cbar.set_label('z')
ax.legend();