Select Git revision
_cluster.py
David Rafferty authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
_cluster.py 6.80 KiB
# -*- coding: utf-8 -*-
#
# Defines cluster functions used by the group operation
#
# Copyright (C) 2013 - Reinout van Weeren
# Copyright (C) 2013 - Francesco de Gasperin
# Modified by David Rafferty as required for integration into LSMTool
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
import os, sys
import numpy as np
import pylab as pl
import itertools
try:
import pyrap.tables as pt
import lofar.stationresponse as lsr
except ImportError:
pass
class Patch():
def __init__(self, name, ra, dec, flux):
self.name = name
self.ra = ra
self.dec = dec
self.flux = flux
class Cluster():
def __init__(self, name, patch):
self.name = name
self.init_patch = patch
self.patches = []
self.update_cluster_coords()
def add_patch(self, patch):
"""
Add a patch to this cluster
"""
self.patches.append(patch)
def total_cluster_flux(self):
"""
Return total cluster flux
"""
cluster_flux = self.init_patch.flux
for patch in self.patches:
cluster_flux += patch.flux
return cluster_flux
def update_cluster_coords(self):
"""
update the self.centroid_ra, self.centroid_dec, self.mid_ra, self.mid_dec
"""
self.centroid_ra = self.init_patch.ra*self.init_patch.flux
self.centroid_dec = self.init_patch.dec*self.init_patch.flux
min_ra = np.inf
max_ra = -np.inf
min_dec = np.inf
max_dec = -np.inf
for patch in self.patches:
self.centroid_ra += patch.ra*patch.flux
self.centroid_dec += patch.dec*patch.flux
if patch.ra < min_ra: min_ra = patch.ra
if patch.ra > max_ra: max_ra = patch.ra
if patch.dec < min_dec: min_dec = patch.dec
if patch.dec > max_dec: max_dec = patch.dec
self.centroid_ra /= self.total_cluster_flux()
self.centroid_dec /= self.total_cluster_flux()
self.mid_ra = min_ra + (max_ra - min_ra)/2.
self.mid_dec = min_dec + (max_dec - min_dec)/2.
def ratohms_string(ra):
rah, ram, ras = ratohms(ra)
return str(rah) + ':' + str(ram) + ':' +str(round(ras,2))
def dectodms_string(dec):
decd, decm, decs = dectodms(dec)
return str(decd) + '.' + str(decm) + '.' +str(round(decs,2))
def compute_patch_center(data, beam_ms = None):
"""
Return the patches names, central (weighted) RA and DEC and total flux
"""
patch_names = np.unique(data['Name'])
patches = []
# get the average time of the obs, OK for 1st order correction
if beam_ms != None:
t = pt.table(beam_ms, ack=False)
tt = t.query('ANTENNA1==0 AND ANTENNA2==1', columns='TIME')
time = tt.getcol("TIME")
time = min(time) + ( max(time) - min(time) ) / 2.
t.close()
sr = lsr.stationresponse(beam_ms, False, True)
for patch_name in patch_names:
comp_ids = np.where(data['Name'] == patch_name)[0]
patch_ra = 0.
patch_dec = 0.
weights_ra = 0.
weights_dec = 0.
patch_flux = 0.
for comp_id in comp_ids:
comp_ra = data['RA'][comp_id]
comp_dec = data['Dec'][comp_id]
comp_flux = np.float(data['I'][comp_id])
# beam correction
if beam_ms != None:
sr.setDirection(comp_ra*np.pi/180.,comp_dec*np.pi/180.)
# use station 0 to compute the beam and get mid channel
beam = sr.evaluateStation(time,0)
r = abs(beam[int(len(beam)/2.)])
beam = ( r[0][0] + r[1][1] ) / 2.
#print "Beam:", beam,
comp_flux *= beam
# calculate the average weighted patch center, and patch flux
patch_flux += comp_flux
patch_ra += comp_ra * comp_flux
patch_dec += comp_dec * comp_flux
weights_ra += comp_flux
weights_dec += comp_flux
patches.append(Patch(patch_name, patch_ra/weights_ra, patch_dec/weights_dec, patch_flux))
return patches
def angsep2(ra1, dec1, ra2, dec2):
"""Returns angular separation between two coordinates (all in degrees)"""
from astropy.coordinates import FK5
import astropy.units as u
coord1 = FK5(ra1, dec1, unit=(u.degree, u.degree))
coord2 = FK5(ra2, dec2, unit=(u.degree, u.degree))
return coord1.separation(coord2)
def create_clusters(patches_orig, Q):
"""
Clusterize all the patches of the skymodel iteratively around the brightest patches
"""
# sort the patches by brightest first
idx = np.argsort([patch.flux for patch in patches_orig])[::-1] # -1 to reverse sort
patches = list(np.array(patches_orig)[idx])
# initialize clusters with the brightest patches
clusters = []
for i, patch in enumerate(patches[0:Q]):
clusters.append(Cluster('CLUSTER_'+str(i), patch))
# Iterate until no changes in which patch belongs to which cluster
count = 1
patches_seq_old = []
while True:
print 'Iteration', count
for cluster in clusters:
cluster.update_cluster_coords()
# reset patches
cluster.patches = []
# select patches after the first Q
for patch in patches[Q:]:
# finding closest cluster for each patch
angulardistance = np.inf
for cluster in clusters:
adis = abs(angsep2(cluster.centroid_ra, cluster.centroid_dec, patch.ra, patch.dec))
if adis < angulardistance:
angulardistance = adis
closest_cluster = cluster
# add this patch to the closest cluster
closest_cluster.add_patch(patch)
patches_seq = []
for cluster in clusters:
patches_seq.extend(cluster.patches)
count += 1
if patches_seq == patches_seq_old: break
patches_seq_old = patches_seq
# Make output patch column
patchNames = [''] * len(patches)
patchNames_orig = [p.name for p in patches_orig]
for c in clusters:
for p in c.patches:
patchNames[patchNames_orig.index(p.name)] = c.name
return np.array(patchNames)