Làm cách nào để tôi đối sánh các tọa độ tương tự bằng cách sử dụng Python?

Question

Bối cảnh:

Tôi đã được cấp Bốn catalog của dữ liệu, là người đầu tiên trong số đó (chúng ta hãy gọi CAT1) cung cấp cho các tọa độ (ở xích kinh và từ chối, RA và Dec) cho các nguồn phát thanh trong các lĩnh vực 1 và 2 , catalô thứ hai (Cat2) cung cấp RA và Dec cho nguồn radio và các nguồn hồng ngoại (IR) trong trường 1, catalô thứ ba (Cat3) cung cấp RA và Dec cho radio và các nguồn IR trong trường 2 và catalog thứ tư (Cat4) cung cấp RA và tháng 12 cho các nguồn quang trong các trường 1 và 2.

Cat1 có khoảng 2000 nguồn cho trường 2, lưu ý rằng một số nguồn thực sự được đo nhiều lần trên các kích thước của chúng, ví dụ: ; nguồn 1, nguồn 2, nguồn 3a, nguồn 3b, nguồn 3c, nguồn 4 ... Cat1 có khoảng 3000 nguồn cho trường 1, một lần nữa với một số nguồn đang ở trong phần. Cat 1 là một tập tin .dat, mà tôi đang mở trong textedit, và chuyển thành .txt để sử dụng với np.genfromtxt.

Cat2 có khoảng 1700 nguồn cho trường 1. Cat3 có khoảng 1700 nguồn cho trường 2. Cat2 và Cat3 là tệp .csv, mà tôi đang mở bằng số.

Cat4 có khoảng 1200 nguồn cho trường 1 và khoảng 700 nguồn cho trường 2. Cat4 là tệp .dat, mà tôi đang mở trong soạn thảo văn bản và được chuyển đổi thành .txt để sử dụng với np.genfromtxt.

Cũng tìm ra cách chuyển đổi Cat1 và Cat4 thành tệp .csv.

Mục tiêu:

Mục đích là để kết hợp bốn catalog vào một cửa hàng duy nhất, cung cấp cho RA và Dec từ Cat2, CAT1 và Cat4 (lĩnh vực 1), sau đó RA và Dec từ CAT3, Cat1 và Cat4 (trường 2), sao cho RA và Dec từ Cat1 và Cat4 gần nhất với RA và Dec từ Cat1 hoặc Cat2, sao cho có thể nói rằng chúng có nhiều khả năng là cùng một nguồn. Mức độ chồng chéo sẽ thay đổi, nhưng tôi đã tạo ra các ô phân tán cho dữ liệu cho thấy có một nguồn Cat1 và Cat4 tương ứng cho mỗi nguồn Cat2 và Cat3, trong độ chính xác của kích thước của điểm đánh dấu. nguồn trong Cat1 và Cat4, chứa nhiều nguồn hơn Cat2 và Cat3.

Bí quyết là vì tọa độ không khớp chính xác, tôi cần có thể xem RA lần đầu và tìm giá trị phù hợp nhất, sau đó xem tháng 12 tương ứng và kiểm tra xem nó có phải là giá trị tương ứng hay nhất không .

ví dụ: Đối với một nguồn tin trong Cat2: RA = 53,13360595, Tháng = -28,0530758

cat1: RA = 53,133496, Tháng = -27,553401 hoặc RA = 53,133873, Tháng = -28,054600

Ở đây, 53.1336 là như nhau giữa 53.1334 và 53.1338, nhưng rõ ràng -28.053 là gần hơn -28.054 hơn -27.553, do đó, tùy chọn thứ hai trong Cat1 là người chiến thắng.

Progress:

Cho đến nay tôi đã xuất hiện trong 15 giá trị đầu tiên trong Cat2 để giá trị trong CAT1 hoàn toàn bằng mắt (command + f để càng nhiều chữ số thập phân càng tốt, sau đó sử dụng phán đoán tốt nhất), nhưng rõ ràng điều này cực kỳ không hiệu quả đối với tất cả 3400 nguồn trên Cat2 và Cat3.Tôi chỉ muốn xem bản thân mình có độ chính xác như mong đợi trong kết hợp, và không may, một số chỉ khớp với vị trí thập phân thứ hai hoặc thứ ba, trong khi một số khác khớp với thứ tư hoặc thứ năm.

Trong sản xuất tán xạ của tôi, tôi đã sử dụng mã:

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

Sau đó chỉ cần vẽ RA_cat1 chống Dec_cat1, và lặp đi lặp lại cho tất cả các catalog của tôi. Vấn đề của tôi bây giờ là tìm kiếm câu trả lời về cách tạo mã có thể khớp với tọa độ của tôi, tôi đã thấy rất nhiều câu trả lời liên quan đến việc chuyển đổi mảng thành danh sách, nhưng khi cố gắng sử dụng

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

Tôi kết thúc bằng danh sách biểu mẫu;

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

thay vì những gì tôi giả định sẽ là hữu ích hơn;

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

Hơn nữa, đối với các câu hỏi tương tự, câu trả lời hữu ích nhất khi chuyển đổi danh sách thành công có vẻ là sử dụng từ điển, nhưng tôi không rõ cách sử dụng từ điển để tạo ra các loại so sánh mà tôi đã mô tả ở trên. Ngoài ra, tôi nên đề cập rằng một khi tôi đã thành công trong nhiệm vụ này, tôi đã được yêu cầu lặp lại quá trình cho một tập dữ liệu lớn hơn đáng kể, tôi không chắc chắn lớn hơn bao nhiêu, nhưng tôi cho rằng có lẽ hàng chục của hàng nghìn tập tọa độ.

Answer 1

Đối với lượng dữ liệu bạn có, bạn có thể tính toán số liệu khoảng cách giữa mỗi cặp điểm. Một cái gì đó như:

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

Để có bộ dữ liệu lớn, bạn có thể muốn sử dụng thuật toán quét đường để chỉ tính số liệu cho các điểm trong cùng một vùng lân cận. Như thế này:

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

Tạo một số dữ liệu thử nghiệm

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 

Dưới đây là các thuật toán:

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

Mở máy tính xách tay của tôi, nó tìm thấy điểm gần trong < thứ mười của một giây.