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

Question

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

的形式

以下是我的数据图片。

我尝试在 numpy 中使用 scipy griddata 和几种不同的插值方法，但其中 none 产生了准确的结果。我相信获得准确结果的一种潜在方法是进行球面插值以创建高分辨率球面网格，然后使用 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 中进行适当的球面插值？抱歉，如果我遗漏了什么，这是我第一次在 Whosebug 上发帖！

Answer 1

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

结构丢失

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

matplotlib 中的函数 plt.tripcolor 直接为您完成此操作：阴影选项可以设置为 'gouraud' 以获得平滑的渲染。我将它设置为 'flat' 以查看从网格化中获得的三角形。

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: 

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

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

Interpolating spherical grid to regular grid?

python

interpolation

numpy

scipy

data-science