In [1]:

    

import os , glob
import matplotlib.pyplot as plt
import sewpy
import pyfits , math
import random
import numpy as np

from scipy.stats import sigmaclip

## Constants
BEAMARCSEC= 3600. 

## directories
rootdir = "/home/stephane/Science/ETGs/RadioGalaxy/calibrators"
wdir    = "%s/products"%(rootdir)
plotdir = "%s/products/plots"%(rootdir)
fitsrootdir = "/home/stephane/Science/ETGs/RadioGalaxy/calibrators/ALMA/analysis/bosscha/tertiary/J1315-5334/"

os.chdir(wdir)


    

    




/home/stephane/.local/lib/python3.7/site-packages/pyfits/__init__.py:22: PyFITSDeprecationWarning: PyFITS is deprecated, please use astropy.io.fits
  PyFITSDeprecationWarning)  # noqa






    




---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
<ipython-input-1-62dbc04d2e5d> in <module>
     17 fitsrootdir = "/home/stephane/Science/ETGs/RadioGalaxy/calibrators/ALMA/analysis/bosscha/tertiary/J1315-5334/"
     18 
---> 19 os.chdir(wdir)

FileNotFoundError: [Errno 2] No such file or directory: '/home/stephane/Science/ETGs/RadioGalaxy/calibrators/products'

In [ ]:

    

## Get the list of all fits images assuming they have to be extracted.

def getFitsFiles(dir):
    list = glob.glob("%s/**/*fits"%(dir), recursive=True)
    
    return(list)


    

In [ ]:

    

def runsextractor(image_file, detect_thresh=3.5, analysis_thresh=3.0):
    params = ['NUMBER', 'FLUX_ISO', 'FLUXERR_ISO', 'FLUX_AUTO', 'FLUXERR_AUTO', 'FLUX_BEST', 'FLUXERR_BEST', 'BACKGROUND', 
              'THRESHOLD', 'FLUX_MAX', 'XMAX_IMAGE', 'YMAX_IMAGE', 'XPEAK_IMAGE', 'YPEAK_IMAGE', 'ALPHAPEAK_J2000', 
              'DELTAPEAK_J2000', 'X_IMAGE', 'Y_IMAGE', 'ALPHA_SKY', 'DELTA_SKY', 'ALPHA_J2000', 'DELTA_J2000']

    config = {"DETECT_THRESH":detect_thresh, "ANALYSIS_THRESH":analysis_thresh}

    sew = sewpy.SEW(workdir="./", sexpath="/usr/bin/sextractor",params=params, config=config)

    out = sew(image_file)
    data = out["table"]
    
    ra, dec, flux, label = data['ALPHA_J2000'], data['DELTA_J2000'], data['FLUX_MAX'], data['NUMBER'].astype('int')
    
    return ra, dec, flux, label


    

In [ ]:

    

# fits inforation

def getARfits(fitsimage):
    
    hdulist = pyfits.open(fitsimage)
    #hdulist.info()
    
    hdr = hdulist[0].header
    nx = int(hdr['NAXIS1'])
    ny = int(hdr['NAXIS2'])
    bmin = float(hdr['BMIN']) * BEAMARCSEC
    bmaj = float(hdr['BMAJ']) * BEAMARCSEC
    ar = math.sqrt(bmin*bmaj)
    lonpole = float(hdr['CRVAL1'])
    latpole = float(hdr['CRVAL2'])
    
    #print(hdr)
    data = hdulist[0].data
    
    return(nx,ny, bmin, bmaj , ar , lonpole , latpole)


    

In [ ]:

    

## RMS estimation

def getRMSfits(fitsimage , nsample=40 , size= 20):
    
    random.seed()
    
    hdulist = pyfits.open(fitsimage)
    
    hdr = hdulist[0].header
    nx = int(hdr['NAXIS1'])
    ny = int(hdr['NAXIS2'])
    
    data = hdulist[0].data
    
    rms = []
    for i in range(nsample):
        ix = random.randint(size,int(nx/2)-size)
        iy = random.randint(size,int(ny/2)-size)
        fx = random.random()
        fy = random.random()
        
        if fx < 0.5:
            ixcenter = ix
        else:
            ixcenter = nx - ix
            
        if fy < 0.5:
            iycenter = iy
        else:
            iycenter = ny - iy
            
        datrms = data[0, 0, ixcenter-size:ixcenter+size, iycenter-size:iycenter+size]
        rms.append(np.std(datrms))
        
    
    c, low, upp = sigmaclip(rms , 3.0 , 3.0)

    rmsestimated = np.median(c)    
    return(rmsestimated)


    

In [ ]:

    

## plots the sources...

def plotSources(plotdir, fitsimage, data , alpha , delta, fov, rms):
    ra  = data['ra']
    dec = data['dec']
    flux = data['flux']
    
    plt.xlim(alpha-fov/2., alpha+fov/2.)
    plt.ylim(delta-fov/2., delta+fov/2.)
    #plt.plot(ra,dec,"r*")
    # plt.plot(alpha, delta, "b+" , markersize= 10)
    plt.scatter(alpha, delta, cmap='hsv', alpha=0.75)
    
    plt.text(alpha + 0.2*fov,delta+ 0.45*fov,"RMS= %3.3f mJy"%(rms*1e3))
    figname= "%s/%s-sources.png"%(plotdir, fitsimage)
    plt.savefig(figname)
    plt.show()


    

In [ ]:

    

def main(plot=True):
    lfits = getFitsFiles(fitsrootdir)
    
    res = []
    
    print("## Running sextractor...")
    for fitsimage in lfits:
        print("### Fits: %s"%(fitsimage))
        
        nx , ny , bmin, bmaj, ar, alpha, delta = getARfits(fitsimage)
        print("## NX= %d , NY= %d"%(nx,ny))
        print("## BMIN= %3.3f arcsec, BMAJ= %3.3f arcsec, AR= %3.3f arcsec"%(bmin, bmaj,ar))
        
        rms = getRMSfits(fitsimage, 40, 30)
        print("## RMS= %3.4f mJy"%(rms*1e3))
        
        dat= {}
        ra, dec , flux , label = runsextractor(fitsimage)
        dat['fits'] = fitsimage
        dat['rms'] = rms
        dat['ra']   = ra
        dat['dec']  = dec
        dat['flux'] = flux
        
        res.append(dat)  
        
        if plot:
            dr= 100./3600.   ### FOV of 100 arcsec
            plotSources(plotdir, fitsimage.split("/")[-1], dat , alpha, delta, dr , rms)
       
    return(res)


    

In [ ]:

    

allres = main()