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

Question

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

以下是我的数据图像。

我尝试过使用 scipy griddata 和 numpy 中的几种不同的插值方法，但它们都不能产生准确的结果。我相信获得准确结果的一种潜在方法是进行球形插值来创建高分辨率球形网格，然后使用 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();