Skip to content
Snippets Groups Projects
Select Git revision
  • 53f9fc90571ece5f81fdd5c392fd041534f7ee73
  • master default protected
  • gec-95-update-api-docs
  • gec-95-add-filter-skymodel
  • gec-105-lint-skymodel-filter
  • add-coverage-badge
  • releases/1.6
  • releases/1.5
  • rap-457-implement-ci-v2
  • test_heroku
  • action-python-publish
  • v1.7.1
  • v1.7.0
  • v1.6.2
  • v1.6.1
  • v1.6.post1
  • v1.6
  • v1.5.1
  • v1.5.0
  • v1.4.11
  • v1.4.10
  • v1.4.9
  • v1.4.8
  • v1.4.7
  • v1.4.6
  • v1.4.5
  • v1.4.4
  • v1.4.3
  • v1.4.2
  • v1.4.1
  • v1.4.0
31 results

_cluster.py

Blame
  • David Rafferty's avatar
    David Rafferty authored
    53f9fc90
    History
    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)