将球形网格插值到常规网格？

Question

我正在尝试使用 Python3 将圆形网格中的值正确插入到常规网格中。与我的 400x400 网格目标相比，数据点稀疏。我的目标是能够获取这些值并将它们准确地显示在地球图像上。我的输入数据是 [x, y, value] 的形式。

以下是我的数据图片。

我已经尝试在 numpy 中使用 scipy griddata 和几种不同的插值方法，但它们都没有产生准确的结果。我相信获得准确结果的一种潜在方法是进行球面插值以创建高分辨率球面网格，然后使用griddata 将其映射到矩形网格，但我不知道为此使用球面插值。以下是几张图片，请忽略照片的方向，因为它们来自不同的时间。

使用 numpy interp2d，我明白了：

我想得到的是与此类似的东西，它应该是非常平滑的：

这里是重现问题的代码。只需要 numpy、matplotlib 和 scipy。没有参数的 get_rotation_array() 函数为任何测试的人创建了一个非常接近示例数据的示例。

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

由于数据不依赖于 x,y 值，而是依赖于与地球中心的角度，因此插值失效。 因此，归根结底，如何在 Python3 中使用这样的数据进行适当的球面插值？抱歉，如果我遗漏了什么，这是我第一次在 StackOverflow 上发帖！

Answer 1

有不同的方法可以做到这一点。我认为主要观点是非结构化数据之间的区分（即，仅给出点的坐标，而不是网格）和结构化数据（即点在网格上）。在您的情况下，数据最初是结构化的（使用 meshgrid 获得的点），但使用循环计算 z 会丢失结构。

使用非结构化数据（以及用于插值），必须首先计算网格（使用Delaunay triangulation）。

matplotlib 中的函数plt.tripcolor 直接为您执行此操作： 可以将着色选项设置为“gouraud”以获得平滑的渲染。我将其设置为“平面”以查看从网格划分中获得的三角形。

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

如果仍然需要笛卡尔网格上的数据，可以使用griddata 对它们进行插值。插值基于类似的 Delaunay 三角剖分。然后，函数pcolormesh可以用来绘制曲面：

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