python - 如何粗化xy坐标

Question

我有一个未均匀采样的 gps 坐标的 xy 数组。这似乎很明显，但我想将它投影到网格上。这是我的上下文脚本:

import numpy as np
from matplotlib.mlab import griddata

gps_track = np.array([[0,0],[1.2,2.3],[1.9,3],[3.2,4.3],[4,2.9],[6.5,3.1]])
x = gps_track[:,0]
y = gps_track[:,1]

# define grid
binsize = 1
xmin, xmax = x.min(), x.max()
ymin, ymax = y.min(), y.max()
xi = np.arange(xmin, xmax+binsize, binsize)
yi = np.arange(ymin, ymax+binsize, binsize)

在给定 (x, y) 原始坐标的情况下，我如何从这里开始获取在 (xi, yi) 网格上插值的 (xnew, ynew) 值？

# grid the data
xnew, ynew = grid(x, y, xi, yi)

我想我会使用类似于 matplotlib 函数 griddata 的东西:

zi = griddata(x, y, z, xi, yi)

我做不到；我想对坐标本身进行网格化，而不是对值 z = f(x,y) 进行网格化。您知道有什么快速解决方案可以做到这一点吗？

编辑:这真的很有帮助。你们好棒!为了更准确地满足我的要求，我采用了您生成的图并用空心圆圈标记了我希望获得的 (xnew, ynew) 样本。它们落在 (xi, yi) 网格上。

Answer 1

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

gps_track = np.array([[0,0],[1.2,2.3],[1.9,3],[3.2,4.3],[4,2.9],[6.5,3.1]])
x = gps_track[:,0]
y = gps_track[:,1]

binsize = 1.0
xmin, xmax = x.min(), x.max()
ymin, ymax = y.min(), y.max()
xi = np.arange(xmin, xmax+binsize, binsize)
yi = np.arange(ymin, ymax+binsize, binsize)

N = 1000
tckp, u = interpolate.splprep([x, y], s=0, k=2, nest=-1)
xx, yy = interpolate.splev(np.linspace(0, 1, N), tckp)

# find integer indices where the xx fall between xi values
x_idx = np.searchsorted(xi, xx)
# find where the x_idx indices change. This is where the curve crosses a grid line 
x_mask = np.diff(x_idx) != 0
# do the same for y
y_idx = np.searchsorted(yi, yy)
y_mask = np.diff(y_idx) != 0

fig, ax = plt.subplots()
line, = ax.plot(xx, yy)
ax.scatter(xx[x_mask], yy[x_mask], color='red')
ax.scatter(xx[y_mask], yy[y_mask], color='red')

# make sure the grid lines fall on xi and yi
ax.set_xticks(xi)
ax.set_yticks(yi)
ax.grid()
plt.show()

产量

---

xi, yi 定义你的网格。

xx, yy 是插值路径上的点。

要查找垂直网格线上的点，首先使用np.searchsorted 查找索引(xi)，其中xx(x 数据)中的每个值都适合排序数组 xi(网格)。

In [88]: x_idx = np.searchsorted(xi, xx)
In [89]: x_idx
Out[89]: 
array([0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
       1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3,
       3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
       4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6,
       6, 6, 6, 6, 7, 7, 7, 7])

然后我们可以取相邻值的差值:

In [91]: np.diff(x_idx)
Out[91]: 
array([1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 1, 0, 0, 0])

非零值表示路径已经穿过网格线。

In [93]: x_mask = (np.diff(x_idx) != 0)

In [94]: x_mask
Out[94]: 
array([ True, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False,  True, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False,  True, False, False, False, False, False, False,
       False, False, False, False, False, False, False,  True, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
        True, False, False, False, False, False, False, False, False,
       False, False, False, False,  True, False, False, False, False,
       False, False, False, False, False,  True, False, False, False], dtype=bool)

然后可以使用 x_mask 找到插值上的点的 x 值最接近网格线的路径:

In [95]: xx[x_mask]
Out[95]: 
array([ -6.93660834e-18,   9.65180188e-01,   1.95810449e+00,
         2.94969940e+00,   3.94167195e+00,   4.92607812e+00,
         5.99593850e+00])

当然我们可以类似地找到相应的 y 值:

In [96]: yy[x_mask]
Out[96]: 
array([ 0.        ,  2.02032307,  3.06460807,  4.25374305,  2.97366674,
        2.4770819 ,  2.79951527])

上面我已经介绍了用于查找垂直网格线交叉点的步骤；在代码中，水平网格线也是如此。