定位区域中的点

我有一个数据框，其中两列是点的坐标。 如果该点位于特定位置，我需要用特定值填充一列（全无）。该位置和标签存储在另一个 df 中

解释起来并不容易，但我希望通过一个例子就可以清楚地理解： DF1

   latitude  longitude  LABEL
0    1.3       2.7      None
1    3.5       3.6      None
2    2.8       3.0      None
3    9.7       1.9      None
4    6.2       5.7      None
5    1.7       3.4      None
6    3.5       1.4      None
7    2.7       6.6      None
8    1.7       2.7      None
9    1.3       1.3      None

DF2

   minlat     maxlat    minlong   maxlong  STRING
0    1.0       2.0        1.0       3.0     AAA
1    3.0       4.0        1.0       2.0     BBB
2    3.0       4.0        3.0       4.0     CCC
3    5.0       7.0        2.0       3.0     DDD

最终结果是：

   latitude  longitude  LABEL
0    1.3       2.7      AAA
1    3.5       3.6      CCC
2    2.8       3.0      None
3    9.7       1.9      None
4    6.2       5.7      None
5    1.7       3.4      None
6    3.5       1.4      BBB
7    2.7       6.6      None
8    1.7       2.7      AAA
9    1.3       1.3      AAA

现在的代码是：

for i in range(len(df2)-1):
DF1.loc[(DF1['latitude']>=DF2.loc[i:i,'minlat'].at[i]) & (DF1['latitude']<DF2.loc[i:i,'maxlat'].at[i]) &
   (DF1['longitude']>=DF2.loc[i:i,'minlong'].at[i]) & (DF1['longitude']<DF2.loc[i:i,'maxlong'].at[i]),'LABEL'] = DF2.loc[i:i,'STRING'].at[i]

屏幕有更好的缩进：

因此，对于 DF2 的每一行，我检查 DF1 的值是否在中间，并分配一个字符串

但这样需要很多时间。你对我能做什么有什么建议吗？ 我的问题是 NUMBER_1 的每个值都必须检查 DF2 的每一行，而不仅仅是具有相同索引的行。

编辑：我正在尝试其他方法：

2)

for i in range(len(xlsx_fact_maneuver_specialareas)-1):
    minLat=DF2.loc[i:i,'minLat'].at[i]
    maxLat=DF2.loc[i:i,'maxLat'].at[i]
    minLong=DF2.loc[i:i,'maxLat'].at[i]
    maxLong=DF2.loc[i:i,'maxLong'].at[i]
    DF1.loc[(DF1['latitude']>=minLat) & (DF1['latitude']<maxLat) &
   (DF1['longitude']>=minLong) & (DF1['longitude']<maxLong),'LABEL'] = DF2.loc[i:i,'STRING'].at[i]

这使我在本地花费的时间更少，但当我在机器上尝试时花费的时间更多。

and

for i in range(len(xlsx_fact_maneuver_specialareas)-1):
    minLat=DF2.loc[i:i,'minLat'].at[i]
    maxLat=DF2.loc[i:i,'maxLat'].at[i]
    minLong=DF2.loc[i:i,'maxLat'].at[i]
    maxLong=DF2.loc[i:i,'maxLong'].at[i]
    DF1 = DF1.assign(
        label =  np.select(
          [(DF1['latitude']>=minLat) & (DF1['latitude']<maxLat) & (DF1['longitude']>=minLong) & (DF1['longitude']<maxLong)],
          [DF2.loc[i:i,'STRING'].at[i]],
          [None]))

这让我在本地花费更多，但在机器上花费更少

这是另一个答案，这个答案使用scipy.spatial.KDTree https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.KDTree.html。在边界框之间没有太多重叠并且它们的“半径”分布不太“疯狂”的情况下，它特别有效（radius.max() / np.median(radius)不要太大，例如1 和 2 之间）。

它适用于p-范数 1（曼哈顿）或范数 2（欧几里得），尽管在实践中p=2更快，因为平均每个点看到的候选点更少（圆圈的面积总和小于菱形的面积）。

为什么这么快？KD-trees https://en.wikipedia.org/wiki/K-d_tree非常适合此类问题。他们通过沿维度和中点在每个节点处分割空间来划分空间。一旦构建完成，由于它们提供了分而治之的方法，查询它们的速度很快。

关键函数如下：

def find(kd_boxes, kd_points, bb, r, p):
    # bb: m bounding boxes, in the form [x0,x1,y0,y1]
    # find one bb for each point (the first one), or m if none
    xy = kd_points.data
    found = np.full(len(xy), len(bb), dtype=int)
    cand = pad_jagged(kd_points.query_ball_tree(kd_boxes, r * 1.001, p=p))

    remaining = None
    for c in range(cand.shape[1]):
        remaining = np.nonzero(cand[:, c] >= 0)[0] if remaining is None else remaining[np.nonzero(cand[remaining, 1] >= 0)[0]]
        if remaining.size == 0:
            break
        i = remaining
        j = cand[i, c]
        ok = (bb[j, 0::2] <= xy[i]).all(1) & (xy[i] <= bb[j, 1::2]).all(1)
        k = np.nonzero(ok)[0]
        found[i[k]] = j[k]
        remaining = remaining[~ok]
    return found

在调用此函数之前，我们计算每个边界框的“半径”，该半径是对角线 p 范数的一半。然后我们使用整体最大半径r作为最大距离KDTree查询。这KDTree query (kd_points.query_ball_tree()）有效地筛选所有边界框中心，并在一次调用中找到半径内的所有边界框中心。这是实际匹配的超集，但速度很快并且大大减少了搜索空间。之后，就是过滤那些点的问题了actually位于边界框中，并跟踪每个点的（第一个）匹配边界框。

作为一种优化（此处未实现），我们可以考虑边界框的大小（radius大批）。如果判断太过不同，则边界框可以分为两组（例如，围绕np.median(radius)），并且可以对每一半进行相同的搜索（如果需要，再次递归）。

对于 OP 示例，经过一些准备以更易于使用的形式（所有 Numpy 数组）获取中心、边界框和半径，该查询返回：

# preparation
xy = df1[['longitude', 'latitude']].values
bb = df2[['minlong', 'maxlong', 'minlat', 'maxlat']].values
centers = bb @ np.array([[1,1,0,0], [0,0,1,1]]).T / 2
radius = np.linalg.norm(bb @ np.array([[-1,1,0,0], [0,0,-1,1]]).T, axis=1) / 2
kd_boxes = KDTree(centers)
kd_points = KDTree(xy)

# query
>>> kd_points.query_ball_tree(kd_boxes, radius.max() * 1.001)
[[0], [2], [2], [], [], [], [1], [], [0], [0]]

使用曼哈顿范数可以获得类似的结果（在本例中是相同的）：

radius = bb @ np.array([-1,1,-1,1]) / 2
>>> kd_points.query_ball_tree(kd_boxes, radius.max() * 1.001, p=1)
[[0], [2], [2], [], [], [], [1], [], [0], [0]]

这两种解决方案都绘制在下图中，其中所有候选边界框以及距离内的实际区域都突出显示r他们的中心：

请注意，在这两种情况下，如何指出#2被错误地包含在 bbox 中#2。第一步，使用 KD 树，只是一个过滤步骤。下一步是检查每个点的候选人actually包含要点。通过该过滤（表达式为ok数组），解决办法是：

让我们看看如果我们有更多的点和边界框（可能有重叠）会发生什么：

def gen_example(n, m, seed=None):
    # let's make a more complex problem, with overlaps
    if seed is not None:
        np.random.seed(seed)
    df1 = pd.DataFrame({
        'latitude': np.random.uniform(0, 1 + 1 / np.sqrt(m), size=n),
        'longitude': np.random.uniform(0, 1 + 1 / np.sqrt(m), size=n),
    })

    df2 = pd.DataFrame({
        'minlat': np.random.uniform(size=m),
        'maxlat': np.random.uniform(.5, 1, size=m) / np.sqrt(m),
        'minlong': np.random.uniform(size=m),
        'maxlong': np.random.uniform(.5, 1, size=m) / np.sqrt(m),
        'STRING': [f'b_{i}' for i in range(m)],
    })
    df2['maxlat'] += df2['minlat']
    df2['maxlong'] += df2['minlong']
    return df1, df2

For df1, df2 = gen_example(32, 12, 0)，对应的图片为：

和候选人（查询的结果，此处为p=2) are:

>>> kd_cand = kd_points.query_ball_tree(kd_boxes, r * 1.001, p=2)
[[], [], [9], [], [], [10], [], [], [], [], [], [9], [10], [], [11], [11], [], [2, 6], [8], [8], [], [4], [7, 9], [4], [3, 5], [2, 4], [0], [], [7, 9], [7, 9], [0, 3, 5], [4]]

由于这是一个锯齿状数组，我们将其转换为矩形数组fill_value=-1:

def pad_jagged(a, fill_value=-1):
    lens = np.array([len(r) for r in a])
    v = np.array([e for r in a for e in r])
    w = np.max(lens)
    return np.insert(
        v,
        np.repeat(np.cumsum(lens), w - lens),
        fill_value
    ).reshape(-1, w)

对于上面的填充数组，这给出：

>>> pad_jagged(kd_cand)
[[-1 -1 -1]
 [ 8 -1 -1]
 [ 9 -1 -1]
 ...
 [ 7  9 -1]
 [ 0  3  5]
 [ 4 -1 -1]]

现在，我们迭代该数组的列，但以贪婪的方式删除先前迭代中的任何成功匹配（这就是remaining).

处理预处理等的其他函数有：

def find_regions(xy, bb, p=2):
    # xy: numpy array (n, 2)
    # bb: numpy array (m, 4): the four columns are xmin, xmax, ymin, ymax
    # for each point in xy, return the index of a region that contains it, or -1
    centers = bb @ np.array([[1,1,0,0], [0,0,1,1]]).T / 2
    assert p in {1,2}, "only 1- or 2-norm supported"
    radius = bb @ np.array([-1,1,-1,1]) / 2 if p == 1 else np.linalg.norm(bb @ np.array([[-1,1,0,0], [0,0,-1,1]]).T, axis=1) / 2
    kd_boxes = KDTree(centers)
    kd_points = KDTree(xy)
    return find(kd_boxes, kd_points, bb, radius.max(), p=p)

def find_region_label(df1, df2, nomatch=None, p=2):
    found = find_regions(
        df1[['longitude', 'latitude']].values,
        df2[['minlong', 'maxlong', 'minlat', 'maxlat']].values,
        p=p,
    )
    lbl = np.r_[df2['STRING'].values, [nomatch]]
    return lbl[found]

在OP示例中：

>>> df1.assign(LABEL=find_region_label(df1, df2))
   latitude  longitude LABEL
0       1.3        2.7   AAA
1       3.5        3.6   CCC
2       2.8        3.0  None
3       9.7        1.9  None
4       6.2        5.7  None
5       1.7        3.4  None
6       3.5        1.4   BBB
7       2.7        6.6  None
8       1.7        2.7   AAA
9       1.3        1.3   AAA

速度测试

现在进行速度测试，其中两个点数n and边界框的数量m更大：

df1, df2 = gen_example(100_000, 10_000, 0)

速度比较@norok2 的 numba 解决方案 https://stackoverflow.com/a/68926522/758174（稍作修改以遵循OP的列名称）：

a = %timeit -o find_region_label(df1, df2)
# 222 ms ± 904 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

b = %timeit -o locate_in_regions_nb(df1, df2)
# 5.38 s ± 40.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

>>> b.average / a.average
24.255

该数据的速度提高了 24 倍。

验证我们是否得到相同的解决方案：

sa = find_region_label(df1, df2)
sb = locate_in_regions_nb(df1, df2)
>>> (sa == sb).all()
True