Jak mogę dopasować podobne współrzędne przy użyciu Pythona?

Question

Tło:

I nadano cztery katalogi danych, z których (nazwijmy CAT1) pierwszy daje współrzędne (w rektascensji i deklinacji, RA i DEC) dla źródeł radiowych w polach 1 i 2 , drugi katalog (Cat2) podaje RA i Dec dla źródeł radiowych i źródeł podczerwieni (IR) w polu 1, trzeci katalog (Cat3) podaje RA i Dec dla źródeł radiowych i źródeł podczerwieni w polu 2, a czwarty katalog (Cat4) daje RA i Dec dla źródeł optycznych w polach 1 i 2.

Cat1 ma około 2000 źródeł dla pola 2, pamiętając, że niektóre źródła są faktycznie mierzone wielokrotnie w swoich wymiarach, np. ; źródło 1, źródło 2, źródło 3a, źródło 3b, źródło 3c, źródło 4 ... Cat1 ma około 3000 źródeł dla pola 1, również z częściami źródłowymi. Cat 1 to plik .dat, który otwieram w textedit i konwertuję na .txt do użycia z np.genfromtxt.

Cat2 ma około 1700 źródeł pola 1. Cat3 ma około 1700 źródeł pola 2. Cat2 i Cat3 .csv są pliki, które mam otwarcia w liczbach.

Cat4 ma około 1200 źródeł dla pola 1 i około 700 źródeł dla pola 2. Cat4 to plik .dat, który otwieram w textedit i konwertuję na .txt do użycia z np.genfromtxt.

Dowiedz się także, jak konwertować Cat1 i Cat4 w plikach .csv.

Cel:

Celem jest, aby połączyć te cztery katalogi w jednym katalogu, który daje RA i Dec z cat2, CAT1 i Cat4 (pole 1), wówczas RA i Dec z Cat3, Cat1 i Cat4 (pole 2), takie, że RA i Dec z Cat1 i Cat4 są najbliżej RA i Dec z Cat1 lub Cat2, tak że można powiedzieć, że najprawdopodobniej są one tym samym źródłem. Poziom nakładania się będzie różny, ale wygenerowałem wykresy rozrzutu dla danych, które pokazują, że istnieje odpowiednie źródło Cat1 i Cat4 dla każdego źródła Cat2 i Cat3, z dokładnością wielkości znacznika, z oczywiście wieloma resztkami źródła Cat1 i Cat4, które zawierają znacznie więcej źródeł niż Cat2 i Cat3.

Sztuczka polega na tym, że ponieważ współrzędne nie są dokładnie takie same, muszę najpierw spojrzeć na RA i znaleźć najlepszą pasującą wartość, a następnie spojrzeć na odpowiednią Dec i sprawdzić, czy jest to również najlepiej odpowiadająca wartość .

przykład w przypadku źródła Cat2: RA = 53,13360595 Dec = -28.0530758

cat1: RA = 53,133496, gru = -27,553401 lub RA = 53,133873, gru = -28,054600

Tutaj 53.1336 jest równe między 53.1334 a 53.1338, ale wyraźnie -28.053 jest bliższe -28.054 niż -27.553, więc druga opcja w Cat1 jest zwycięzcą.

Postęp:

Do tej pory dopasowane pierwszych 15 wartości w cat2 do wartości w CAT1 czysto oko (Command + F do jak największej liczby miejsc po przecinku, jak to możliwe, a następnie przy użyciu najlepszych wyroku), ale wyraźnie jest to niezwykle nieefektywne dla wszystkich 3400 źródeł w Cat2 i Cat3.Chciałem tylko przekonać się, jakiego rodzaju dokładności oczekiwać w dopasowaniu, i niestety, niektóre pasują tylko do drugiego lub trzeciego miejsca po przecinku, podczas gdy inne pasują do czwartego lub piątego.

W produkcji moje wykresy rozrzutu, użyłem kodu:

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

Wtedy po prostu wykreślono RA_cat1 przeciwko Dec_cat1 i powtarza się dla wszystkich moich katalogach.

Mój problem polega na tym, że szukając odpowiedzi na pytanie, w jaki sposób wytworzyć kod, który będzie w stanie dopasować swoje współrzędne, widziałem wiele odpowiedzi, które wymagają konwersji tablic na listy, ale gdy próbuję to zrobić za pomocą

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

Skończyłem z listą formularza;

[RA1, RA2, RA3, ..., RA3000], [dEC1, DEC2, Dec3, ..., Dec3000]

zamiast co zakładam byłby bardziej pomocny;

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

Co więcej, w przypadku podobnych pytań najbardziej użytecznymi odpowiedziami po pomyślnym przekonwertowaniu listy wydają się być używanie słowników, ale nie jestem pewien, jak używać słownika do tworzenia porównań, które opisałem powyżej.

Dodatkowo, powinienem wspomnieć, że gdy już mi się to udało, poproszono mnie o powtórzenie procesu dla znacznie większego zestawu danych, nie jestem pewien, o ile większy, ale przypuszczam, że może dziesiątki tysięcy zestawów współrzędnych.

Answer 1

Dla ilości posiadanych danych można obliczyć miarę odległości między każdą parą punktów. Coś jak:

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

dla dużego zestawu danych możesz użyć algorytmu linii przesuwu tylko obliczyć miary dla punktów, które są w tej samej okolicy. Tak:

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

wygenerować pewne dane testowe

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 

Oto algorytm:

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

na moim laptopie, stwierdzi ścisły punkt w < dziesiątej sekundy.