Как сопоставить аналогичные координаты с помощью Python?

Question

фона:

Я получил четыре каталога данных, первый из которых (назовём Cat1) дает координаты (в правом восхождения и склонения, RA и Dec) для радиоисточников в полях 1 и 2 , второй каталог (Cat2) дает RA и Dec для радиоисточников и инфракрасных (IR) источников в поле 1, третий каталог (Cat3) дает RA и Dec для радио и ИК-источников в поле 2, а четвертый каталог (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 имеет приблизительно 1700 источников для поля 1. Cat3 имеет приблизительно 1700 источников для поля 2. Cat2 и Cat3 являются .csv-файлами, которые я открываю в Numbers.

Cat4 имеет приблизительно 1200 источников для поля 1 и приблизительно 700 источников для поля 2. Cat4 - это .dat-файл, который я открываю в текстовом формате и преобразовывается в .txt для использования с np.genfromtxt.

Также выяснено, как преобразовать Cat1 и Cat4 в CSV-файлы.

Цель:

Цель состоит в том, чтобы объединить эти четыре каталога в один каталог, который дает RA и Dec от CAT2, Cat1 и Cat4 (поле 1), то RA и Dec от cat3, Cat1 и Cat4 (поле 2), так что RA и Dec от Cat1 и Cat4 находятся ближе всего к RA и Dec от Cat1 или Cat2, так что можно сказать, что они, скорее всего, будут тем же источником. Уровень перекрытия будет меняться, но я создал диаграммы рассеяния для данных, которые показывают, что есть соответствующий источник Cat1 и Cat4 для каждого источника Cat2 и Cat3 с точностью до размера маркера графика, причем, конечно, много оставшихся источников в Cat1 и Cat4, которые содержат гораздо больше источников, чем Cat2 и Cat3.

Фокус в том, что, поскольку координаты не совпадают, мне нужно сначала взглянуть на RA и найти наилучшее совпадающее значение, а затем посмотреть на соответствующий Dec и проверить, что это также наилучшее соответствующее значение ,

например, для источника в CAT2: RA = 53,13360595, декабрь = -28.0530758

Cat1: RA = 53,133496, Dec = -27,553401 или RA = 53,133873, Dec = -28,054600

Здесь 53,1336 равно равно 53,1334 и 53,1338, но, очевидно, -28,053 приближается к -28,054, чем -27,553, поэтому второй вариант в Cat1 является победителем.

Прогресс:

До сих пор я соответствовал первые 15 значений в CAT2 значениям в Cat1 чисто на глаз (Ctrl + F, как много знаков после запятой, как это возможно, то, используя лучшие решения), но очевидно, это крайне неэффективно для всех 3400 источников в Cat2 и Cat3.Я просто хотел убедиться в себе, какую точность можно ожидать при сопоставлении, и, к сожалению, некоторые соответствуют только второму или третьему знаку после запятой, в то время как другие соответствуют четвертому или пятому.

В производстве мои диаграммы рассеяния, я использовал код:

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

Кроме того, для подобных вопросов наиболее полезные ответы после успешного преобразования списка, по-видимому, состоят в использовании словарей, но я не понимаю, как использовать словарь для создания подобных сравнений, описанных выше.

Кроме того, я должен упомянуть, что, как только я был успешным в этой задаче, меня попросили повторить процесс для значительно большего набора данных, я не уверен, насколько больше, но я предполагаю, что, возможно, десятки тысяч координатных наборов.

Answer 1

Для количества данных, которые у вас есть, вы можете рассчитать метрику расстояния между каждой парой точек. Что-то вроде:

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

На моем ноутбуке, он находит близкую точку в < десятой секунды.