Come faccio a far corrispondere coordinate simili usando Python?

Question

Background:

mi è stato dato quattro cataloghi di dati, il primo dei quali (chiamiamolo Cat1) dà le coordinate (in ascensione retta e declinazione, AR e Dec) per sorgenti radio in campi 1 e 2 , il secondo catalogo (Cat2) fornisce RA e Dec per le sorgenti radio e le sorgenti a infrarossi (IR) nel campo 1, il terzo catalogo (Cat3) fornisce RA e Dec per le sorgenti radio e IR nel campo 2 e il quarto catalogo (Cat4) fornisce RA e Dec per le sorgenti ottiche nei campi 1 e 2.

Cat1 ha circa 2000 sorgenti per il campo 2, tenendo presente che alcune delle sorgenti sono effettivamente misurate più volte attraverso le loro dimensioni, ad es. ; source 1, source 2, source 3a, source 3b, source 3c, source 4 ... Cat1 ha circa 3000 sorgenti per il campo 1, sempre con alcune fonti in parti. Cat 1 è un file .dat, che sto aprendo in textedit e convertito in .txt per l'uso con np.genfromtxt.

Cat2 ha circa 1700 fonti per il campo 1. Cat3 ha circa 1700 fonti per il campo 2. Cat2 e Cat3 sono file .csv, che sto aprendo in Numbers.

Cat4 ha circa 1200 sorgenti per il campo 1 e circa 700 fonti per il campo 2. Cat4 è un file .dat, che sto aprendo in textedit e convertito in .txt per l'uso con np.genfromtxt.

Ha anche calcolato come convertire Cat1 e Cat4 in file .csv.

Obiettivo:

L'obiettivo è quello di combinare queste quattro cataloghi in un unico catalogo, che dà la RA e Dec da Cat2, Cat1 e Cat4 (campo 1), allora la RA e Dec da Cat3, Cat1 e Cat4 (campo 2), in modo che RA e Dec da Cat1 e Cat4 siano i più vicini a RA e Dec da Cat1 o Cat2, in modo tale che si possa dire che è altamente probabile che siano la stessa fonte. Il livello di sovrapposizione varierà, ma ho prodotto grafici a dispersione per i dati che mostrano che esiste una sorgente Cat1 e Cat4 corrispondente per ciascuna fonte Cat2 e Cat3, entro l'accuratezza della dimensione dell'indicatore di trama, con ovviamente molti avanzi fonti in Cat1 e Cat4, che contengono molte più fonti di Cat2 e Cat3.

Il trucco è che poiché le coordinate non corrispondono esattamente, devo essere in grado di guardare prima RA e trovare il valore di corrispondenza migliore, quindi guardare il Dec corrispondente e controllare che sia anche il valore corrispondente migliore .

esempio, per una fonte Cat2: RA = 53,13,360595 millions, Dic = -28,0530758

Cat1: RA = 53,133,496 mila, Dic = -27,553401 o RA = 53,133,873 mila, Dic = -28,054600

Qui, 53.1336 è ugualmente tra 53.1334 e 53.1338, ma chiaramente -28.053 è più vicino a -28.054 di -27.553, quindi la seconda opzione in Cat1 è la vincente.

Progress:

Finora ho abbinato i primi 15 valori in Cat2 ai valori in Cat1 puro occhio (Comando + F per il maggior numero di decimali possibili, quindi utilizzando miglior giudizio), ma chiaramente questo è estremamente inefficiente per tutte le 3400 sorgenti tra Cat2 e Cat3.Volevo solo vedere di persona che tipo di accuratezza di essere in attesa della corrispondenza, e, purtroppo, un po 'di partite solo di seconda o terza cifra decimale, mentre altri match per quarto o quinto.

Nel produrre i miei grafici a dispersione, ho usato il codice:

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

Poi basta tramato contro RA_cat1 Dec_cat1, e ripetuto per tutti i miei cataloghi.

Il mio problema ora è che nella ricerca di risposte su come produrre un codice che sarà in grado di abbinare le mie coordinate, ho visto molte risposte che implicano la conversione degli array in liste, ma quando si tenta di farlo usando

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

Finisco con un elenco del modulo;

[RA1, RA2, RA3, ..., RA3000], [Dec1, Dec2, DEC3, ..., Dec3000]

invece di quello che presumo sarebbe più utile;

[RA1, Dec1], [RA2, Dec2], ..., [RA3000, Dec3000].

Inoltre, per le domande simili, le risposte più utili una volta che la conversione della lista ha avuto successo sembrano essere utilizzare i dizionari, ma non sono chiaro su come utilizzare un dizionario per produrre i tipi di confronti che ho descritto sopra.

Inoltre, dovrei menzionare che una volta che ho avuto successo in questo compito, mi è stato chiesto di ripetere il processo per un insieme considerevolmente più grande di dati, non sono sicuro quanto più grande, ma sto assumendo forse decine di migliaia di coordinate insiemi.

Answer 1

Per la quantità di dati che si possiede, si può calcolare una distanza metrica tra ogni coppia di punti. Qualcosa di simile:

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

Per un grande insieme di dati si consiglia di utilizzare un algoritmo di linea spazzata solo calcolare la metrica per i punti che si trovano nello stesso quartiere. Come questo:

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

generare alcuni dati di prova

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 

Ecco l'algoritmo:

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

Sul mio computer portatile, si trova il punto vicino a < decimo di secondo.