如何粗粒xy坐标

I have an xy array of gps coordinates that are not evenly sampled. That seems obvious but I would like to project it onto a grid. Here is my script for the context:

我有一个gy数组的gps坐标，没有均匀采样。这似乎很明显，但我想把它投射到网格上。这是我的上下文脚本：

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)

How do I proceed from here to get the (xnew, ynew) values interpolated on the (xi, yi) grid given the (x, y) original coordinates?

在给定（x，y）原始坐标的情况下，如何从这里开始在（xi，yi）网格上插入（xnew，ynew）值？

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

I thought I would use something similar to the matplotlib function griddata:

我以为我会使用类似于matplotlib函数griddata的东西：

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

Which I can't do; I want to grid the coordinates themselves, not the value z = f(x,y). Do you know of any fast solution that would do it?

哪个我不能做;我想网格坐标本身，而不是值z = f（x，y）。你知道任何快速解决方案吗？

EDIT: It really helps. You guys rock! To be more accurate in my request, I have taken the plot you've generated and marked the (xnew, ynew) samples I am expecting to get with hollow circles. They fall on the (xi, yi) grid.

编辑：这真的很有帮助。你们好棒！为了在我的请求中更准确，我已经拍摄了你生成的情节并标记了我期望用空心圆圈得到的（xnew，ynew）样本。它们落在（xi，yi）网格上。

2 个解决方案

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

yields

产量

---

xi, yi define your grid.

xi，yi定义你的网格。

xx, yy are points on the interpolated path.

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

To find points on the vertical grid lines, first use np.searchsorted to find the index (of xi) where each value in xx (the x-data) would fit into the sorted array xi (the grid).

要在垂直网格线上找到点，首先使用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])

We can then take the difference of neighboring values:

然后我们可以得出相邻值的差异：

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

Where the value is non-zero indicates that the path has crossed a grid line.

值不为零表示路径已越过网格线。

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)

The x_mask can then be used to find the x-values of points on the interpolated path closest to a grid line:

然后可以使用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])

and of course we can find the corresponding y-values similarly:

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

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

Above I've walked through the steps used to find the crossing points of the vertical grid lines; in the code the same was done for the horizontal grid lines.

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

#2

I hope the revision below helps... I have added the required comments.

我希望下面的修订有助于...我添加了所需的评论。

from __future__ import division
import numpy as np
from matplotlib.mlab import griddata
from scipy.interpolate import InterpolatedUnivariateSpline 
from matplotlib import pyplot as plt

# Data
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)

# Get mean values 
x_mean = np.mean(x)
y_mean = np.mean(y)

# Create the fit data for spline
x_fit_data = x/x_mean
y_fit_data = y/y_mean

# A spline of 3. degree
factor = 1000 # arbitrary sampling factor
g = np.arange(0, len(x))
fx = InterpolatedUnivariateSpline(g, x_fit_data, k=3)
fy = InterpolatedUnivariateSpline(g, y_fit_data, k=3)

# Get new data points
gg = np.linspace(0, len(x), len(x)*factor)
x_new = fx(gg) * x_mean
y_new = fy(gg) * y_mean

# Plot
fig = plt.figure()
plt.axis('equal') # optional. To set equal aspect ratio
ax = fig.add_subplot(111)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_xlim(xmin,xmax)
ax.set_ylim(ymin,ymax)
ax.plot(x_new, y_new)
ax.plot(x, y, 'ro') # if you also want to plot the sampling points
ax.grid()
plt.show()

I think you should split your question into two: interpolation and plotting.

我认为你应该将你的问题分成两个：插值和绘图。