如何制作不规则数据的网格？

Question

我有经度、纬度和数据的 numpy 数组。 我想使用 numpy、scipy 和 matplotlib 将此数据绘制为光栅图像。

import numpy as np
from matplotlib.mlab import griddata
import matplotlib.pyplot as plt

longitudes = np.array([[139.79391479492188, 140.51760864257812, 141.19119262695312, 141.82083129882812, 142.41165161132812],
        [139.79225158691406, 140.51416015625, 141.18606567382812, 141.8140869140625, 142.40338134765625],
        [139.78591918945312, 140.50637817382812, 141.17694091796875, 141.80377197265625, 142.3919677734375],
        [139.78387451171875, 140.50253295898438, 141.17147827148438, 141.79678344726562, 142.38360595703125],
        [139.77781677246094, 140.4949951171875, 141.16250610351562, 141.78646850585938, 142.37196350097656]],dtype=float)


latitudes =  np.array([[55.61929702758789, 55.621070861816406, 55.61888122558594, 55.613487243652344, 55.60547637939453],
             [55.53120040893555, 55.532840728759766, 55.53053665161133, 55.525047302246094, 55.5169677734375],
             [55.44305419921875, 55.444580078125, 55.44219207763672, 55.43663024902344, 55.42848587036133],
             [55.35470199584961, 55.356109619140625, 55.353614807128906, 55.34796905517578, 55.33975601196289],
             [55.26683807373047, 55.268131256103516, 55.26553726196289, 55.25981140136719, 55.25152587890625]],dtype=float)

data =  np.array([[10, 10, 10, 10, 10],
        [20, 20, 20, 20, 20],
        [30, 30, 30, 30, 30],
        [40, 40, 40, 40, 40],
        [50, 50, 50, 50, 50]],dtype=float)


x = longitudes.ravel()
y = latitudes.ravel()
z = data.ravel()

xMin, xMax = np.min(x), np.max(x)
yMin, yMax = np.min(y), np.max(y)

xi = np.linspace(xMin, xMax, 0.005)  ##choosen spacing of 0.005
yi = np.linspace(yMin, yMax, 0.005)  ##choosen spacing of 0.005

数据并不完全是一个网格。其实我无法想象如何提前做到这一点：

zi_matplotlib = griddata(x, y, z, xi, yi, interp='linear')

from scipy.interpolate import griddata  ##Using scipy method

zi_scipy = griddata((x, y), z, (xi, yi), method='nearest')    

plt.imshow(????)

请有任何想法和解决方案。

Answer 1

您可以使用插值将扭曲的网格转换为常规网格。插值拟合原始数据点并返回一个可以在您选择的任何点进行评估的函数，在这种情况下，您将选择一个规则的点网格。

这是一个例子：

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

# your data here, as posted in the question

f = interp2d(lon, lat, data, kind="cubic", bounds_error=False)

dlon, dlat = 1.2, .2
xlon = np.linspace(min(lon.flat), max(lon.flat), 20)
xlat = np.linspace(min(lat.flat), max(lat.flat), 20)

# the next few lines are because there seems to be a bug in interp2d
#  instead one would just want to use   r = interp2d(X.flat, Y.flat) (where X,Y are as below)
#  but for the version of scipy I'm using ('0.13.3'), this throws an exception.
r = np.zeros((len(xlon), len(xlat)))
for i, rlat in enumerate(xlat):
    for j, rlon in enumerate(xlon):
        r[i,j] = f(rlon, rlat)

X, Y = np.meshgrid(xlon, xlat)
plt.imshow(r, interpolation="nearest", origin="lower", extent=[min(xlon), max(xlon), min(xlat), max(xlat)], aspect=6.)

plt.scatter(lon.flat, lat.flat, color='k')
plt.show()

在这里，我留下了相当粗糙的网格 (20x20) 并使用了interpolation="nearest"，因此您仍然可以看到代表每个插值的彩色方块，当然，在常规网格上完成（使用两个 @987654324 创建@电话）。还要注意使用或origin="lower"，它将图像和散点图设置为具有相同的方向。

要解释这一点，主要问题是从左到右改变值。这是因为数据在水平点集上被指定为常数，但是因为这些指定的点被扭曲，插值在它们移动时会缓慢变化。例如，右侧最低的散射点应与左侧最高的散射点具有大致相同的颜色。此外，这表明最左边的两对之间没有太大的颜色变化，但在最右边的两对之间变化很大，翘曲最大。

请注意，可以对任何值进行插值，而不仅仅是常规网格，根据原始问题，该网格仅用于 imshow。另请注意，我使用了bounds_error=False，因此我可以在原始数据集之外稍微评估几个点，但要非常小心，因为原始数据之外的点将很快变得不合理，因为在它们所在的区域之外评估立方很合适。

Answer 2

假设longitudes 和latitudes 等间距，您可以直接使用imshow，因为它具有插值功能：

import numpy as np
import matplotlib.pyplot as plt

longitudes = np.array([[139.79391479492188, 140.51760864257812, 141.19119262695312, 141.82083129882812, 142.41165161132812],
        [139.79225158691406, 140.51416015625, 141.18606567382812, 141.8140869140625, 142.40338134765625],
        [139.78591918945312, 140.50637817382812, 141.17694091796875, 141.80377197265625, 142.3919677734375],
        [139.78387451171875, 140.50253295898438, 141.17147827148438, 141.79678344726562, 142.38360595703125],
        [139.77781677246094, 140.4949951171875, 141.16250610351562, 141.78646850585938, 142.37196350097656]],dtype=float)


latitudes =  np.array([[55.61929702758789, 55.621070861816406, 55.61888122558594, 55.613487243652344, 55.60547637939453],
             [55.53120040893555, 55.532840728759766, 55.53053665161133, 55.525047302246094, 55.5169677734375],
             [55.44305419921875, 55.444580078125, 55.44219207763672, 55.43663024902344, 55.42848587036133],
             [55.35470199584961, 55.356109619140625, 55.353614807128906, 55.34796905517578, 55.33975601196289],
             [55.26683807373047, 55.268131256103516, 55.26553726196289, 55.25981140136719, 55.25152587890625]],dtype=float)

data =  np.array([[10, 10, 10, 10, 10],
        [20, 20, 20, 20, 20],
        [30, 30, 30, 30, 30],
        [40, 40, 40, 40, 40],
        [50, 50, 50, 50, 50]],dtype=float)

extent = (longitudes[0,0], longitudes[0,-1], latitudes[0,0], latitudes[-1,0])
plt.imshow(data, interpolation='bilinear', extent=extent, aspect='auto')
plt.show()

我知道这并没有完全回答您的问题。但我认为这是解决根本问题的简单方法。

---

编辑

我刚刚意识到您的数据实际上并不完全是一个网格，而是几乎。您必须决定是否仍要使用我的解决方案...

Answer 3

这是一个使用您的数据的 3d 散点图示例，用各自的彩色标记将每组纬度/经度数据分解成其自己的系列。

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt


longitudes = np.array([[139.79391479492188, 140.51760864257812, 141.19119262695312, 141.82083129882812, 142.41165161132812],
        [139.79225158691406, 140.51416015625, 141.18606567382812, 141.8140869140625, 142.40338134765625],
        [139.78591918945312, 140.50637817382812, 141.17694091796875, 141.80377197265625, 142.3919677734375],
        [139.78387451171875, 140.50253295898438, 141.17147827148438, 141.79678344726562, 142.38360595703125],
        [139.77781677246094, 140.4949951171875, 141.16250610351562, 141.78646850585938, 142.37196350097656]],dtype=float)


latitudes =  np.array([[55.61929702758789, 55.621070861816406, 55.61888122558594, 55.613487243652344, 55.60547637939453],
             [55.53120040893555, 55.532840728759766, 55.53053665161133, 55.525047302246094, 55.5169677734375],
             [55.44305419921875, 55.444580078125, 55.44219207763672, 55.43663024902344, 55.42848587036133],
             [55.35470199584961, 55.356109619140625, 55.353614807128906, 55.34796905517578, 55.33975601196289],
             [55.26683807373047, 55.268131256103516, 55.26553726196289, 55.25981140136719, 55.25152587890625]],dtype=float)

data =  np.array([[10, 10, 10, 10, 10],
        [20, 20, 20, 20, 20],
        [30, 30, 30, 30, 30],
        [40, 40, 40, 40, 40],
        [50, 50, 50, 50, 50]],dtype=float)

colors = ['r','g','b','k','k']
markers = ['o','o','o','o','^']

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

for i in range(5):
    ax.scatter(longitudes[i], latitudes[i], data[i], c=colors[i], marker=markers[i])

ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.set_zlabel('Data')

plt.show()

这会产生类似的图像