如何使用 Python 匹配相似的坐标？

背景：

我收到了四个数据目录，其中第一个目录（我们称之为 Cat1）给出了场 1 和 2 中无线电源的坐标（赤经和赤纬、RA 和 Dec），第二个目录（Cat2）给出了 RA和 Dec 适用于领域 1 中的无线电源和红外 (IR) 源，第三个目录 (Cat3) 给出了领域 2 中的无线电和红外源的 RA 和 Dec，第四个目录 (Cat4) 给出了光源的 RA 和 Dec在字段 1 和 2 中。

Cat1 具有大约 2000 个场 2 源，请记住，某些源实际上在其维度上进行了多次测量，例如；源 1、源 2、源 3a、源 3b、源 3c、源 4... Cat1 大约有 3000 个场 1 源，其中一些源还是部分的。 Cat 1 是一个 .dat 文件，我在 textedit 中打开该文件，并将其转换为 .txt 以与 np.genfromtxt 一起使用。

Cat2 具有字段 1 的大约 1700 个源。 Cat3 拥有大约 1700 个字段 2 的源。 Cat2 和 Cat3 是 .csv 文件，我在 Numbers 中打开它们。

Cat4 大约有 1200 个字段 1 的源，大约 700 个字段 2 的源。Cat4 是一个 .dat 文件，我在 textedit 中打开该文件，并将其转换为 .txt 以与 np.genfromtxt 一起使用。

还弄清楚了如何转换 .csv 文件中的 Cat1 和 Cat4。

Aim:

目标是将这四个目录合并为一个目录，该目录给出来自 Cat2、Cat1 和 Cat4（字段 1）的 RA 和 Dec，然后给出来自 Cat3、Cat1 和 Cat4（字段 2）的 RA 和 Dec，这样 RA来自 Cat1 和 Cat4 的 RA 和 Dec 与来自 Cat1 或 Cat2 的 RA 和 Dec 最接近，因此可以说它们很可能是同一源。 重叠程度会有所不同，但我已经为数据生成了散点图，显示每个 Cat2 和 Cat3 源都有相应的 Cat1 和 Cat4 源，在绘图标记大小的精度范围内，当然还有许多剩余源在 Cat1 和 Cat4 中，它们比 Cat2 和 Cat3 包含更多的源。

诀窍在于，因为坐标不完全匹配，所以我需要能够首先查看 RA 并找到最佳匹配值，然后查看相应的 Dec，并检查它是否也是最佳对应值。

例如，对于 Cat2 中的源：RA = 53.13360595，Dec = -28.0530758

类别 1：RA = 53.133496，12 月 = -27.553401 或 RA = 53.133873，十二月 = -28.054600

在这里，53.1336 等于 53.1334 和 53.1338 之间，但显然 -28.053 比 -27.553 更接近 -28.054，因此 Cat1 中的第二个选项是获胜者。

进步：

到目前为止，我已经纯粹通过肉眼将 Cat2 中的前 15 个值与 Cat1 中的值进行了匹配（command+f 尽可能多的小数位，然后使用最佳判断），但显然这对于​​跨 Cat2 和 Cat3 的所有 3400 个源来说效率极低。我只是想亲眼看看匹配的准确性如何，不幸的是，有些匹配仅到小数点后第二或第三位，而另一些则匹配到第四或第五位。

在生成散点图时，我使用了代码：

cat1 = np.genfromtext('filepath/cat1.txt', delimiter = '   ')
RA_cat1 = cat1[:,][:,0]
Dec_cat1 = cat1[:,][:,1]

然后简单地绘制 RA_cat1 与 Dec_cat1 的关系，并对我的所有目录重复此操作。

我现在的问题是，在寻找如何生成能够匹配我的坐标的代码的答案时，我看到了很多涉及将数组转换为列表的答案，但是当尝试使用

cat1list = np.array([RA_cat1, Dec_cat1])
cat1list.tolist()

我最终得到了一份表格清单；

[RA1、RA2、RA3、...、RA3000]、[Dec1、Dec2、Dec3、...、Dec3000]

而不是我认为更有帮助的；

[RA1，Dec1]，[RA2，Dec2]，...，[RA3000，Dec3000]。

此外，对于类似的问题，列表转换成功后最有用的答案似乎是使用字典，但我不清楚如何使用字典来产生我上面描述的那种比较。

此外，我应该提到，一旦我成功完成了这项任务，我就被要求对相当大的数据集重复该过程，我不确定有多大，但我假设可能有数万个坐标套。

对于您拥有的数据量，您可以计算每对点之间的距离度量。就像是：

def close_enough(p1, p2):
    # You may need to scale the RA difference with dec. 
    return (p1.RA - p2.RA)**2 + (p1.Dec - p2.Dec)**2) < 0.01

candidates = [(p1,p2) for p1,p2 in itertools.combinations(points, 2)
              if close_enough(p1,p2)]

对于大型数据集，您可能希望使用线扫描算法来仅计算同一邻域中的点的度量。像这样：

import itertools as it
import operator as op
import sortedcontainers     # handy library on Pypi
import time

from collections import namedtuple
from math import cos, degrees, pi, radians, sqrt
from random import sample, uniform

Observation = namedtuple("Observation", "dec ra other")

生成一些测试数据

number_of_observations = 5000
field1 = [Observation(uniform(-25.0, -35.0),     # dec
                      uniform(45.0, 55.0),       # ra
                      uniform(0, 10))            # other data
          for shop_id in range(number_of_observations)]

# add in near duplicates
number_of_dups = 1000
dups = []
for obs in sample(field1, number_of_dups):
    dDec = uniform(-0.0001, 0.0001)
    dRA  = uniform(-0.0001, 0.0001)
    dups.append(Observation(obs.dec + dDec, obs.ra + dRA, obs.other))

data = field1 + dups

算法如下：

# Note: dec is first in Observation, so data is sorted by .dec then .ra.
data.sort()

# Parameter that determines the size of a sliding declination window
# and therefore how close two observations need to be to be considered
# observations of the same object.
dec_span = 0.0001

# Result. A list of observation pairs close enough to be considered 
# observations of the same object.
candidates = []

# Sliding declination window.  Within the window, observations are
# ordered by .ra.
window = sortedcontainers.SortedListWithKey(key=op.attrgetter('ra'))

# lag_obs is the 'southernmost' observation within the sliding declination window.
observation = iter(data)
lag_obs = next(observation)

# lead_obs is the 'northernmost' observation in the sliding declination window.
for lead_obs in data:

    # Dec of lead_obs represents the leading edge of window.
    window.add(lead_obs)

    # Remove observations further than the trailing edge of window.
    while lead_obs.dec - lag_obs.dec > dec_span:
        window.discard(lag_obs)
        lag_obs = next(observation)

    # Calculate 'east-west' width of window_size at dec of lead_obs
    ra_span = dec_span / cos(radians(lead_obs.dec))
    east_ra = lead_obs.ra + ra_span
    west_ra = lead_obs.ra - ra_span

    # Check all observations in the sliding window within
    # ra_span of lead_obs.
    for other_obs in window.irange_key(west_ra, east_ra):

        if other_obs != lead_obs:
            # lead_obs is at the top center of a box 2 * ra_span wide by 
            # 1 * ra_span tall.  other_obs is is in that box. If desired, 
            # put additional fine-grained 'closeness' tests here. 
            # For example:
            #    average_dec = (other_obs.dec + lead_obs.dec) / 2
            #    delta_dec = other_obs.dec - lead_obs.dec
            #    delta_ra  = other_obs.ra - lead_obs.ra)/cos(radians(average_dec))
            # e.g. if delta_dec**2 + delta_ra**2 < threshold:
            candidates.append((lead_obs, other_obs))

在我的笔记本电脑上，它可以在不到十分之一秒的时间内找到接近点。