2016-02-20 6 views
6

背景:私はフィールド1及び2の無線ソースのデータの4社のカタログ、そのうちの最初の(のは、CAT1を呼び出してみましょう()、赤経および赤緯にRAとDEC)の座標を与えるが与えられているPythonを使用して同様の座標をどのように一致させることができますか?

第2カタログ(CAT2)は、フィールド1のラジオソースおよび赤外線(IR)ソースのRAおよびDECを、第3カタログ(Cat3)はフィールド2のラジオおよびIRソースのRAおよびDecを、 Cat1)は、フィールド1と2の光源のRAとDecを示します。

Cat1には、フィールド2の約2000のソースがあります。 ;ソース1、ソース2、ソース3a、ソース3b、ソース3c、ソース4 ... Cat1はフィールド1に約3000のソースを持ちます。 Cat 1は.datファイルで、これはtexteditで開き、np.genfromtxtで使用する.txtに変換されます。

Cat2のフィールド1のソースは約1700です。 Cat3のフィールド2のソースは約1700です。 Cat2とCat3はNumbersで開いています.csvファイルです。

Cat4には、フィールド1のソースが約1200、フィールド2のソースが約700あります。Cat4はtexteditで開いている.datファイルで、np.genfromtxtで使用するために.txtに変換されます。

Cat1とCat4を.csvファイルに変換する方法もわかりました。

目的:

目標はCAT3からCat2の、CAT1とCAT4(フィールド1)、その後、RAおよび12月からRAおよび12月を与える1つのカタログにこれら四つのカタログを組み合わせることで、 Cat1とCat4のRAとDecがRAに最も近く、Cat1またはCat2のDecが同じソースとなる可能性が高いと言えるように、Cat1とCat4(フィールド2) 重複のレベルは変わりますが、プロットマーカーのサイズの精度内で、それぞれのCat2およびCat3ソースに対応するCat1およびCat4ソースがあることを示すデータの散布図を作成しました。 Cat1とCat4のソースには、Cat2とCat3よりも多くのソースが含まれています。

トリックは、座標が正確に一致しないため、まずRAを見て最適な一致する値を見つけて、それに対応するDecを見て、それが最良の対応値であることを確認する必要があります。

例えば、Cat2のソースの場合:RA = 53.13360595、12月= -28.0530758

CAT1:RA = 53.133496、12月= -27.553401 またはRA = 53.133873、ここでは12月= -28.054600

、 53.1336は53.1334と53.1338の間で等しくなりますが、明らかに-28.053は-27.553より-28.054に近く、Cat1の2番目の選択肢が勝者です。

進捗状況:

は、これまでのところ、私は(そして、最善の判断を使用して、できるだけ多くの小数点以下の桁にFコマンド+)目で純粋CAT1の値にCat2の最初の15の値と一致したが、明確にこれは、Cat2とCat3のすべての3400のソースに対して非常に非効率的です。私はちょうどマッチングにどのような正確さが期待されているのか自分自身で見たいと思っていました。残念なことに、小数点第2位または第3位に一致するものもあれば、第4または第5にマッチするものもあります。私の散布図を生成するには

が、私はコードを使用:次に

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

を単にDec_cat1に対してRA_cat1をプロットし、すべての私のカタログのために繰り返しました。

私の問題は、座標を一致させるコードをどのように生成するかについての回答を検索する際に、配列をリストに変換することに関連する多くの答えを見てきましたが、

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

私は最終的にフォームのリストを作成します。

[RA1、RA2、RA3、...、RA3000]、[DEC1、DEC2、DEC3、...、Dec3000]

代わりに、私はより参考になると仮定し何の。

[RA1、DEC1]、[RA2、DEC2]、...、[RA3000、Dec3000]。

また、同様の質問のために、リストの変換が成功した場合の最も有用な答えは、辞書を使用するように見えますが、辞書を使用して上記で説明した種類の比較を生成する方法は不明です。

さらに、私がこの作業に成功すると、かなり大きなデータセットの処理を繰り返すように求められていますが、どれほど大きいか分かりませんが、数千の座標セットの

答えて

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

、それは第二の<第十で近いポイントを検索します。

関連する問題