diff --git a/.gitignore b/.gitignore
index d3a18ae07e5d8898c3e0987b797d85ff2f769df9..969b2d256e20aca8f48ab1e5a40fa1f2eda83f07 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,5 @@
test_data/*.MS
+Docker/scripts
# Byte-compiled / optimized / DLL files
__pycache__/
diff --git a/Docker/build_docker.sh b/Docker/build_docker.sh
new file mode 100755
index 0000000000000000000000000000000000000000..45d397f27ae40e9026acf46daffcd89eeebf60a4
--- /dev/null
+++ b/Docker/build_docker.sh
@@ -0,0 +1,10 @@
+#!/bin/bash
+BRANCH=production
+REPO_URL=https://github.com/lofar-astron/prefactor
+
+# FETCHES THE SCRIPTS ONLY
+svn checkout ${REPO_URL}/branches/${BRANCH}/scripts
+SCRIPT_PATH=$(realpath ${BASH_SOURCE[0]})
+PATH=$(dirname ${SCRIPT_PATH})
+
+docker build $PATH -t prefactor:latest
diff --git a/Docker/scripts/Ateamclipper.py b/Docker/scripts/Ateamclipper.py
deleted file mode 100755
index da0f12ea2efbc4a1b86166f34243800c17a75121..0000000000000000000000000000000000000000
--- a/Docker/scripts/Ateamclipper.py
+++ /dev/null
@@ -1,64 +0,0 @@
-#!/usr/bin/env python
-
-## changelog
-# W.Williams 2014/11/03 add - to give input/output statistics per channel
-# W.Williams 2014/11/03 fix - statistics per correlation
-# A.Drabent 2019/07/24 write fraction of flagged data into output file (for prefactor3)
-
-import numpy
-import pyrap.tables as pt
-import os, sys
-
-msname = str(sys.argv[1])
-
-cliplevelhba = 5.0
-cliplevellba = 50.0
-
-t = pt.table(msname, readonly=False)
-data = t.getcol('MODEL_DATA')
-flag = t.getcol('FLAG')
-freq_tab= pt.table(msname + '::SPECTRAL_WINDOW')
-freq = freq_tab.getcol('REF_FREQUENCY')
-
-if freq[0] > 100e6:
- cliplevel = cliplevelhba
-if freq[0] < 100e6:
- cliplevel = cliplevellba
-
-
-print('------------------------------')
-print('SB Frequency [MHz]', freq[0]/1e6)
-for chan in range(0,numpy.size(data[0,:,0])):
- print('chan %i : %.5f%% input XX flagged' %( chan, 100.*numpy.sum(flag[:,chan,0] == True)/numpy.size(flag[:,chan,0]) ))
- print('chan %i : %.5f%% input YY flagged' %( chan, 100.*numpy.sum(flag[:,chan,3] == True)/numpy.size(flag[:,chan,3]) ))
-input_flags_xx = 100. * numpy.sum(flag[:,:,0] == True)/numpy.size(flag[:,:,0])
-input_flags_yy = 100. * numpy.sum(flag[:,:,3] == True)/numpy.size(flag[:,:,3])
-print('Total : %.5f%% input XX flagged' %( input_flags_xx ))
-print('Total : %.5f%% input YY flagged' %( input_flags_yy ))
-print('')
-print('Cliplevel used [Jy]', cliplevel)
-print('\n\n')
-
-for pol in range(0,numpy.size(data[0,0,:])):
- for chan in range(0,numpy.size(data[0,:,0])):
- print('Doing polarization,chan', pol, chan)
- idx = numpy.where(abs(data[:,chan,pol]) > cliplevel)
- flag[idx,chan,0] = True
- flag[idx,chan,1] = True
- flag[idx,chan,2] = True
- flag[idx,chan,3] = True
-
-print('')
-for chan in range(0,numpy.size(data[0,:,0])):
- print('chan %i : %.5f%% output XX flagged' %( chan, 100.*numpy.sum(flag[:,chan,0] == True)/numpy.size(flag[:,chan,0]) ))
- print('chan %i : %.5f%% output YY flagged' %( chan, 100.*numpy.sum(flag[:,chan,3] == True)/numpy.size(flag[:,chan,3]) ))
-output_flags_xx = 100. * numpy.sum(flag[:,:,0] == True)/numpy.size(flag[:,:,0])
-output_flags_yy = 100. * numpy.sum(flag[:,:,3] == True)/numpy.size(flag[:,:,3])
-print('Total : %.5f%% input XX flagged' %( output_flags_xx ))
-print('Total : %.5f%% input YY flagged' %( output_flags_yy ))
-print('')
-t.putcol('FLAG', flag)
-t.close()
-freq_tab.close()
-
-os.system('echo ' + str(freq[0]) + ' ' + str(output_flags_xx - input_flags_xx) + ' ' + str(output_flags_yy - input_flags_yy) + ' >> Ateamclipper.txt')
\ No newline at end of file
diff --git a/Docker/scripts/BLsmooth.py b/Docker/scripts/BLsmooth.py
deleted file mode 100755
index 0bec38e690232e9241c9df0d727b68283a0a8e19..0000000000000000000000000000000000000000
--- a/Docker/scripts/BLsmooth.py
+++ /dev/null
@@ -1,168 +0,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-#
-# Copyright (C) 2019 - Francesco de Gasperin
-#
-# 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
-
-# Usage: BLavg.py vis.MS
-# Load a MS, average visibilities according to the baseline lenght,
-# i.e. shorter BLs are averaged more, and write a new MS
-
-import os, sys, time
-import optparse, itertools
-import logging
-import numpy as np
-from scipy.ndimage.filters import gaussian_filter1d as gfilter
-import pyrap.tables as pt
-logging.basicConfig(level=logging.DEBUG)
-
-def addcol(ms, incol, outcol):
- if outcol not in ms.colnames():
- logging.info('Adding column: '+outcol)
- coldmi = ms.getdminfo(incol)
- coldmi['NAME'] = outcol
- ms.addcols(pt.makecoldesc(outcol, ms.getcoldesc(incol)), coldmi)
- if outcol != incol:
- # copy columns val
- logging.info('Set '+outcol+'='+incol)
- pt.taql("update $ms set "+outcol+"="+incol)
-
-opt = optparse.OptionParser(usage="%prog [options] MS", version="%prog 0.1")
-opt.add_option('-f', '--ionfactor', help='Gives an indication on how strong is the ionosphere [default: 0.2]', type='float', default=0.2)
-opt.add_option('-s', '--bscalefactor', help='Gives an indication on how the smoothing varies with BL-lenght [default: 0.5]', type='float', default=0.5)
-opt.add_option('-i', '--incol', help='Column name to smooth [default: DATA]', type='string', default='DATA')
-opt.add_option('-o', '--outcol', help='Output column [default: SMOOTHED_DATA]', type="string", default='SMOOTHED_DATA')
-opt.add_option('-w', '--weight', help='Save the newly computed WEIGHT_SPECTRUM, this action permanently modify the MS! [default: False]', action="store_true", default=False)
-opt.add_option('-r', '--restore', help='If WEIGHT_SPECTRUM_ORIG exists then restore it before smoothing [default: False]', action="store_true", default=False)
-opt.add_option('-b', '--nobackup', help='Do not backup the old WEIGHT_SPECTRUM in WEIGHT_SPECTRUM_ORIG [default: do backup if -w]', action="store_true", default=False)
-opt.add_option('-a', '--onlyamp', help='Smooth only amplitudes [default: smooth real/imag]', action="store_true", default=False)
-opt.add_option('-S', '--smooth', help='Performs smoothing (otherwise column will be only copied)', type="string", default=True)
-(options, msfile) = opt.parse_args()
-
-if msfile == []:
- opt.print_help()
- sys.exit(0)
-msfile = msfile[0]
-
-if not os.path.exists(msfile):
- logging.error("Cannot find MS file.")
- sys.exit(1)
-
-# open input/output MS
-ms = pt.table(msfile, readonly=False, ack=False)
-
-freqtab = pt.table(msfile + '::SPECTRAL_WINDOW', ack=False)
-freq = freqtab.getcol('REF_FREQUENCY')[0]
-freqtab.close()
-wav = 299792458. / freq
-timepersample = ms.getcell('INTERVAL',0)
-
-# check if ms is time-ordered
-times = ms.getcol('TIME_CENTROID')
-if not all(times[i] <= times[i+1] for i in range(len(times)-1)):
- logging.critical('This code cannot handle MS that are not time-sorted.')
- sys.exit(1)
-
-# create column to smooth
-addcol(ms, options.incol, options.outcol)
-# if smoothing should not be performed
-if options.smooth == 'False':
- sys.exit(0)
-
-# retore WEIGHT_SPECTRUM
-if 'WEIGHT_SPECTRUM_ORIG' in ms.colnames() and options.restore:
- addcol(ms, 'WEIGHT_SPECTRUM_ORIG', 'WEIGHT_SPECTRUM')
-# backup WEIGHT_SPECTRUM
-elif options.weight and not options.nobackup:
- addcol(ms, 'WEIGHT_SPECTRUM', 'WEIGHT_SPECTRUM_ORIG')
-
-# iteration on antenna1
-for ms_ant1 in ms.iter(["ANTENNA1"]):
- ant1 = ms_ant1.getcol('ANTENNA1')[0]
- a_ant2 = ms_ant1.getcol('ANTENNA2')
- logging.debug('Working on antenna: %s' % ant1)
-
- a_uvw = ms_ant1.getcol('UVW')
- a_data = ms_ant1.getcol(options.outcol)
- a_weights = ms_ant1.getcol('WEIGHT_SPECTRUM')
- a_flags = ms_ant1.getcol('FLAG')
-
- for ant2 in set(a_ant2):
- if ant1 == ant2: continue # skip autocorr
- idx = np.where(a_ant2 == ant2)
-
- uvw = a_uvw[idx]
- data = a_data[idx]
- weights = a_weights[idx]
- flags = a_flags[idx]
-
- # compute the FWHM
- uvw_dist = np.sqrt(uvw[:, 0]**2 + uvw[:, 1]**2 + uvw[:, 2]**2)
- dist = np.mean(uvw_dist) / 1.e3
- if np.isnan(dist): continue # fix for missing anstennas
-
- stddev = options.ionfactor * (25.e3 / dist)**options.bscalefactor * (freq / 60.e6) # in sec
- stddev = stddev/timepersample # in samples
- logging.debug("%s - %s: dist = %.1f km: sigma=%.2f samples." % (ant1, ant2, dist, stddev))
-
- if stddev == 0: continue # fix for flagged antennas
- if stddev < 0.5: continue # avoid very small smoothing
-
- flags[ np.isnan(data) ] = True # flag NaNs
- weights[flags] = 0 # set weight of flagged data to 0
- del flags
-
- # Multiply every element of the data by the weights, convolve both the scaled data and the weights, and then
- # divide the convolved data by the convolved weights (translating flagged data into weight=0). That's basically the equivalent of a
- # running weighted average with a Gaussian window function.
-
- # set bad data to 0 so nans do not propagate
- data = np.nan_to_num(data*weights)
-
- # smear weighted data and weights
- if options.onlyamp:
- dataAMP = gfilter(np.abs(data), stddev, axis=0)
- dataPH = np.angle(data)
- else:
- dataR = gfilter(np.real(data), stddev, axis=0)#, truncate=4.)
- dataI = gfilter(np.imag(data), stddev, axis=0)#, truncate=4.)
-
- weights = gfilter(weights, stddev, axis=0)#, truncate=4.)
-
- # re-create data
- if options.onlyamp:
- data = dataAMP * ( np.cos(dataPH) + 1j*np.sin(dataPH) )
- else:
- data = (dataR + 1j * dataI)
- data[(weights != 0)] /= weights[(weights != 0)] # avoid divbyzero
-
- #print np.count_nonzero(data[~flags]), np.count_nonzero(data[flags]), 100*np.count_nonzero(data[flags])/np.count_nonzero(data)
- #print "NANs in flagged data: ", np.count_nonzero(np.isnan(data[flags]))
- #print "NANs in unflagged data: ", np.count_nonzero(np.isnan(data[~flags]))
- #print "NANs in weights: ", np.count_nonzero(np.isnan(weights))
-
- a_data[idx] = data
- if options.weight: a_weights[idx] = weights
-
- #logging.info('Writing %s column.' % options.outcol)
- ms_ant1.putcol(options.outcol, a_data)
-
- if options.weight:
- #logging.warning('Writing WEIGHT_SPECTRUM column.')
- ms_ant1.putcol('WEIGHT_SPECTRUM', a_weights)
-
-ms.close()
-logging.info("Done.")
diff --git a/Docker/scripts/InitSubtract_deep_sort_and_compute.py b/Docker/scripts/InitSubtract_deep_sort_and_compute.py
deleted file mode 100755
index 0776b06ac3bc13b9676798e50e02301c140c29ad..0000000000000000000000000000000000000000
--- a/Docker/scripts/InitSubtract_deep_sort_and_compute.py
+++ /dev/null
@@ -1,482 +0,0 @@
-#! /usr/bin/env python
-"""
-Script to sort a list of MSs into frequency-bands, and compute additional values needed for initsubtract
-"""
-import casacore.tables as pt
-import sys, os
-import numpy as np
-from lofarpipe.support.data_map import DataMap
-from lofarpipe.support.data_map import DataProduct
-import argparse
-from argparse import RawTextHelpFormatter
-
-class Band(object):
- """
- The Band object contains parameters needed for each band
- """
- def __init__(self, MSfiles):
- self.files = MSfiles
- self.msnames = [ MS.split('/')[-1] for MS in self.files ]
- self.numMS = len(self.files)
- # Get the frequency info and set name
- sw = pt.table(self.files[0]+'::SPECTRAL_WINDOW', ack=False)
- self.freq = sw.col('REF_FREQUENCY')[0]
- self.nchan = sw.col('NUM_CHAN')[0]
- self.chan_freqs_hz = sw.col('CHAN_FREQ')[0]
- self.chan_width_hz = sw.col('CHAN_WIDTH')[0][0]
- sw.close()
- self.name = str(int(self.freq/1e6))
- # Get the station diameter
- ant = pt.table(self.files[0]+'::ANTENNA', ack=False)
- self.diam = float(ant.col('DISH_DIAMETER')[0])
- ant.close()
-
- def get_image_sizes(self, cellsize_highres_deg=None, cellsize_lowres_deg=None,
- fieldsize_highres=2.5, fieldsize_lowres=6.5):
- """
- Sets sizes for initsubtract images
-
- The image sizes are scaled from the mean primary-beam FWHM. For
- the high-res image, we use 2.5 * FWHM; for low-res, we use 6.5 * FHWM.
-
- Parameters
- ----------
- cellsize_highres_deg : float, optional
- cellsize for the high-res images in deg
- cellsize_lowres_deg : float, optional
- cellsize for the low-res images in deg
- fieldsize_highres : float, optional
- How many FWHM's shall the high-res images be.
- fieldsize_lowres : float, optional
- How many FWHM's shall the low-res images be.
- """
- if cellsize_highres_deg:
- self.cellsize_highres_deg = cellsize_highres_deg
- if cellsize_lowres_deg:
- self.cellsize_lowres_deg = cellsize_lowres_deg
- if not hasattr(self, 'mean_el_rad'):
- for MS_id in range(self.numMS):
- # calculate mean elevation
- tab = pt.table(self.files[MS_id], ack=False)
- el_values = pt.taql("SELECT mscal.azel1()[1] AS el from "
- + self.files[MS_id] + " limit ::10000").getcol("el")
- if MS_id == 0:
- global_el_values = el_values
- else:
- global_el_values = np.hstack( (global_el_values, el_values ) )
- self.mean_el_rad = np.mean(global_el_values)
-
- # Calculate mean FOV
- sec_el = 1.0 / np.sin(self.mean_el_rad)
- self.fwhm_deg = 1.1 * ((3.0e8 / self.freq) / self.diam) * 180. / np.pi * sec_el
- self.imsize_high_res = self.get_optimum_size(self.fwhm_deg
- /self.cellsize_highres_deg * fieldsize_highres)
- self.imsize_low_res = self.get_optimum_size(self.fwhm_deg
- /self.cellsize_lowres_deg * fieldsize_lowres)
- return (self.imsize_high_res, self.imsize_low_res)
-
- def get_optimum_size(self, size):
- """
- Gets the nearest optimum image size
-
- Taken from the casa source code (cleanhelper.py)
-
- Parameters
- ----------
- size : int
- Target image size in pixels
-
- Returns
- -------
- optimum_size : int
- Optimum image size nearest to target size
-
- """
- import numpy
-
- def prime_factors(n, douniq=True):
- """ Return the prime factors of the given number. """
- factors = []
- lastresult = n
- sqlast=int(numpy.sqrt(n))+1
- if n == 1:
- return [1]
- c=2
- while 1:
- if (lastresult == 1) or (c > sqlast):
- break
- sqlast=int(numpy.sqrt(lastresult))+1
- while 1:
- if(c > sqlast):
- c=lastresult
- break
- if lastresult % c == 0:
- break
- c += 1
- factors.append(c)
- lastresult = lastresult // c
- if (factors==[]): factors=[n]
- return numpy.unique(factors).tolist() if douniq else factors
-
- n = int(size)
- if (n%2 != 0):
- n+=1
- fac=prime_factors(n, False)
- for k in range(len(fac)):
- if (fac[k] > 7):
- val=fac[k]
- while (numpy.max(prime_factors(val)) > 7):
- val +=1
- fac[k]=val
- newlarge=numpy.product(fac)
- for k in range(n, newlarge, 2):
- if ((numpy.max(prime_factors(k)) < 8)):
- return k
- return newlarge
-
- def get_averaging_steps(self):
- """
- Sets the averaging step sizes
-
- Note: the frequency step must be an even divisor of the number of
- channels
-
- """
- # Get time per sample and number of samples
- t = pt.table(self.files[0], readonly=True, ack=False)
- for t2 in t.iter(["ANTENNA1","ANTENNA2"]):
- if (t2.getcell('ANTENNA1',0)) < (t2.getcell('ANTENNA2',0)):
- self.timestep_sec = t2[1]['TIME']-t2[0]['TIME'] # sec
- break
- t.close()
- # generate a (numpy-)array with the divisors of nchan
- tmp_divisors = []
- for step in range(self.nchan,0,-1):
- if (self.nchan % step) == 0:
- tmp_divisors.append(step)
- freq_divisors = np.array(tmp_divisors)
-
- # Average to 0.2 MHz per channel and 16 sec per time slot
- initsubtract_freqstep = max(1, min(int(round(0.2 * 1e6 / self.chan_width_hz)), self.nchan))
- initsubtract_freqstep = freq_divisors[np.argmin(np.abs(freq_divisors-initsubtract_freqstep))]
- initsubtract_timestep = max(1, int(round(16.0 / self.timestep_sec)))
-
- return (initsubtract_freqstep, initsubtract_timestep)
-
- def get_nwavelengths(self, cellsize_deg, timestep_sec,):
- """
- Returns nwavelengths for WSClean BL-based averaging
- The value depends on the integration time given the specified maximum
- allowed smearing. We scale it from the imaging cell size assuming normal
- sampling as:
- max baseline in nwavelengths = 1 / theta_rad ~= 1 / (cellsize_deg * 3 * pi / 180)
- nwavelengths = max baseline in nwavelengths * 2 * pi *
- integration time in seconds / (24 * 60 * 60) / 4
- Parameters
- ----------
- cellsize_deg : float
- Pixel size of image in degrees
- timestep_sec : float
- Length of one timestep in seconds
- """
- max_baseline = 1 / (3 * cellsize_deg * np.pi / 180)
- wsclean_nwavelengths_time = int(max_baseline * 2*np.pi * timestep_sec /
- (24 * 60 * 60) / 4)
- return wsclean_nwavelengths_time
-
-
-def main(ms_input, outmapname=None, mapfile_dir=None, cellsize_highres_deg=0.00208, cellsize_lowres_deg=0.00694,
- fieldsize_highres=2.5, fieldsize_lowres=6.5, image_padding=1., y_axis_stretch=1.):
- """
- Check a list of MS files for missing frequencies
-
- Parameters
- ----------
- ms_input : list or str
- List of MS filenames, or string with list, or path to a mapfile
- outmapname: str
- Name of output mapfile
- mapfile_dir : str
- Directory for output mapfile
- cellsize_highres_deg : float, optional
- cellsize for the high-res images in deg
- cellsize_lowres_deg : float, optional
- cellsize for the low-res images in deg
- fieldsize_highres : float, optional
- How many FWHM's shall the high-res images be.
- fieldsize_lowres : float, optional
- How many FWHM's shall the low-res images be.
- image_padding : float, optional
- How much padding shall we add to the padded image sizes.
- y_axis_stretch : float, optional
- How much shall the y-axis be stretched or compressed.
-
- Returns
- -------
- result : dict
- Dict with the name of the generated mapfiles
-
- """
-
- if not outmapname or not mapfile_dir:
- raise ValueError('sort_times_into_freqGroups: outmapname and mapfile_dir are needed!')
- if type(ms_input) is str:
- if ms_input.startswith('[') and ms_input.endswith(']'):
- ms_list = [f.strip(' \'\"') for f in ms_input.strip('[]').split(',')]
- else:
- map_in = DataMap.load(ms_input)
- map_in.iterator = DataMap.SkipIterator
- ms_list = []
- for fname in map_in:
- if fname.startswith('[') and fname.endswith(']'):
- for f in fname.strip('[]').split(','):
- ms_list.append(f.strip(' \'\"'))
- else:
- ms_list.append(fname.strip(' \'\"'))
- elif type(ms_input) is list:
- ms_list = [str(f).strip(' \'\"') for f in ms_input]
- else:
- raise TypeError('sort_into_freqBands: type of "ms_input" unknown!')
-
- cellsize_highres_deg = float(cellsize_highres_deg)
- cellsize_lowres_deg = float(cellsize_lowres_deg)
- fieldsize_highres = float(fieldsize_highres)
- fieldsize_lowres = float(fieldsize_lowres)
- image_padding = float(image_padding)
- y_axis_stretch = float(y_axis_stretch)
-
- msdict = {}
- for ms in ms_list:
- # group all MSs by frequency
- sw = pt.table(ms+'::SPECTRAL_WINDOW', ack=False)
- msfreq = int(sw.col('REF_FREQUENCY')[0])
- sw.close()
- if msfreq in msdict:
- msdict[msfreq].append(ms)
- else:
- msdict[msfreq] = [ms]
- bands = []
- bandfreqs = []
- print("InitSubtract_deep_sort_and_compute.py: Putting files into bands.")
- for MSkey in msdict.keys():
- bands.append( Band(msdict[MSkey]) )
- bandfreqs.append( Band(msdict[MSkey]).freq )
-
- ## min freq gives largest image size for deep image
- bandfreqs = np.array(bandfreqs)
- minfreq = np.min(bandfreqs)
- bandmin = np.argmin(bandfreqs)
- ## need to map the output from wsclean channels to the right frequencies
- ## just put the bands in the right freq order
- wsclean_channum = np.argsort(bandfreqs)
- bands = np.array(bands)
- bands = bands[wsclean_channum]
-
- #minfreq = 1e9
- #for ib, band in enumerate(bands):
- #if band.freq < minfreq:
- #minfreq = band.freq
- #bandmin = ib
-
- group_map = MultiDataMap()
- file_single_map = DataMap([])
- high_size_map = DataMap([])
- low_size_map = DataMap([])
- high_paddedsize_map = DataMap([])
- low_paddedsize_map = DataMap([])
- numfiles = 0
- nbands = np.int(len(bands)/10)
- if nbands > 8:
- nchansout_clean1 = np.int(nbands/4)
- elif nbands > 4:
- nchansout_clean1 = np.int(nbands/2)
- else:
- nchansout_clean1 = np.int(nbands)
- if nchansout_clean1 < 2:
- nchansout_clean1 = 2
-
- (freqstep, timestep) = bands[0].get_averaging_steps()
- int_time_sec = bands[0].timestep_sec * timestep # timestep_sec gets added to band object in get_averaging_steps()
- nwavelengths_high = bands[0].get_nwavelengths(cellsize_highres_deg, int_time_sec)
- nwavelengths_low = bands[0].get_nwavelengths(cellsize_lowres_deg, int_time_sec)
-
- print("InitSubtract_deep_sort_and_compute.py: analyzing data...")
- for band in bands:
- group_map.append(MultiDataProduct('localhost', band.files, False))
- numfiles += len(band.files)
- for filename in band.files:
- file_single_map.append(DataProduct('localhost', filename, False))
- (imsize_high_res, imsize_low_res) = band.get_image_sizes(cellsize_highres_deg, cellsize_lowres_deg,
- fieldsize_highres, fieldsize_lowres)
- imsize_high_res_stretch = band.get_optimum_size(int(imsize_high_res*y_axis_stretch))
- high_size_map.append(DataProduct('localhost', str(imsize_high_res)+" "+str(imsize_high_res_stretch), False))
- imsize_low_res_stretch = band.get_optimum_size(int(imsize_low_res*y_axis_stretch))
- low_size_map.append(DataProduct('localhost', str(imsize_low_res)+" "+str(imsize_low_res_stretch), False))
- imsize_high_pad = band.get_optimum_size(int(imsize_high_res*image_padding))
- imsize_high_pad_stretch = band.get_optimum_size(int(imsize_high_res*image_padding*y_axis_stretch))
- high_paddedsize_map.append(DataProduct('localhost', str(imsize_high_pad)+" "+str(imsize_high_pad_stretch), False))
- imsize_low_pad = band.get_optimum_size(int(imsize_low_res*image_padding))
- imsize_low_pad_stretch = band.get_optimum_size(int(imsize_low_res*image_padding*y_axis_stretch))
- low_paddedsize_map.append(DataProduct('localhost', str(imsize_low_pad)+" "+str(imsize_low_pad_stretch), False))
-
- if band.freq == minfreq:
- deep_imsize_high_res = imsize_high_res
- deep_imsize_high_res_stretch = imsize_high_res_stretch
- deep_imsize_high_pad = imsize_high_pad
- deep_imsize_high_pad_stretch = imsize_high_pad_stretch
- deep_imsize_low_res = imsize_low_res
- deep_imsize_low_res_stretch = imsize_low_res_stretch
- deep_imsize_low_pad = imsize_low_pad
- deep_imsize_low_pad_stretch = imsize_low_pad_stretch
-
- deep_high_size_map = DataMap([DataProduct('localhost', str(deep_imsize_high_res)+" "+str(deep_imsize_high_res_stretch), False)])
- deep_high_paddedsize_map = DataMap([DataProduct('localhost', str(deep_imsize_high_pad)+" "+str(deep_imsize_high_pad_stretch), False)])
- deep_low_size_map = DataMap([DataProduct('localhost', str(deep_imsize_low_res)+" "+str(deep_imsize_low_res_stretch), False)])
- deep_low_paddedsize_map = DataMap([DataProduct('localhost', str(deep_imsize_low_pad)+" "+str(deep_imsize_low_pad_stretch), False)])
- nbands_map = DataMap([DataProduct('localhost', str(nbands), False)])
- nchansout_clean1_map = DataMap([DataProduct('localhost', str(nchansout_clean1), False)])
-
- # get mapfiles for freqstep and timestep with the length of single_map
- freqstep_map = DataMap([])
- timestep_map = DataMap([])
- nwavelengths_high_map = DataMap([])
- nwavelengths_low_map = DataMap([])
-
- for index in range(numfiles):
- # set time and frequency averaging values (in sec and Hz)
- freqstep_map.append(DataProduct('localhost', str(freqstep*bands[0].chan_width_hz), False))
- timestep_map.append(DataProduct('localhost', str(timestep*bands[0].timestep_sec), False))
- nwavelengths_high_map.append(DataProduct('localhost', str(nwavelengths_high), False))
- nwavelengths_low_map.append(DataProduct('localhost', str(nwavelengths_low), False))
-
- groupmapname = os.path.join(mapfile_dir, outmapname)
- group_map.save(groupmapname)
- file_single_mapname = os.path.join(mapfile_dir, outmapname+'_single')
- file_single_map.save(file_single_mapname)
-
- high_sizename = os.path.join(mapfile_dir, outmapname+'_high_size')
- high_size_map.save(high_sizename)
- low_sizename = os.path.join(mapfile_dir, outmapname+'_low_size')
- low_size_map.save(low_sizename)
- high_padsize_name = os.path.join(mapfile_dir, outmapname+'_high_padded_size')
- high_paddedsize_map.save(high_padsize_name)
- low_padsize_name = os.path.join(mapfile_dir, outmapname+'_low_padded_size')
- low_paddedsize_map.save(low_padsize_name)
-
- deep_high_sizename = os.path.join(mapfile_dir, outmapname+'_deep_high_size')
- deep_high_size_map.save(deep_high_sizename)
- deep_low_sizename = os.path.join(mapfile_dir, outmapname+'_deep_low_size')
- deep_low_size_map.save(deep_low_sizename)
- deep_high_padsize_name = os.path.join(mapfile_dir, outmapname+'_deep_high_padded_size')
- deep_high_paddedsize_map.save(deep_high_padsize_name)
- deep_low_padsize_name = os.path.join(mapfile_dir, outmapname+'_deep_low_padded_size')
- deep_low_paddedsize_map.save(deep_low_padsize_name)
-
- nbands_mapname = os.path.join(mapfile_dir, outmapname+'_nbands')
- nbands_map.save(nbands_mapname)
- nchansout_clean1_mapname = os.path.join(mapfile_dir, outmapname+'_nchansout_clean1')
- nchansout_clean1_map.save(nchansout_clean1_mapname)
-
- freqstepname = os.path.join(mapfile_dir, outmapname+'_freqstep')
- freqstep_map.save(freqstepname)
- timestepname = os.path.join(mapfile_dir, outmapname+'_timestep')
- timestep_map.save(timestepname)
- nwavelengths_high_name = os.path.join(mapfile_dir, outmapname+'_nwavelengths_high')
- nwavelengths_high_map.save(nwavelengths_high_name)
- nwavelengths_low_name = os.path.join(mapfile_dir, outmapname+'_nwavelengths_low')
- nwavelengths_low_map.save(nwavelengths_low_name)
-
- result = {'groupmap': groupmapname, 'single_mapfile' : file_single_mapname,
- 'high_size_mapfile' : high_sizename, 'low_size_mapfile' : low_sizename,
- 'high_padsize_mapfile' : high_padsize_name, 'low_padsize_mapfile' : low_padsize_name,
- 'deep_high_size_mapfile' : deep_high_sizename, 'deep_low_size_mapfile' : deep_low_sizename,
- 'deep_high_padsize_mapfile' : deep_high_padsize_name, 'deep_low_padsize_mapfile' : deep_low_padsize_name,
- 'nbands' : nbands_mapname, 'nchansout_clean1' : nchansout_clean1_mapname,
- 'freqstep' : freqstepname, 'timestep' : timestepname, 'nwavelengths_high_mapfile': nwavelengths_high_name,
- 'nwavelengths_low_mapfile': nwavelengths_low_name}
- return result
-
-
-class MultiDataProduct(DataProduct):
- """
- Class representing multiple files in a DataProduct.
- """
- def __init__(self, host=None, file=None, skip=True):
- super(MultiDataProduct, self).__init__(host, file, skip)
- if not file:
- self.file = list()
- else:
- self._set_file(file)
-
- def __repr__(self):
- """Represent an instance as a Python dict"""
- return (
- "{{'host': '{0}', 'file': {1}, 'skip': {2}}}".format(self.host, self.file, str(self.skip))
- )
-
- def __str__(self):
- """Print an instance as 'host:[filelist]'"""
- return ':'.join((self.host, str(self.file)))
-
- def _set_file(self, data):
- try:
- # Try parsing as a list
- if isinstance(data, list):
- self.file = data
- if isinstance(data, DataProduct):
- self._from_dataproduct(data)
- if isinstance(data, DataMap):
- self._from_datamap(data)
-
- except TypeError:
- raise DataProduct("No known method to set a filelist from %s" % str(file))
-
- def _from_dataproduct(self, prod):
- self.host = prod.host
- self.file = prod.file
- self.skip = prod.skip
-
- def _from_datamap(self, inmap):
- filelist = {}
- for item in inmap:
- if not item.host in filelist:
- filelist[item.host] = []
- filelist[item.host].append(item.file)
- self.file = filelist['i am']
-
- def append(self, item):
- self.file.append(item)
-
-
-class MultiDataMap(DataMap):
- """
- Class representing a specialization of data-map, a collection of data
- products located on the same node, skippable as a set and individually
- """
- @DataMap.data.setter
- def data(self, data):
- if isinstance(data, DataMap):
- mdpdict = {}
- data.iterator = DataMap.SkipIterator
- for item in data:
- if not item.host in mdpdict:
- mdpdict[item.host] = []
- mdpdict[item.host].append(item.file)
- mdplist = []
- for k, v in mdpdict.iteritems():
- mdplist.append(MultiDataProduct(k, v, False))
- self._set_data(mdplist, dtype=MultiDataProduct)
- elif isinstance(data, MultiDataProduct):
- self._set_data(data, dtype=MultiDataProduct)
- elif not data:
- pass
- else:
- self._set_data(data, dtype=MultiDataProduct)
-
- def split_list(self, number):
- mdplist = []
- for item in self.data:
- for i in range(0, len(item.file), number):
- chunk = item.file[i:i+number]
- mdplist.append(MultiDataProduct(item.host, chunk, item.skip))
- self._set_data(mdplist)
diff --git a/Docker/scripts/InitSubtract_sort_and_compute.py b/Docker/scripts/InitSubtract_sort_and_compute.py
deleted file mode 100755
index 1bd1443662bcb6abc197c2bd7df5661bdb93d751..0000000000000000000000000000000000000000
--- a/Docker/scripts/InitSubtract_sort_and_compute.py
+++ /dev/null
@@ -1,479 +0,0 @@
-#! /usr/bin/env python
-"""
-Script to sort a list of MSs into frequency-bands, and compute additional values needed for initsubtract
-"""
-import pyrap.tables as pt
-import os
-import numpy as np
-from lofarpipe.support.data_map import DataMap
-from lofarpipe.support.data_map import DataProduct
-
-class Band(object):
- """
- The Band object contains parameters needed for each band
- """
- def __init__(self, MSfiles):
- self.files = MSfiles
- self.msnames = [ MS.split('/')[-1] for MS in self.files ]
- self.numMS = len(self.files)
- # Get the frequency info and set name
- sw = pt.table(self.files[0]+'::SPECTRAL_WINDOW', ack=False)
- self.freq = sw.col('REF_FREQUENCY')[0]
- self.nchan = sw.col('NUM_CHAN')[0]
- self.chan_freqs_hz = sw.col('CHAN_FREQ')[0]
- self.chan_width_hz = sw.col('CHAN_WIDTH')[0][0]
- sw.close()
- self.name = str(int(self.freq/1e6))
- # Get the station diameter
- ant = pt.table(self.files[0]+'::ANTENNA', ack=False)
- self.diam = float(ant.col('DISH_DIAMETER')[0])
- ant.close()
-
-
-
- def get_image_sizes(self, cellsize_highres_deg=None, cellsize_lowres_deg=None,
- fieldsize_highres=2.5, fieldsize_lowres=6.5):
- """
- Sets sizes for initsubtract images
-
- The image sizes are scaled from the mean primary-beam FWHM. For
- the high-res image, we use 2.5 * FWHM; for low-res, we use 6.5 * FHWM.
-
- Parameters
- ----------
- cellsize_highres_deg : float, optional
- cellsize for the high-res images in deg
- cellsize_lowres_deg : float, optional
- cellsize for the low-res images in deg
- fieldsize_highres : float, optional
- How many FWHM's shall the high-res images be.
- fieldsize_lowres : float, optional
- How many FWHM's shall the low-res images be.
- """
- if cellsize_highres_deg:
- self.cellsize_highres_deg = cellsize_highres_deg
- if cellsize_lowres_deg:
- self.cellsize_lowres_deg = cellsize_lowres_deg
- if not hasattr(self, 'mean_el_rad'):
- for MS_id in xrange(self.numMS):
- # calculate mean elevation
- tab = pt.table(self.files[MS_id], ack=False)
- el_values = pt.taql("SELECT mscal.azel1()[1] AS el from "
- + self.files[MS_id] + " limit ::10000").getcol("el")
- if MS_id == 0:
- global_el_values = el_values
- else:
- global_el_values = np.hstack( (global_el_values, el_values ) )
- self.mean_el_rad = np.mean(global_el_values)
-
- # Calculate mean FOV
- sec_el = 1.0 / np.sin(self.mean_el_rad)
- self.fwhm_deg = 1.1 * ((3.0e8 / self.freq) / self.diam) * 180. / np.pi * sec_el
- self.imsize_high_res = self.get_optimum_size(self.fwhm_deg
- /self.cellsize_highres_deg * fieldsize_highres)
- self.imsize_low_res = self.get_optimum_size(self.fwhm_deg
- /self.cellsize_lowres_deg * fieldsize_lowres)
- return (self.imsize_high_res, self.imsize_low_res)
-
- def get_optimum_size(self, size):
- """
- Gets the nearest optimum image size
-
- Taken from the casa source code (cleanhelper.py)
-
- Parameters
- ----------
- size : int
- Target image size in pixels
-
- Returns
- -------
- optimum_size : int
- Optimum image size nearest to target size
-
- """
- import numpy
-
- def prime_factors(n, douniq=True):
- """ Return the prime factors of the given number. """
- factors = []
- lastresult = n
- sqlast=int(numpy.sqrt(n))+1
- if n == 1:
- return [1]
- c=2
- while 1:
- if (lastresult == 1) or (c > sqlast):
- break
- sqlast=int(numpy.sqrt(lastresult))+1
- while 1:
- if(c > sqlast):
- c=lastresult
- break
- if lastresult % c == 0:
- break
- c += 1
- factors.append(c)
- lastresult /= c
- if (factors==[]): factors=[n]
- return numpy.unique(factors).tolist() if douniq else factors
-
- n = int(size)
- if (n%2 != 0):
- n+=1
- fac=prime_factors(n, False)
- for k in range(len(fac)):
- if (fac[k] > 7):
- val=fac[k]
- while (numpy.max(prime_factors(val)) > 7):
- val +=1
- fac[k]=val
- newlarge=numpy.product(fac)
- for k in range(n, newlarge, 2):
- if ((numpy.max(prime_factors(k)) < 8)):
- return k
- return newlarge
-
- def get_averaging_steps(self):
- """
- Sets the averaging step sizes
-
- Note: the frequency step must be an even divisor of the number of
- channels
-
- """
- # Get time per sample and number of samples
- t = pt.table(self.files[0], readonly=True, ack=False)
- for t2 in t.iter(["ANTENNA1","ANTENNA2"]):
- if (t2.getcell('ANTENNA1',0)) < (t2.getcell('ANTENNA2',0)):
- self.timestep_sec = t2[1]['TIME']-t2[0]['TIME'] # sec
- break
- t.close()
- # generate a (numpy-)array with the divisors of nchan
- tmp_divisors = []
- for step in range(self.nchan,0,-1):
- if (self.nchan % step) == 0:
- tmp_divisors.append(step)
- freq_divisors = np.array(tmp_divisors)
-
- # For initsubtract, average to 0.5 MHz per channel and 20 sec per time
- # slot. Since each band is imaged separately and the smearing and image
- # sizes both scale linearly with frequency, a single frequency and time
- # step is valid for all bands
- initsubtract_freqstep = max(1, min(int(round(0.5 * 1e6 / self.chan_width_hz)), self.nchan))
- initsubtract_freqstep = freq_divisors[np.argmin(np.abs(freq_divisors-initsubtract_freqstep))]
- initsubtract_timestep = max(1, int(round(20.0 / self.timestep_sec)))
-
- return (initsubtract_freqstep, initsubtract_timestep)
-
- def nwavelengths(self, cellsize_highres_deg, cellsize_lowres_deg, initsubtract_timestep):
- int_time_sec = self.timestep_sec * initsubtract_timestep
- max_baseline_in_nwavelenghts_h = 1.0/(cellsize_highres_deg*3.0*np.pi/180.0)
- max_baseline_in_nwavelenghts_l = 1.0/(cellsize_lowres_deg*3.0*np.pi/180.0)
- self.nwavelengths_high = max_baseline_in_nwavelenghts_h*2.0*np.pi*int_time_sec/(24.0*60.0*60.0)/2
- self.nwavelengths_low = max_baseline_in_nwavelenghts_l*2.0*np.pi*int_time_sec/(24.0*60.0*60.0)/2
- return (self.nwavelengths_high, self.nwavelengths_low)
-
-
-def input2bool(invar):
- if invar is None:
- return None
- if isinstance(invar, bool):
- return invar
- elif isinstance(invar, str):
- if invar.upper() == 'TRUE' or invar == '1':
- return True
- elif invar.upper() == 'FALSE' or invar == '0':
- return False
- else:
- raise ValueError('input2bool: Cannot convert string "'+invar+'" to boolean!')
- elif isinstance(invar, int) or isinstance(invar, float):
- return bool(invar)
- else:
- raise TypeError('input2bool: Unsupported data type:'+str(type(invar)))
-
-
-def main(ms_input, outmapname=None, mapfile_dir=None, cellsize_highres_deg=0.00208, cellsize_lowres_deg=0.00694,
- fieldsize_highres=2.5, fieldsize_lowres=6.5, image_padding=1., y_axis_stretch=1.,
- calc_y_axis_stretch=False, apply_y_axis_stretch_highres=True, apply_y_axis_stretch_lowres=True):
- """
- Check a list of MS files for missing frequencies
-
- Parameters
- ----------
- ms_input : list or str
- List of MS filenames, or string with list, or path to a mapfile
- outmapname: str
- Name of output mapfile
- mapfile_dir : str
- Directory for output mapfile
- cellsize_highres_deg : float, optional
- cellsize for the high-res images in deg
- cellsize_lowres_deg : float, optional
- cellsize for the low-res images in deg
- fieldsize_highres : float, optional
- How many FWHM's shall the high-res images be.
- fieldsize_lowres : float, optional
- How many FWHM's shall the low-res images be.
- image_padding : float, optional
- How much padding shall we add to the padded image sizes.
- y_axis_stretch : float, optional
- How much shall the y-axis be stretched or compressed.
- calc_y_axis_stretch : bool, optional
- Adjust the image sizes returned by this script for the mean elevation.
- If True, the value of y_axis_stretch above is ignored
- apply_y_axis_stretch_highres : bool, optional
- Apply the y-axis stretch to the high-res image sizes
- apply_y_axis_stretch_lowres : bool, optional
- Apply the y-axis stretch to the low-res image sizes
-
- Returns
- -------
- result : dict
- Dict with the name of the generated mapfiles
-
- """
- if not outmapname or not mapfile_dir:
- raise ValueError('sort_times_into_freqGroups: outmapname and mapfile_dir are needed!')
- if type(ms_input) is str:
- if ms_input.startswith('[') and ms_input.endswith(']'):
- ms_list = [f.strip(' \'\"') for f in ms_input.strip('[]').split(',')]
- else:
- map_in = DataMap.load(ms_input)
- map_in.iterator = DataMap.SkipIterator
- ms_list = []
- for fname in map_in:
- if fname.startswith('[') and fname.endswith(']'):
- for f in fname.strip('[]').split(','):
- ms_list.append(f.strip(' \'\"'))
- else:
- ms_list.append(fname.strip(' \'\"'))
- elif type(ms_input) is list:
- ms_list = [str(f).strip(' \'\"') for f in ms_input]
- else:
- raise TypeError('sort_into_freqBands: type of "ms_input" unknown!')
-
- cellsize_highres_deg = float(cellsize_highres_deg)
- cellsize_lowres_deg = float(cellsize_lowres_deg)
- fieldsize_highres = float(fieldsize_highres)
- fieldsize_lowres = float(fieldsize_lowres)
- image_padding = float(image_padding)
- y_axis_stretch = float(y_axis_stretch)
- calc_y_axis_stretch = input2bool(calc_y_axis_stretch)
- apply_y_axis_stretch_highres = input2bool(apply_y_axis_stretch_highres)
- apply_y_axis_stretch_lowres = input2bool(apply_y_axis_stretch_lowres)
-
- msdict = {}
- for ms in ms_list:
- # group all MSs by frequency
- sw = pt.table(ms+'::SPECTRAL_WINDOW', ack=False)
- msfreq = int(sw.col('REF_FREQUENCY')[0])
- sw.close()
- if msfreq in msdict:
- msdict[msfreq].append(ms)
- else:
- msdict[msfreq] = [ms]
- bands = []
- print "InitSubtract_sort_and_compute.py: Putting files into bands."
- for MSkey in msdict.keys():
- bands.append( Band(msdict[MSkey]) )
-
- group_map = MultiDataMap()
- file_single_map = DataMap([])
- high_size_map = DataMap([])
- low_size_map = DataMap([])
- high_paddedsize_map = DataMap([])
- low_paddedsize_map = DataMap([])
- numfiles = 0
- for i, band in enumerate(bands):
- print "InitSubtract_sort_and_compute.py: Working on Band:",band.name
- group_map.append(MultiDataProduct('localhost', band.files, False))
- numfiles += len(band.files)
- for filename in band.files:
- file_single_map.append(DataProduct('localhost', filename, False))
- (imsize_high_res, imsize_low_res) = band.get_image_sizes(cellsize_highres_deg, cellsize_lowres_deg,
- fieldsize_highres, fieldsize_lowres)
-
- # Calculate y_axis_stretch if desired
- if calc_y_axis_stretch:
- if i == 0:
- y_axis_stretch = 1.0 / np.sin(band.mean_el_rad)
- print "InitSubtract_sort_and_compute.py: Using y-axis stretch of:",y_axis_stretch
- y_axis_stretch_lowres = y_axis_stretch
- y_axis_stretch_highres = y_axis_stretch
- else:
- y_axis_stretch_lowres = 1.0
- y_axis_stretch_highres = 1.0
-
- # Adjust sizes so that we get the correct ones below
- if not apply_y_axis_stretch_highres:
- y_axis_stretch_highres = 1.0
- if not apply_y_axis_stretch_lowres:
- y_axis_stretch_lowres = 1.0
- imsize_low_res /= y_axis_stretch_lowres
- imsize_high_res /= y_axis_stretch_highres
-
- imsize_high_res = band.get_optimum_size(int(imsize_high_res))
- imsize_high_res_stretch = band.get_optimum_size(int(imsize_high_res*y_axis_stretch_highres))
- high_size_map.append(DataProduct('localhost', str(imsize_high_res)+" "+str(imsize_high_res_stretch), False))
-
- imsize_low_res = band.get_optimum_size(int(imsize_low_res))
- imsize_low_res_stretch = band.get_optimum_size(int(imsize_low_res*y_axis_stretch_lowres))
- low_size_map.append(DataProduct('localhost', str(imsize_low_res)+" "+str(imsize_low_res_stretch), False))
-
- imsize_high_pad = band.get_optimum_size(int(imsize_high_res*image_padding))
- imsize_high_pad_stretch = band.get_optimum_size(int(imsize_high_res*image_padding*y_axis_stretch_highres))
- high_paddedsize_map.append(DataProduct('localhost', str(imsize_high_pad)+" "+str(imsize_high_pad_stretch), False))
-
- imsize_low_pad = band.get_optimum_size(int(imsize_low_res*image_padding))
- imsize_low_pad_stretch = band.get_optimum_size(int(imsize_low_res*image_padding*y_axis_stretch_lowres))
- low_paddedsize_map.append(DataProduct('localhost', str(imsize_low_pad)+" "+str(imsize_low_pad_stretch), False))
-
- print "InitSubtract_sort_and_compute.py: Computing averaging steps."
- (freqstep, timestep) = bands[0].get_averaging_steps()
- (nwavelengths_high, nwavelengths_low) = bands[0].nwavelengths(cellsize_highres_deg,
- cellsize_lowres_deg, timestep)
-
- # get mapfiles for freqstep and timestep with the length of single_map
- freqstep_map = DataMap([])
- timestep_map = DataMap([])
- nwavelengths_high_map = DataMap([])
- nwavelengths_low_map = DataMap([])
- for index in xrange(numfiles):
- freqstep_map.append(DataProduct('localhost', str(freqstep), False))
- timestep_map.append(DataProduct('localhost', str(timestep), False))
- freqstep_map.append(DataProduct('localhost', str(freqstep), False))
- timestep_map.append(DataProduct('localhost', str(timestep), False))
- nwavelengths_high_map.append(DataProduct('localhost', str(nwavelengths_high), False))
- nwavelengths_low_map.append(DataProduct('localhost', str(nwavelengths_low), False))
-
- groupmapname = os.path.join(mapfile_dir, outmapname)
- group_map.save(groupmapname)
- file_single_mapname = os.path.join(mapfile_dir, outmapname+'_single')
- file_single_map.save(file_single_mapname)
- high_sizename = os.path.join(mapfile_dir, outmapname+'_high_size')
- high_size_map.save(high_sizename)
- low_sizename = os.path.join(mapfile_dir, outmapname+'_low_size')
- low_size_map.save(low_sizename)
- high_padsize_name = os.path.join(mapfile_dir, outmapname+'_high_padded_size')
- high_paddedsize_map.save(high_padsize_name)
- low_padsize_name = os.path.join(mapfile_dir, outmapname+'_low_padded_size')
- low_paddedsize_map.save(low_padsize_name)
- freqstepname = os.path.join(mapfile_dir, outmapname+'_freqstep')
- freqstep_map.save(freqstepname)
- timestepname = os.path.join(mapfile_dir, outmapname+'_timestep')
- timestep_map.save(timestepname)
- nwavelengths_high_name = os.path.join(mapfile_dir, outmapname+'_nwavelengths_high')
- nwavelengths_high_map.save(nwavelengths_high_name)
- nwavelengths_low_name = os.path.join(mapfile_dir, outmapname+'_nwavelengths_low')
- nwavelengths_low_map.save(nwavelengths_low_name)
-
- result = {'groupmap': groupmapname, 'single_mapfile' : file_single_mapname,
- 'high_size_mapfile' : high_sizename, 'low_size_mapfile' : low_sizename,
- 'high_padsize_mapfile' : high_padsize_name, 'low_padsize_mapfile' : low_padsize_name,
- 'freqstep' : freqstepname, 'timestep' : timestepname, 'nwavelengths_high_mapfile': nwavelengths_high_name,
- 'nwavelengths_low_mapfile': nwavelengths_low_name}
- return result
-
-
-class MultiDataProduct(DataProduct):
- """
- Class representing multiple files in a DataProduct.
- """
- def __init__(self, host=None, file=None, skip=True):
- super(MultiDataProduct, self).__init__(host, file, skip)
- if not file:
- self.file = list()
- else:
- self._set_file(file)
-
- def __repr__(self):
- """Represent an instance as a Python dict"""
- return (
- "{{'host': '{0}', 'file': '{1}', 'skip': {2}}}".format(self.host,
- '[{}]'.format(','.join(self.file)), str(self.skip))
- )
-
- def __str__(self):
- """Print an instance as 'host:[filelist]'"""
- return ':'.join((self.host, str(self.file)))
-
- def _set_file(self, data):
- try:
- # Try parsing as a list
- if isinstance(data, list):
- self.file = data
- if isinstance(data, DataProduct):
- self._from_dataproduct(data)
- if isinstance(data, DataMap):
- self._from_datamap(data)
-
- except TypeError:
- raise DataProduct("No known method to set a filelist from %s" % str(file))
-
- def _from_dataproduct(self, prod):
- print 'setting filelist from DataProduct'
- self.host = prod.host
- self.file = prod.file
- self.skip = prod.skip
-
- def _from_datamap(self, inmap):
- print 'setting filelist from DataMap'
- filelist = {}
- for item in inmap:
- if not item.host in filelist:
- filelist[item.host] = []
- filelist[item.host].append(item.file)
- self.file = filelist['i am']
-
- def append(self, item):
- self.file.append(item)
-
-
-class MultiDataMap(DataMap):
- """
- Class representing a specialization of data-map, a collection of data
- products located on the same node, skippable as a set and individually
- """
- def __init__(self, data=list(), iterator=iter):
- super(MultiDataMap, self).__init__(data, iterator)
-
- @classmethod
- def load(cls, filename):
- """Load a data map from file `filename`. Return a DataMap instance."""
- with open(filename) as f:
- datamap = eval(f.read())
- for i, d in enumerate(datamap):
- file_entry = d['file']
- if file_entry.startswith('[') and file_entry.endswith(']'):
- file_list = [e.strip(' \'\"') for e in file_entry.strip('[]').split(',')]
- datamap[i] = {'host': d['host'], 'file': file_list, 'skip': d['skip']}
- return cls(datamap)
-
- @DataMap.data.setter
- def data(self, data):
- if isinstance(data, DataMap):
- mdpdict = {}
- data.iterator = DataMap.SkipIterator
- for item in data:
- if not item.host in mdpdict:
- mdpdict[item.host] = []
- mdpdict[item.host].append(item.file)
- mdplist = []
- for k, v in mdpdict.iteritems():
- mdplist.append(MultiDataProduct(k, v, False))
- self._set_data(mdplist, dtype=MultiDataProduct)
- elif isinstance(data, MultiDataProduct):
- self._set_data(data, dtype=MultiDataProduct)
- elif not data:
- pass
- else:
- self._set_data(data, dtype=MultiDataProduct)
-
- def split_list(self, number):
- mdplist = []
- for item in self.data:
- for i in xrange(0, len(item.file), number):
- chunk = item.file[i:i+number]
- mdplist.append(MultiDataProduct(item.host, chunk, item.skip))
- self._set_data(mdplist)
diff --git a/Docker/scripts/add_missing_stations.py b/Docker/scripts/add_missing_stations.py
deleted file mode 100755
index 2801223f8e03faae7b0bcaf7151453073a05885f..0000000000000000000000000000000000000000
--- a/Docker/scripts/add_missing_stations.py
+++ /dev/null
@@ -1,144 +0,0 @@
-#!/usr/bin/env python
-# -* coding: utf-8 -*-
-
-"""
-Adds phases (=0) and amplitudes (=1) to any missing station if they appear in an h5parm, but not in a particular soltab.
-
-Created on Tue Mar 14 2019
-
-@author: Alexander Drabent
-"""
-
-import argparse
-import numpy as np
-import os
-from losoto.h5parm import h5parm
-from losoto.lib_operations import *
-import logging
-
-def main(h5parmfile, refh5 = None, solset='sol000', refsolset='sol000', soltab_in='phase000', soltab_out='GSMphase', filter='[CR]S*&', bad_antennas='[CR]S*&'):
-
-
- if refh5 == None:
- refh5 = h5parmfile
- pass
-
- if not os.path.exists(h5parmfile) or not os.path.exists(refh5):
- logging.error("H5parm file %s doesn't exist!" % h5parmfile)
- return(1)
-
- if soltab_in == soltab_out:
- logging.error("Output soltab has to be different from input soltab!")
- return(1)
-
- ### Open up the h5parm, get an example value
- data = h5parm(h5parmfile, readonly = False)
- refdata = h5parm(refh5, readonly = False)
- solset = data.getSolset(solset)
- refsolset = refdata.getSolset(refsolset)
-
- ### Get antenna information of the solset
- ref_station_names = sorted(refsolset.getAnt().keys())
- bad_antennas_list = bad_antennas.lstrip(filter).replace('!','').replace('*','').replace('&','').split(';')
- new_station_names = [ ref_station_name for ref_station_name in ref_station_names if ref_station_name not in bad_antennas_list ]
-
- ### Load antenna list of input soltab
- logging.info('Processing solution table %s'%(soltab_in))
- soltab = solset.getSoltab(soltab_in)
- old_station_names = list(soltab.ant[:])
- soltab_type = soltab.getType()
- out_axes_vals = []
- out_axes = soltab.getAxesNames()
- out_axes.insert(1, out_axes.pop(out_axes.index('ant'))) ## put antenna table to the second place
-
- ### resort the val and weights array according to the new axes order
- for vals, weights, coord, selection in soltab.getValuesIter(returnAxes=out_axes, weight=True):
- vals = reorderAxes( vals, soltab.getAxesNames(), out_axes)
- weights = reorderAxes( weights, soltab.getAxesNames(), out_axes )
- pass
-
- ### setting the proper soltab_axis
- for axis in out_axes:
- if axis == 'time':
- out_axes_vals.append(soltab.time)
- elif axis == 'freq':
- out_axes_vals.append(soltab.freq)
- elif axis == 'ant':
- out_axes_vals.append(new_station_names)
- elif axis == 'pol':
- out_axes_vals.append(soltab.pol)
- elif axis == 'dir':
- out_axes_vals.append(soltab.dir)
- else:
- logging.error('Unknown axis in soltab: ' + str(axis))
- data.close()
- refdata.close()
- return 1
-
- ### just copy if number of antennas is the same
- if len(new_station_names) == len(old_station_names):
- logging.warning('Station list in soltab ' + str(soltab_in) + ' matches the station list in the selected solset. Data will just be copied.')
- new_soltab = solset.makeSoltab(soltype=soltab_type, soltabName=soltab_out, axesNames=out_axes, axesVals=out_axes_vals, vals=vals, weights=weights)
-
- ### add missing stations with zero phase/uniform amplitudes and uniform weights
- elif len(new_station_names) > len(old_station_names):
-
- # define proper shape of array
- logging.info('Adding missing antennas to the soltab: ' + str(soltab_out))
- dimension = list(np.shape(vals))
- dimension[1] = len(new_station_names)
-
- if soltab_type == 'amplitude':
- new_vals = np.ones(tuple(dimension))
- else:
- new_vals = np.zeros(tuple(dimension))
- new_weights = np.ones(tuple(dimension))
-
- for i, new_station in enumerate(new_station_names):
- if new_station in old_station_names:
- ant_index = list(soltab.ant).index(new_station)
- new_vals[:,i] = vals[:,ant_index]
- new_weights[:,i] = weights[:,ant_index]
-
- new_soltab = solset.makeSoltab(soltype=soltab_type, soltabName=soltab_out, axesNames=out_axes, axesVals=out_axes_vals, vals=new_vals, weights=new_weights)
-
- else:
- logging.error('There are less antennas in the solset than in the soltab ' + str(soltab_in))
- data.close()
- refdata.close()
- return(1)
-
- data.close()
- refdata.close()
-
- return(0)
-
-if __name__ == "__main__":
- parser = argparse.ArgumentParser(description='Adds phases and amplitudes to any missing station if they appear in an h5parm, but not in a particular soltab.')
-
- parser.add_argument('h5parm', type=str,
- help='H5parm to which this action should be performed .')
- parser.add_argument('--refh5', type=str,
- help='External H5parm from which the full list of antennas is used from.')
- parser.add_argument('--solset', type=str, default='sol000',
- help='Input calibration solutions')
- parser.add_argument('--refsolset', type=str, default='sol000',
- help='Input calibration solutions of the reference h5parm file')
- parser.add_argument('--soltab_in', type=str, default='phase000',
- help='Input solution table')
- parser.add_argument('--soltab_out', type=str, default='GSMphase',
- help='Output solution table (has to be different from input solution table)')
-
-
- args = parser.parse_args()
-
- format_stream = logging.Formatter("%(asctime)s\033[1m %(levelname)s:\033[0m %(message)s","%Y-%m-%d %H:%M:%S")
- format_file = logging.Formatter("%(asctime)s %(levelname)s: %(message)s","%Y-%m-%d %H:%M:%S")
- logging.root.setLevel(logging.INFO)
-
- log = logging.StreamHandler()
- log.setFormatter(format_stream)
- logging.root.addHandler(log)
-
- main(h5parmfile=args.h5parm, refh5=args.refh5, solset=args.solset, refsolset=args.refsolset, soltab_in=args.soltab_in, soltab_out=args.soltab_out)
-
diff --git a/Docker/scripts/check_Ateam_separation.py b/Docker/scripts/check_Ateam_separation.py
deleted file mode 100755
index 51ca65a27eda34bfeffb5c2363fd60a47f460004..0000000000000000000000000000000000000000
--- a/Docker/scripts/check_Ateam_separation.py
+++ /dev/null
@@ -1,182 +0,0 @@
-#!/usr/bin/python
-# -*- coding: utf-8 -*-
-
-#
-# Written by Bas van der Tol (vdtol@strw.leidenuniv.nl), March 2011.
-# Adapted for prefactor by Alexander Drabent (alex@tls-tautenburg.de), February 2019.
-
-#from pylab import *
-import matplotlib
-matplotlib.use('Agg') # Force matplotlib to not use any Xwindows backend.
-import pylab
-import pyrap.quanta as qa
-import pyrap.tables as pt
-import pyrap.measures as pm
-import sys
-import numpy
-import os
-
-targets = [ {'name' : 'CasA', 'ra' : 6.123487680622104, 'dec' : 1.0265153995604648},
- {'name' : 'CygA', 'ra' : 5.233686575770755, 'dec' : 0.7109409582180791},
- {'name' : 'TauA', 'ra' : 1.4596748493730913, 'dec' : 0.38422502335921294},
- {'name' : 'HerA', 'ra' : 4.4119087330382163, 'dec' : 0.087135562905816893},
- {'name' : 'VirA', 'ra' : 3.276086511413598, 'dec' : 0.21626589533567378},
- {'name' : 'Sun'},
- {'name' : 'Jupiter'},
- {'name' : 'Moon'}]
-
-########################################################################
-def input2strlist_nomapfile(invar):
- """
- from bin/download_IONEX.py
- give the list of MSs from the list provided as a string
- """
- str_list = None
- if type(invar) is str:
- if invar.startswith('[') and invar.endswith(']'):
- str_list = [f.strip(' \'\"') for f in invar.strip('[]').split(',')]
- else:
- str_list = [invar.strip(' \'\"')]
- elif type(invar) is list:
- str_list = [str(f).strip(' \'\"') for f in invar]
- else:
- raise TypeError('input2strlist: Type '+str(type(invar))+' unknown!')
- return str_list
-
-########################################################################
-def main(ms_input, min_separation = 30, outputimage = None):
-
- """
- Print seperation of the phase reference center of an input MS
-
-
- Parameters
- ----------
- ms_input : str
- String from the list (map) of the calibrator MSs
-
- Returns
- -------
- 0 --just for printing
- """
-
- msname = input2strlist_nomapfile(ms_input)[0]
-
-
- # Create a measures object
- me = pm.measures()
-
- # Open the measurement set and the antenna and pointing table
- ms = pt.table(msname)
-
- # Get the position of the first antenna and set it as reference frame
- ant_table = pt.table(msname + '::ANTENNA')
- ant_no = 0
- pos = ant_table.getcol('POSITION')
- x = qa.quantity( pos[ant_no,0], 'm' )
- y = qa.quantity( pos[ant_no,1], 'm' )
- z = qa.quantity( pos[ant_no,2], 'm' )
- position = me.position( 'wgs84', x, y, z )
- me.doframe( position )
- ant_table.close()
-
- # Get the first pointing of the first antenna
- field_table = pt.table(msname + '::FIELD')
- field_no = 0
- direction = field_table.getcol('PHASE_DIR')
- ra = direction[ ant_no, field_no, 0 ]
- dec = direction[ ant_no, field_no, 1 ]
- targets.insert(0, {'name' : 'Pointing', 'ra' : ra, 'dec' : dec})
- field_table.close()
-
- # Get a ordered list of unique time stamps from the measurement set
- time_table = pt.taql('select TIME from $1 orderby distinct TIME', tables = [ms])
- time = time_table.getcol('TIME')
- time1 = time/3600.0
- time1 = time1 - pylab.floor(time1[0]/24)*24
-
- ra_qa = qa.quantity( targets[0]['ra'], 'rad' )
- dec_qa = qa.quantity( targets[0]['dec'], 'rad' )
- pointing = me.direction('j2000', ra_qa, dec_qa)
-
- separations = []
-
- print 'SEPERATION from A-Team sources'
- print '------------------------------'
- print 'The minimal accepted distance to an A-Team source is: ' + str(min_separation) + ' deg.'
- for target in targets:
-
- t = qa.quantity(time[0], 's')
- t1 = me.epoch('utc', t)
- me.doframe(t1)
-
- if 'ra' in target.keys():
- ra_qa = qa.quantity( target['ra'], 'rad' )
- dec_qa = qa.quantity( target['dec'], 'rad' )
- direction = me.direction('j2000', ra_qa, dec_qa)
- pass
- else :
- direction = me.direction(target['name'])
- pass
-
- separations.append(me.separation(pointing, direction))
-
- # Loop through all time stamps and calculate the elevation of the pointing
- el = []
- for t in time:
- t_qa = qa.quantity(t, 's')
- t1 = me.epoch('utc', t_qa)
- me.doframe(t1)
- a = me.measure(direction, 'azel')
- elevation = a['m1']
- el.append(elevation['value']/pylab.pi*180)
- pass
-
- el = numpy.array(el)
- pylab.plot(time1, el)
-
- if target['name'] != 'Pointing':
- print target['name'] + ': ' + str(me.separation(pointing, direction))
- if int(float(min_separation)) > int(float(str(me.separation(pointing, direction)).split(' deg')[0])):
- print 'WARNING: The A-Team source ' + target['name'] + ' is closer than ' + str(min_separation) + ' deg to the phase reference center. Calibration might not perform as expected.'
- pass
- pass
-
- pass
- print '------------------------------'
- pylab.title('Pointing Elevation')
- pylab.title('Elevation')
- pylab.ylabel('Elevation (deg)');
- pylab.xlabel('Time (h)');
- pylab.legend( [ target['name'] + '(' + separation.to_string() + ')' for target, separation in zip(targets, separations) ])
-
- if outputimage != None:
- inspection_directory = os.path.dirname(outputimage)
- if not os.path.exists(inspection_directory):
- os.makedirs(inspection_directory)
- print "Directory" , inspection_directory , "created."
- pass
- else:
- print("Directory", inspection_directory, "already exists.")
- pass
- print 'Plotting A-Team elevation to: ' + outputimage
- pylab.savefig(outputimage)
- pass
- return 0
- pass
-
-
-########################################################################
-if __name__ == '__main__':
- import argparse
- parser = argparse.ArgumentParser(description='Print seperation of the phase reference center of an input MS to an A-team source')
-
- parser.add_argument('MSfile', type=str, nargs='+', help='One (or more MSs).')
- parser.add_argument('--min_separation', type=int, default=30, help='minimal accepted distance to an A-team source on the sky in degrees (will raise a WARNING). Default: 30')
- parser.add_argument('--outputimage', type=str, default=None, help='location of the elevation plot of the A-Team sources.')
-
- args = parser.parse_args()
-
- main(args.MSfile, args.min_separation, args.outputimage)
-
- pass
diff --git a/Docker/scripts/check_unflagged_fraction.py b/Docker/scripts/check_unflagged_fraction.py
deleted file mode 100755
index 3b1084a163a3c169c399df7fdd8cde3f1297b31e..0000000000000000000000000000000000000000
--- a/Docker/scripts/check_unflagged_fraction.py
+++ /dev/null
@@ -1,72 +0,0 @@
-#!/usr/bin/env python
-"""
-Flag MS if unflagged fraction is too low.
-"""
-import argparse
-import os
-import numpy as np
-
-def find_unflagged_fraction(ms_file):
- """
- Finds the fraction of data that is unflagged
- Parameters
- ----------
- ms_file : str
- Filename of input MS
- Returns
- -------
- unflagged_fraction : float
- Fraction of unflagged data
- """
- import subprocess
-
- # Call taql. Note that we do not use pt.taql(), as pt.taql() can cause
- # hanging/locking issues on some systems
- p = subprocess.Popen("taql 'CALC sum([select nfalse(FLAG) from {0}]) / "
- "sum([select nelements(FLAG) from {0}])'".format(ms_file),
- shell=True, stdout=subprocess.PIPE)
- r = p.communicate()
- unflagged_fraction = float(r[0])
-
- return unflagged_fraction
-
-def main(ms_file, min_fraction=0.01, print_fraction=False):
- """
- Flag MS if unflagged fraction is too low.
-
- Parameters
- ----------
- ms_file : str
- Name (path) of input MS
- min_fraction : float , optional
- minimum fraction of unflagged data needed to keep this MS
- print_fraction : bool, optional
- print the actual fration of unflagged data
-
- Returns
- -------
- result : dict
- Dict with the name of the input MS or "None"
-
- """
- min_fraction = float(min_fraction)
- unflagged_fraction = find_unflagged_fraction(ms_file)
- if print_fraction:
- print("File %s has %.2f%% unflagged data."%(os.path.basename(ms_file),unflagged_fraction*100.))
- if unflagged_fraction < min_fraction:
- print('check_unflagged_fraction.py: Unflagged fraction of {0} is: {1}, '
- 'removing file.'.format(os.path.basename(ms_file), str(unflagged_fraction)))
- return {'flagged': 'None', 'unflagged_fraction': unflagged_fraction}
- else:
- return {'flagged': ms_file, 'unflagged_fraction': unflagged_fraction}
-
-if __name__ == '__main__':
- descriptiontext = "Check a MS for a minimum fraction of unflagged data.\n"
-
- parser = argparse.ArgumentParser(description=descriptiontext, formatter_class=argparse.RawTextHelpFormatter)
- parser.add_argument('inputms', help='name of the input MS')
- parser.add_argument('-f', '--min_fraction', help='Minimum fraction of unflagged data needed to keep this MS '
- '(default 0.01 = "keep if at least 1%% is not flagged")', type=float, default=0.01)
- args = parser.parse_args()
-
- erg = main(args.inputms, args.min_fraction,print_fraction=True)
diff --git a/Docker/scripts/concat_MS.py b/Docker/scripts/concat_MS.py
deleted file mode 100755
index 4b589b74cdd23855a98c6eb33bd18fdb17a5d9c2..0000000000000000000000000000000000000000
--- a/Docker/scripts/concat_MS.py
+++ /dev/null
@@ -1,132 +0,0 @@
-#!/usr/bin/env python
-import sys
-import glob
-import re
-import pyrap.tables as pt
-import numpy
-import os
-from lofarpipe.support.data_map import DataMap, DataProduct
-
-global_limit = 637871244
-########################################################################
-def input2strlist_nomapfile(invar):
- """
- from bin/download_IONEX.py
- give the list of MSs from the list provided as a string
- """
- str_list = None
- if type(invar) is str:
- if invar.startswith('[') and invar.endswith(']'):
- str_list = [f.strip(' \'\"') for f in invar.strip('[]').split(',')]
- else:
- str_list = [invar.strip(' \'\"')]
- elif type(invar) is list:
- str_list = [str(f).strip(' \'\"') for f in invar]
- else:
- raise TypeError('input2strlist: Type '+str(type(invar))+' unknown!')
- return str_list
-
-########################################################################
-def getsystemmemory():
-
- memory = int(os.popen('cat /proc/meminfo | grep MemAvailable').readlines()[0].split(':')[-1].split()[0])
- return memory
- pass
-
-########################################################################
-def getfilesize(MS):
-
- size = int(os.popen('du -cks ' + MS).readlines()[0].split()[0])
- return size
-
- pass
-
-########################################################################
-def main(ms_input, ms_output, min_length, overhead = 0.8, filename=None, mapfile_dir=None):
-
- """
- Virtually concatenate subbands
-
- Parameters
- ----------
- ms_input : str
- String from the list (map) of the calibrator MSs
- ms_output : str
- String from the outut concatenated MS
- min_length : int
- minimum amount of subbands to concatenate in frequency necessary to perform the wide-band flagging in the RAM. It data is too big aoflag will use indirect-read.
- overhead : float
- Only use this fraction of the available memory for deriving the amount of data to be concatenated.
- filename: str
- Name of output mapfile
- mapfile_dir : str
- Directory for output mapfile
-
- """
- system_memory = getsystemmemory()
- filelist = input2strlist_nomapfile(ms_input)
- file_size = getfilesize(filelist[0])
- overhead = float(overhead)
-
- print "Detected available system memory is: " + str(int(((system_memory / 1024. / 1024.) + 0.5))) + " GB"
- #print "Detected file size is: " + str(int(((file_size / 1024. / 1024.) + 0.5))) + " GB"
- if overhead * system_memory > global_limit:
- system_memory = global_limit
- overhead = 1.0
- print "Number of files to concat will be limited to the global limit of: " + str(int(((global_limit / 1024. / 1024.) + 0.5))) + " GB"
- pass
-
- #i = 1
- max_space = overhead * system_memory / file_size
- max_length = len(filelist) / int((len(filelist) / max_space) + 1.)
-
- #while max_length * file_size > overhead * system_memory:
- #i += 1
- #max_length = len(filelist) / ((len(filelist) / max_space) + i)
- #pass
-
- if max_length >= int(min_length):
- memory = '-memory-read'
- pass
- else:
- memory = '-indirect-read'
- max_length = len(filelist)
- if max_length * file_size > global_limit:
- max_space = global_limit / file_size
- max_length = len(filelist) / int((len(filelist) / max_space) + 1.)
- print "Number of files to concat was limited to the global limit of: " + str(int(((global_limit / 1024. / 1024.) + 0.5))) + " GB"
- print "WARNING: The number of concatenated files will thus be lower than the min_length of: " + str(min_length)
- pass
- pass
-
- print "Applying an overhead of: " + str(overhead)
- print "The max_length value is: " + str(max_length)
- set_ranges = list(numpy.arange(0, len(filelist) + 1, int(max_length)))
- set_ranges[-1] = len(filelist)
-
- map_out = DataMap([])
- for i in numpy.arange(len(set_ranges) - 1):
- f = ms_output + '_' + str(i)
- pt.msconcat(filelist[set_ranges[i]:set_ranges[i + 1]], f)
- map_out.data.append(DataProduct('localhost', f, False))
-
-
-########################################################################
-if __name__ == '__main__':
- import argparse
- parser = argparse.ArgumentParser(description='Virtually concatenate subbands')
-
- parser.add_argument('MSfile', type=str, nargs='+',
- help='One (or more MSs) that we want to concatenate.')
- parser.add_argument('MSout', type=str,
- help='Output MS file')
- parser.add_argument('--min_length', type=str,
- help='Minimum amount of subbands to concatenate in frequency.')
- parser.add_argument('--overhead', type=float, default=0.8,
- help='Only use this fraction of the available memory for deriving the amount of data to be concatenated.')
-
-
-
- args = parser.parse_args()
-
- main(args.MSfile,args.MSout,args.min_length,args.overhead)
diff --git a/Docker/scripts/convert_npys_to_h5parm.py b/Docker/scripts/convert_npys_to_h5parm.py
deleted file mode 100755
index 499684b13aaf61f409c00a7a118fa0adf29de47c..0000000000000000000000000000000000000000
--- a/Docker/scripts/convert_npys_to_h5parm.py
+++ /dev/null
@@ -1,134 +0,0 @@
-#!/usr/bin/env python
-import numpy as np
-import pyrap.tables as pt
-from losoto.h5parm import h5parm
-from losoto.lib_operations import *
-import logging
-from os.path import join
-
-def makesolset(MS, data, solset_name):
- solset = data.makeSolset(solset_name)
-
- logging.info('Collecting information from the ANTENNA table.')
- antennaTable = pt.table(MS + "::ANTENNA", ack=False)
- antennaNames = antennaTable.getcol('NAME')
- antennaPositions = antennaTable.getcol('POSITION')
- antennaTable.close()
- antennaTable = solset.obj._f_get_child('antenna')
- antennaTable.append(zip(*(antennaNames,antennaPositions)))
-
- logging.info('Collecting information from the FIELD table.')
- fieldTable = pt.table(MS + "::FIELD", ack=False)
- phaseDir = fieldTable.getcol('PHASE_DIR')
- pointing = phaseDir[0, 0, :]
- fieldTable.close()
-
- sourceTable = solset.obj._f_get_child('source')
- # add the field centre, that is also the direction for Gain and Common*
- sourceTable.append([('pointing',pointing)])
-
- return antennaNames
- pass
-
-########################################################################
-def main(MSfiles, h5parmfile, store_basename='caldata_transfer', store_directory='.', solset_out = 'calibrator'):
-
- # select on MS
- mslist = MSfiles.lstrip('[').rstrip(']').replace(' ','').replace("'","").split(',')
-
- if len(mslist) == 0:
- logging.error("Did not find any existing directory in input MS list!")
- return(1)
- pass
- else:
- MSfile = mslist[0]
- pass
-
- # Open up the h5parm, get an example value
- data = h5parm(h5parmfile, readonly = False)
-
- # name (path) for parmdb to be written
- reference_station_names = makesolset(MSfile, data, solset_out)
- OutSolset = data.getSolset(solset_out)
-
- # load the numpy arrays written by the previous scripts
- # (filenames constructed in the same way as in these scripts)
- freqs_ampl = np.load(join(store_directory, 'freqs_for_amplitude_array.npy'))
- amps_array = np.load(join(store_directory, store_basename + '_amplitude_array.npy'))
- clock_array = np.load(join(store_directory, 'fitted_data_dclock_' + store_basename + '_1st.npy'))
- tec_array = np.load(join(store_directory, 'fitted_data_dTEC_' + store_basename + '_1st.npy'))
- freqs_phase = np.load(join(store_directory, 'freqs_for_phase_array.npy'))
- phases_array = np.load(join(store_directory, store_basename + '_phase_array.npy'))
- station_names = np.load(join(store_directory, store_basename + '_station_names.npy'))
-
- if len(reference_station_names) > len(station_names):
- missing_stations = set(reference_station_names) - set(station_names)
- logging.warning('CAUTION: For these stations now calibration solutions are available: ' + str(missing_stations))
- pass
-
- # create weights
- dimensions_amp = np.shape(amps_array)
- dimensions_phase = np.shape(phases_array)
- dimensions_clock = np.shape(clock_array)
- dimensions_tec = np.shape(tec_array)
- weights_amp = np.ndarray(shape = dimensions_amp)
- weights_phase = np.ndarray(shape = dimensions_phase)
- weights_clock = np.ndarray(shape = dimensions_clock)
- weights_tec = np.ndarray(shape = dimensions_tec)
- weights_amp.fill(1)
- weights_phase.fill(1)
- weights_clock.fill(1)
- weights_tec.fill(1)
-
- OutSolset.makeSoltab(soltype='clock', soltabName='clock',
- axesNames=['ant'],
- axesVals=[station_names],
- vals=clock_array[-1,:], weights=weights_clock[-1,:])
- OutSolset.makeSoltab(soltype='tec', soltabName='tec',
- axesNames=['ant'],
- axesVals=[station_names],
- vals=tec_array[-1,:], weights=weights_tec[-1,:])
- OutSolset.makeSoltab(soltype='amplitude', soltabName='amplitude',
- axesNames=['ant', 'freq', 'pol'],
- axesVals=[station_names, freqs_ampl, ['XX','YY']],
- vals=amps_array[:,-1,:,:], weights=weights_amp[:,-1,:,:])
- OutSolset.makeSoltab(soltype='phase', soltabName='phase',
- axesNames=['freq', 'ant'],
- axesVals=[freqs_phase, station_names],
- vals=phases_array, weights=weights_phase)
-
- data.close()
- return(0)
-
-
-########################################################################
-if __name__ == '__main__':
- import argparse
- parser = argparse.ArgumentParser(description='Create a parmDB with values from the calibrator.')
-
- parser.add_argument('MSfiles', type=str,
- help='One or more MSs to use as a reference.')
- parser.add_argument('h5parmfile', type=str,
- help='Output H5parm file.')
- parser.add_argument('--basename', type=str, default='caldata_transfer',
- help='Base-name of the numpy files with the calibrator values. (default: \"caldata_transfer\")')
- parser.add_argument('--storedir', type=str, default='.',
- help='Directory of the numpy files with the calibrator values. (default: \".\")')
- parser.add_argument('--solset_out', type=str, default='calibrator',
- help='Name of the solution set')
-
- format_stream = logging.Formatter("%(asctime)s\033[1m %(levelname)s:\033[0m %(message)s","%Y-%m-%d %H:%M:%S")
- format_file = logging.Formatter("%(asctime)s %(levelname)s: %(message)s","%Y-%m-%d %H:%M:%S")
- logging.root.setLevel(logging.INFO)
-
- log = logging.StreamHandler()
- log.setFormatter(format_stream)
- logging.root.addHandler(log)
-
- args = parser.parse_args()
-
- MSfiles = args.MSfiles
- logging.error(MSfiles)
- h5parmfile = args.h5parmfile
-
- main(MSfiles, h5parmfile, store_basename=args.basename, store_directory=args.storedir, solset_out=args.solset_out)
diff --git a/Docker/scripts/createRMh5parm.py b/Docker/scripts/createRMh5parm.py
deleted file mode 100755
index 6e1e5e8c4e20864cb3d762dc8673cbb02d213e4a..0000000000000000000000000000000000000000
--- a/Docker/scripts/createRMh5parm.py
+++ /dev/null
@@ -1,226 +0,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-"""
-Extract RM values for a LOFAR observation. Fill h5parm
-
-Created on Tue Aug 7 11:46:57 2018
-
-@author: mevius
-"""
-from losoto.h5parm import h5parm
-from RMextract import getRM
-from RMextract import PosTools
-import RMextract.getIONEX as ionex
-import pyrap.tables as pt
-import os
-import numpy as np
-import sys
-import logging
-
-
-def makesolset(MS, data, solset_name):
- solset = data.makeSolset(solset_name)
-
- logging.info('Collecting information from the ANTENNA table.')
- antennaTable = pt.table(MS + "::ANTENNA", ack=False)
- antennaNames = antennaTable.getcol('NAME')
- antennaPositions = antennaTable.getcol('POSITION')
- antennaTable.close()
- antennaTable = solset.obj._f_get_child('antenna')
- antennaTable.append(zip(*(antennaNames,antennaPositions)))
-
- logging.info('Collecting information from the FIELD table.')
- fieldTable = pt.table(MS + "::FIELD", ack=False)
- phaseDir = fieldTable.getcol('PHASE_DIR')
- pointing = phaseDir[0, 0, :]
- fieldTable.close()
-
- sourceTable = solset.obj._f_get_child('source')
- # add the field centre, that is also the direction for Gain and Common*
- sourceTable.append([('pointing',pointing)])
-
-
-
-
-def main(MSfiles, h5parmdb, solset_name = "sol000",timestepRM=300,
- ionex_server="ftp://ftp.aiub.unibe.ch/CODE/",
- ionex_prefix='CODG',ionexPath="./",earth_rot=0,proxyServer=None,proxyPort=None,proxyType=None,proxyUser=None,proxyPass=None):
- '''Add rotation measure to existing h5parmdb
-
- Args:
- MSfiles (string) : path + filename of Measurement Set
- h5parmdb (string) : name of existing h5parm
- solset_name (string) : optional name of solset in h5parmdb,
- if not set, first one will be chosen
- timestep (float) : timestep in seconds
- ionex_server (string) : ftp server for IONEX files
- ionex_prefix (string) : prefix of IONEX files
- ionexPath (string) : location where IONEX files are stored
- earth_rot (float) : parameter to determine how much earth rotation is taken \
- into account when interpolating IONEX data. (0 is none, 1 is full)
- '''
-
- try:
- mslist = MSfiles.lstrip('[').rstrip(']').replace(' ','').replace("'","").split(',')
- except AttributeError:
- mslist = MSfiles
-
- if len(mslist) == 0:
- raise ValueError("Did not find any existing directory in input MS list!")
- pass
- else:
- MS = mslist[0]
- pass
-
- if not os.path.exists(h5parmdb):
- logging.error('Could not find h5parmdb: ' + h5parmdb)
- return(1)
-
- data = h5parm(h5parmdb, readonly=False)
- if not solset_name in data.getSolsetNames():
- makesolset(MS,data,solset_name)
- solset = data.getSolset(solset_name)
- soltabs = solset.getSoltabs()
- station_names = sorted(solset.getAnt().keys())
-
- if 'RMextract' in [s.name for s in soltabs]:
- logging.warning('Soltab RMextract exists already. Skipping...')
- return(0)
-
- #extract the timerange information
- (timerange,timestep,pointing,stat_names,stat_pos) = PosTools.getMSinfo(MS)
- start_time = timerange[0]
- end_time = timerange[1]
- timerange[0] = start_time - timestep
- timerange[1] = end_time + timestep
- times,timerange = PosTools.getIONEXtimerange(timerange, timestep)
- if len(times[-1]) == 0 or times[-1][-1] < timerange[1]:
- timestmp = list(times[-1])
- timestmp.append(timerange[1]) #add one extra step to make sure you have a value for all times in the MS in case timestep hase been changed
- times[-1] = np.array(timestmp)
-
- for time_array in times[::-1]: #check in reverse order, since datamight not be available for latest days
- starttime = time_array[0]
- date_parms = PosTools.obtain_observation_year_month_day_fraction(starttime)
- if not proxyServer:
- if not "http" in ionex_server: #ftp server use ftplib
- ionexf = ionex.getIONEXfile(time=date_parms,
- server = ionex_server,
- prefix = ionex_prefix,
- outpath = ionexPath)
- else:
- ionexf = ionex.get_urllib_IONEXfile(time = date_parms,
- server = ionex_server,
- prefix = ionex_prefix,
- outpath = ionexPath)
-
- else:
- ionexf = ionex.get_urllib_IONEXfile(time = date_parms,
- server = ionex_server,
- prefix = ionex_prefix,
- outpath = ionexPath,
- proxy_server = proxyServer,
- proxy_type = proxyType,
- proxy_port = proxyPort,
- proxy_user = proxyUser,
- proxy_pass = proxyPass)
-
- if ionexf == -1:
- if not "igsiono.uwm.edu.pl" in ionex_server:
- logging.info("cannot get IONEX data, try fast product server instead")
- if not proxyServer:
- ionexf = ionex.get_urllib_IONEXfile(time = date_parms,
- server = "https://igsiono.uwm.edu.pl",
- prefix = "igrg",
- outpath = ionexPath)
- else:
- ionexf = ionex.get_urllib_IONEXfile(time = date_parms,
- server = "https://igsiono.uwm.edu.pl",
- prefix = "igrg",
- outpath = ionexPath,
- proxy_server = proxyServer,
- proxy_type = proxyType,
- proxy_port = proxyPort,
- proxy_user = proxyUser,
- proxy_pass = proxyPass)
- if ionexf == -1:
- logging.error("IONEX data not available, even not from fast product server")
- return -1
-
-
- if not proxyServer:
- rmdict = getRM.getRM(MS,
- server = ionex_server,
- prefix = ionex_prefix,
- ionexPath = ionexPath,
- timestep = timestepRM,
- earth_rot = earth_rot)
- else:
- rmdict = getRM.getRM(MS,
- server = ionex_server,
- prefix = ionex_prefix,
- ionexPath = ionexPath,
- timestep = timestepRM,
- earth_rot = earth_rot,
- proxy_server = proxyServer,
- proxy_type = proxyType,
- proxy_port = proxyPort,
- proxy_user = proxyUser,
- proxy_pass = proxyPass)
-
-
- if not rmdict:
- if not ionex_server:
- raise ValueError("One or more IONEX files is not found on disk and download is disabled!\n"
- "(You can run \"bin/download_IONEX.py\" outside the pipeline if needed.)")
- else:
- raise ValueError("Couldn't get RM information from RMextract! (But I don't know why.)")
-
- logging.info('Adding rotation measure values to: ' + solset_name + ' of ' + h5parmdb)
- rm_vals=np.array([rmdict["RM"][stat].flatten() for stat in station_names])
- new_soltab = solset.makeSoltab(soltype='rotationmeasure', soltabName='RMextract',
- axesNames=['ant', 'time'], axesVals=[station_names, rmdict['times']],
- vals=rm_vals,
- weights=np.ones_like(rm_vals, dtype=np.float16))
-
- return(0)
-
- data.close()
-
-
- ########################################################################
-if __name__ == '__main__':
- import argparse
- parser = argparse.ArgumentParser(description='Adds CommonRotationAngle to an H5parm from RMextract.')
-
- parser.add_argument('MSfiles', type=str, nargs='+',
- help='MS for which the parmdb should be created.')
- parser.add_argument('h5parm', type=str,
- help='H5parm to which the results of the CommonRotationAngle is added.')
- parser.add_argument('--server', type=str, default='ftp://ftp.aiub.unibe.ch/CODE/',
- help='URL of the server to use. (default: ftp://ftp.aiub.unibe.ch/CODE/)')
- parser.add_argument('--prefix', type=str, default='CODG',
- help='Prefix of the IONEX files. (default: \"CODG\")')
- parser.add_argument('--ionexpath', '--path', type=str, default='./',
- help='Directory where to store the IONEX files. (default: \"./\")')
- parser.add_argument('--solsetName', '--solset', type=str, default='sol000',
- help='Name of the h5parm solution set (default: sol000)')
- parser.add_argument('-t','--timestep', help=
- 'timestep in seconds. for values <=0 (default) the timegrid of the MS is used ',
- dest="timestep",type=float, default=300.)
- parser.add_argument('-e','--smart_interpol', type=float, default=0, help=
- 'float parameter describing how much of Earth rotation is taken in to account \
- in interpolation of the IONEX files. 1.0 means time interpolation assumes \
- ionosphere rotates in opposite direction of the Earth. 0.0 (default) means \
- no rotation applied',dest="earth_rot",)
-
- args = parser.parse_args()
-
- MS = args.MSfiles
- h5parmdb = args.h5parm
- logging.info("Working on: %s %s" % (MS, h5parmdb))
- main(MS, h5parmdb, ionex_server=args.server, ionex_prefix=args.prefix,
- ionexPath=args.ionexpath, solset_name=args.solsetName,
- timestepRM=args.timestep, earth_rot=args.earth_rot)
-
diff --git a/Docker/scripts/download_skymodel_target.py b/Docker/scripts/download_skymodel_target.py
deleted file mode 100755
index 3475c19c2b5bd791b9d59d41366ca46f9343f62c..0000000000000000000000000000000000000000
--- a/Docker/scripts/download_skymodel_target.py
+++ /dev/null
@@ -1,160 +0,0 @@
-#!/usr/bin/env python
-import os,sys
-import glob
-import pyrap.tables as pt
-import numpy as np
-import time
-import lsmtool
-
-########################################################################
-def grab_coord_MS(MS):
- """
- Read the coordinates of a field from one MS corresponding to the selection given in the parameters
-
- Parameters
- ----------
- MS : str
- Full name (with path) to one MS of the field
-
- Returns
- -------
- RA, Dec : "tuple"
- coordinates of the field (RA, Dec in deg , J2000)
- """
-
- # reading the coordinates ("position") from the MS
- # NB: they are given in rad,rad (J2000)
- [[[ra,dec]]] = pt.table(MS+'::FIELD', readonly=True, ack=False).getcol('PHASE_DIR')
-
- # RA is stocked in the MS in [-pi;pi]
- # => shift for the negative angles before the conversion to deg (so that RA in [0;2pi])
- if ra<0:
- ra=ra+2*np.pi
-
- # convert radians to degrees
- ra_deg = ra/np.pi*180.
- dec_deg = dec/np.pi*180.
-
- # and sending the coordinates in deg
- return ra_deg,dec_deg
-
-
-########################################################################
-def input2strlist_nomapfile(invar):
- """ from bin/download_IONEX.py
- give the list of MSs from the list provided as a string
- """
-
- str_list = None
- if type(invar) is str:
- if invar.startswith('[') and invar.endswith(']'):
- str_list = [f.strip(' \'\"') for f in invar.strip('[]').split(',')]
- else:
- str_list = [invar.strip(' \'\"')]
- elif type(invar) is list:
- str_list = [str(f).strip(' \'\"') for f in invar]
- else:
- raise TypeError('input2strlist: Type '+str(type(invar))+' unknown!')
- return str_list
-
-
-########################################################################
-def main(ms_input, SkymodelPath, Radius="5.", DoDownload="True", Source="TGSS"):
- """
- Download the skymodel for the target field
-
- Parameters
- ----------
- ms_input : str
- String from the list (map) of the target MSs
- SkymodelPath : str
- Full name (with path) to the skymodel; if YES is true, the skymodel will be downloaded here
- Radius : string with float (default = "5.")
- Radius for the TGSS/GSM cone search in degrees
- DoDownload : str ("Force" or "True" or "False")
- Download or not the TGSS skymodel or GSM.
- "Force": download skymodel from TGSS or GSM, delete existing skymodel if needed.
- "True" or "Yes": use existing skymodel file if it exists, download skymodel from
- TGSS or GSM if it does not.
- "False" or "No": Do not download skymodel, raise an exception if skymodel
- file does not exist.
-
- """
-
- FileExists = os.path.isfile(SkymodelPath)
- if (not FileExists and os.path.exists(SkymodelPath)):
- raise ValueError("download_tgss_skymodel_target: Path: \"%s\" exists but is not a file!"%(SkymodelPath))
- download_flag = False
- if not os.path.exists(os.path.dirname(SkymodelPath)):
- os.makedirs(os.path.dirname(SkymodelPath))
- if DoDownload.upper() == "FORCE":
- if FileExists:
- os.remove(SkymodelPath)
- download_flag = True
- elif DoDownload.upper() == "TRUE" or DoDownload.upper() == "YES":
- if FileExists:
- print "USING the exising skymodel in "+ SkymodelPath
- return
- else:
- download_flag = True
- elif DoDownload.upper() == "FALSE" or DoDownload.upper() == "NO":
- if FileExists:
- print "USING the exising skymodel in "+ SkymodelPath
- return
- else:
- raise ValueError("download_tgss_skymodel_target: Path: \"%s\" does not exist and skymodel download is disabled!"%(SkymodelPath))
-
- # If we got here, then we are supposed to download the skymodel.
- assert download_flag is True # Jaja, belts and suspenders...
- print "DOWNLOADING skymodel for the target into "+ SkymodelPath
-
- # Reading a MS to find the coordinate (pyrap)
- [RATar,DECTar]=grab_coord_MS(input2strlist_nomapfile(ms_input)[0])
-
- # Downloading the skymodel, skip after five tries
- errorcode = 1
- tries = 0
- while errorcode != 0 and tries < 5:
- if Source == 'TGSS':
- errorcode = os.system("wget -O "+SkymodelPath+ " \'http://tgssadr.strw.leidenuniv.nl/cgi-bin/gsmv3.cgi?coord="+str(RATar)+","+str(DECTar)+"&radius="+str(Radius)+"&unit=deg&deconv=y\' ")
- pass
- elif Source == 'GSM':
- errorcode = os.system("wget -O "+SkymodelPath+ " \'https://lcs165.lofar.eu/cgi-bin/gsmv1.cgi?coord="+str(RATar)+","+str(DECTar)+"&radius="+str(Radius)+"&unit=deg&deconv=y\' ")
- pass
- time.sleep(5)
- tries += 1
- pass
-
- if not os.path.isfile(SkymodelPath):
- raise IOError("download_tgss_skymodel_target: Path: \"%s\" does not exist after trying to download the skymodel."%(SkymodelPath))
-
- # Treat all sources as one group (direction)
- skymodel = lsmtool.load(SkymodelPath)
- skymodel.group('single')
- skymodel.write(clobber=True)
-
- return
-
-
-########################################################################
-if __name__ == '__main__':
- import argparse
- parser = argparse.ArgumentParser(description=' Download the TGSS or GSM skymodel for the target field')
-
- parser.add_argument('MSfile', type=str, nargs='+',
- help='One (or more MSs) for which a TGSS/GSM skymodel will be download.')
- parser.add_argument('SkyTar', type=str,
- help='Full name (with path) to the skymodel; the TGSS/GSM skymodel will be downloaded here')
- parser.add_argument('--Radius', type=float, default=5.,
- help='Radius for the TGSS/GSM cone search in degrees')
- parser.add_argument('--Source', type=str, default='TGSS',
- help='Choose source for skymodel: TGSS or GSM')
- parser.add_argument('--DoDownload', type=str, default="True",
- help='Download or not the TGSS skymodel or GSM ("Force" or "True" or "False").')
-
- args = parser.parse_args()
- radius=5
- if args.Radius:
- radius=args.Radius
-
- main(args.MSfile, args.SkyTar, str(radius), args.DoDownload, args.Source)
diff --git a/Docker/scripts/find_skymodel_cal.py b/Docker/scripts/find_skymodel_cal.py
deleted file mode 100755
index 9f3ae967dbc63d08f9f4d40ed37cd002d771ba67..0000000000000000000000000000000000000000
--- a/Docker/scripts/find_skymodel_cal.py
+++ /dev/null
@@ -1,161 +0,0 @@
-#!/usr/bin/env python
-import os
-import glob
-import pyrap.tables as pt
-import math
-import lsmtool
-import numpy
-
-def grab_pointing(MS):
- """
- Read the name of the calibrator from one MS corresponding to the selection given in the parameters
-
- Parameters
- ----------
- MS : str
- Full name (with path) to one MS of the field
-
-
- Returns
- -------
- nameCal : str
- Name of the calibrator as provided in the MS
- NB: we suppose that all the calibrators' observations have this field filled in (MS/Observation, column LOFAR_TARGET)
- """
-
- [ra, dec] = pt.table(MS+'::FIELD', readonly=True, ack=False).getcol('PHASE_DIR')[0][0] * 180 / math.pi
- return ra, dec
-
-
-
-########################################################################
-
-def check_skymodel(skymodel, ra, dec, max_separation_arcmin = 1.0):
- """
- Searches for an appropriate sky model
- """
- s = lsmtool.load(skymodel)
- dist_deg = s.getDistance(ra, dec)
- if any(dist_deg * 60.0 < max_separation_arcmin):
- patch_position = int(numpy.where(dist_deg * 60 < max_separation_arcmin)[0][0])
- patch_name = s.getPatchNames()[patch_position]
- return (True, patch_name)
- pass
- else:
- return (False, '')
- pass
- pass
-
-########################################################################
-def find_skymodel(ra, dec, PathSkyMod, extensionSky = ".skymodel", max_separation_arcmin = 1.0):
- """
- Find in the provided folder the correponding skymodel for the given source
-
- Parameters
- ----------
- NameSouce : str
- Name of the source for which we want to find the skymodel
- PathSkyMod : str
- Path to the folder containing the skymodels in use
- (with Pattern, will enable to read and get the info from a MS)
- [extensionSky] : str
- Default: ".skymodel"
- extension of the skymodel files
-
- NB: we suppose that the name of the calibrator is included in the name of the skymodel
-
- Returns
- -------
- list_skymodel[0] : str
- Full name (with path) to the matching skymodel
- """
-
- if os.path.isfile(PathSkyMod):
- print("Checking the skymodel provided: " + PathSkyMod)
- skymodels = [PathSkyMod]
- else:
- skymodels = glob.glob(PathSkyMod + "/*" + extensionSky)
-
- for skymodel in skymodels:
- check = check_skymodel(skymodel, ra, dec, max_separation_arcmin)
- if check[0]:
- print "The following skymodel will be used for the calibrator: " + skymodel.split("/")[-1] + " (in " + PathSkyMod + ")"
- return (skymodel, check[-1])
- pass
- else:
- pass
- pass
-
- raise TypeError('find_skymodel: SKYMODEL FOR THE CALIBRATOR NOT FOUND IN ' + PathSkyMod)
- pass
-
-########################################################################
-def input2strlist_nomapfile(invar):
- """
- from bin/download_IONEX.py
- give the list of MSs from the list provided as a string
- """
- str_list = None
- if type(invar) is str:
- if invar.startswith('[') and invar.endswith(']'):
- str_list = [f.strip(' \'\"') for f in invar.strip('[]').split(',')]
- else:
- str_list = [invar.strip(' \'\"')]
- elif type(invar) is list:
- str_list = [str(f).strip(' \'\"') for f in invar]
- else:
- raise TypeError('input2strlist: Type '+str(type(invar))+' unknown!')
- return str_list
-
-########################################################################
-def main(ms_input, DirSkymodelCal, extensionSky=".skymodel", max_separation_arcmin = 1.0):
-
- """
- Find automatically the skymodel to use for the Calibrator
-
-
- Parameters
- ----------
- ms_input : str
- String from the list (map) of the calibrator MSs
- DirSkymodelCal : str
- Path to the skymodel file, or to the folder where the skymodels are stored
- [extensionSky] : str
- Default: ".skymodel"
- extension of the skymodel files
-
-
- Returns
- -------
- {'SkymodelCal':skymodelCal} : "dict"
- Path to the skymodel of the calibrator
- """
-
- if os.path.isfile(DirSkymodelCal) or os.path.isdir(DirSkymodelCal):
- ra, dec = grab_pointing(input2strlist_nomapfile(ms_input)[0])
- skymodelCal, skymodelName = find_skymodel(ra, dec, DirSkymodelCal, extensionSky, float(max_separation_arcmin))
- return { 'SkymodelCal' : skymodelCal, 'SkymodelName': skymodelName}
- else:
- raise ValueError("find_skymodel_cal: The path \"%s\" is neither a file nor a directory!"%(DirSkymodelCal))
-
-
-########################################################################
-if __name__ == '__main__':
- import argparse
- parser = argparse.ArgumentParser(description='Find automatically between skymodels the one to use (for the Calibrator)')
-
- parser.add_argument('MSfile', type=str, nargs='+',
- help='One (or more MSs) for which we search the matching skymodel.')
- parser.add_argument('DirSky', type=str,
- help='Path to the skymodel file, or to the folder where the skymodels are stored.')
- parser.add_argument('--extsky', type=str,
- help='extension of the skymodel files. (default: \".skymodel\")')
-
- args = parser.parse_args()
- extensionSky='.skymodel'
- if args.extsky:
- extensionSky=args.extsky
-
- main(args.MSfile,args.DirSky, extensionSky)
-
- pass
diff --git a/Docker/scripts/fits2sky.py b/Docker/scripts/fits2sky.py
deleted file mode 100755
index 1359429038d419ceb0d50d88cbcb84653dfa92aa..0000000000000000000000000000000000000000
--- a/Docker/scripts/fits2sky.py
+++ /dev/null
@@ -1,195 +0,0 @@
-#! /usr/bin/env python
-"""
-Script to make a sky model from fits model images
-"""
-import argparse
-from argparse import RawTextHelpFormatter
-from astropy.io import fits
-from astropy import wcs
-import numpy as np
-import scipy.interpolate
-import casacore.tables as pt
-import sys
-import os
-import glob
-
-
-def ra2hhmmss(deg):
- """Convert RA coordinate (in degrees) to HH MM SS"""
-
- from math import modf
- if deg < 0:
- deg += 360.0
- x, hh = modf(deg/15.)
- x, mm = modf(x*60)
- ss = x*60
-
- return (int(hh), int(mm), ss)
-
-
-def dec2ddmmss(deg):
- """Convert DEC coordinate (in degrees) to DD MM SS"""
-
- from math import modf
- sign = (-1 if deg < 0 else 1)
- x, dd = modf(abs(deg))
- x, ma = modf(x*60)
- sa = x*60
-
- return (int(dd), int(ma), sa, sign)
-
-
-def convert_radec_str(ra, dec):
- """Takes ra, dec in degrees and returns makesourcedb strings"""
- ra = ra2hhmmss(ra)
- sra = str(ra[0]).zfill(2)+':'+str(ra[1]).zfill(2)+':'+str("%.3f" % (ra[2])).zfill(6)
- dec = dec2ddmmss(dec)
- decsign = ('-' if dec[3] < 0 else '+')
- sdec = decsign+str(dec[0]).zfill(2)+'.'+str(dec[1]).zfill(2)+'.'+str("%.3f" % (dec[2])).zfill(6)
-
- return sra, sdec
-
-
-def main(fits_models, ms_file, skymodel, fits_masks, min_flux_jy=0.005, interp='linear'):
- """
- Make a makesourcedb sky model for input MS from WSClean fits model images
-
- Parameters
- ----------
- fits_models : str
- Filename of FITS model images. Can be a list of files (e.g.,
- '[mod1.fits,mod2.fits,...]')
- ms_file : str
- Filename of MS for which sky model is to be made. Can be a list of files
- (e.g., '[ms1,ms2,...]', in which case they should all have the same
- frequency
- skymodel : str
- Filename of the output makesourcedb sky model
- fits_masks : str
- Filename of FITS mask images. Can be a list of files (e.g.,
- '[msk1.fits,msk2.fits,...]')
- min_flux_jy : float, optional
- Minimum value of flux in Jy of a source to include in output model
- interp : str, optional
- Interpolation method. Can be any supported by scipy.interpolate.interp1d:
- 'linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic'
-
- """
- min_flux_jy = float(min_flux_jy)
-
- # Get filenames of model images and masks
- if '[' in fits_models and ']' in fits_models:
- fits_models = fits_models.strip('[] ').split(',')
- fits_models = [f.strip('\'\" ') for f in fits_models]
- else:
- fits_models = [fits_models.strip('\'\" ')]
- if '[' in fits_masks and ']' in fits_masks:
- fits_masks = fits_masks.strip('[]').split(',')
- fits_masks = [f.strip('\'\" ') for f in fits_masks]
- else:
- fits_masks = [fits_masks.strip('\'\" ')]
-
- # Read (first) MS file and get the frequency info
- if '[' in ms_file and ']' in ms_file:
- files = ms_file.strip('[] ').split(',')
- #files = [f.strip() for f in files]
- ms_file = files[0].strip('\'\" ')
- sw = pt.table(ms_file+'::SPECTRAL_WINDOW', ack=False)
- ms_freq = sw.col('REF_FREQUENCY')[0]
- ms_freq_low = sw.col('CHAN_FREQ')[0][0]
- ms_freq_high = sw.col('CHAN_FREQ')[0][-1]
- sw.close()
-
- # Get frequencies and data of model images and masks
- freqs = []
- model_images = []
- mask_images = []
- for f in fits_models:
- hdr = fits.getheader(f, 0)
- freqs.append(hdr['CRVAL3']) # Hz
- model_images.append(fits.getdata(f, 0))
- for f in fits_masks:
- mask_images.append(fits.getdata(f, 0))
-
- # Sort by freq
- sorted_ind = np.argsort(freqs)
- freqs = np.array(freqs)[sorted_ind]
- fits_models = np.array(fits_models)[sorted_ind]
- model_images = np.array(model_images)[sorted_ind]
- mask_images = np.array(mask_images)[sorted_ind]
-
- # Check if there is a model at the ms frequency. If so, just use that one
- ind = np.where( np.logical_and(freqs >= ms_freq_low, freqs <= ms_freq_high) )
- if len(ind[0]) == 1:
- freqs = freqs[ind]
- fits_models = fits_models[ind]
- model_images = model_images[ind]
- mask_images = mask_images[ind]
-
- # Set the WCS reference
- hdr = fits.getheader(fits_models[0], 0)
- w = wcs.WCS(hdr)
-
- # Find nonzero pixels in stacked image
- stacked_model = np.zeros(model_images[0].shape)
- stacked_mask = np.zeros(mask_images[0].shape)
- for im, ma in zip(model_images, mask_images):
- stacked_model += im
- stacked_mask += ma
- nonzero_ind = np.where((stacked_model != 0.0) & (stacked_mask > 0))
-
- # Interpolate the fluxes to the frequency of the MS
- nsources = len(nonzero_ind[0])
- fluxes = []
- names = []
- ras = []
- decs = []
- for i in range(nsources):
- flux_array = []
- freq_array = []
- for ind, (im, ma) in enumerate(zip(model_images, mask_images)):
- index = [nonzero_ind[j][i] for j in range(4)]
- if ma[tuple(index)] > 0:
- flux_array.append(im[tuple(index)])
- freq_array.append(freqs[ind])
-
- if len(flux_array) > 0:
- # Only create a source entry if MS frequency is within +/- 4 MHz of
- # the sampled frequency range to prevent large extrapolations when
- # the primary beam may cause rapid spectral changes
- if ms_freq > freq_array[0]-4e6 and ms_freq < freq_array[-1]+4e6:
- # If MS frequency lies outside range, just use nearest freq
- if ms_freq < freq_array[0]:
- flux = flux_array[0]
- elif ms_freq > freq_array[-1]:
- flux = flux_array[-1]
- else:
- # Otherwise interpolate
- flux = scipy.interpolate.interp1d(freq_array, flux_array, kind=interp)(ms_freq)
- if flux > min_flux_jy:
- index.reverse() # change to WCS coords
- ras.append(w.wcs_pix2world(np.array([index]), 0, ra_dec_order=True)[0][0])
- decs.append(w.wcs_pix2world(np.array([index]), 0, ra_dec_order=True)[0][1])
- names.append('cc{}'.format(i))
- fluxes.append(flux)
-
- # Write sky model
- with open(skymodel, 'w') as outfile:
- outfile.write('FORMAT = Name, Type, Ra, Dec, I, Q, U, V, ReferenceFrequency\n')
- for name, ra, dec, flux in zip(names, ras, decs, fluxes):
- ra_str, dec_str = convert_radec_str(ra, dec)
- outfile.write('{0}, POINT, {1}, {2}, {3}, 0.0, 0.0, 0.0, {4}\n'
- .format(name, ra_str, dec_str, flux, ms_freq))
-
-
-if __name__ == '__main__':
- descriptiontext = "Make a makesourcedb sky model from WSClean fits model images.\n"
-
- parser = argparse.ArgumentParser(description=descriptiontext, formatter_class=RawTextHelpFormatter)
- parser.add_argument('fits_models', help='Model images')
- parser.add_argument('msfile', help='MS file')
- parser.add_argument('skymodel', help='Filename of output sky model')
- parser.add_argument('fits_masks', help='Mask images')
-
- args = parser.parse_args()
- main(args.fits_models, args.msfile, args.skymodel, args.fits_masks)
diff --git a/Docker/scripts/getStructure_from_phases.py b/Docker/scripts/getStructure_from_phases.py
deleted file mode 100755
index fcfe1b5b4081bc468681b1a753082c6201c038ef..0000000000000000000000000000000000000000
--- a/Docker/scripts/getStructure_from_phases.py
+++ /dev/null
@@ -1,226 +0,0 @@
-#!/usr/bin/env python
-
-##########################################################################
-#
-# Written my Maaijke Mevius -- see Radio Science publication DOI: 10.1002/2016RS006028
-#
-##########################################################################
-# History
-# 20/07/2016 Outlier rejection and argparse added -- Tim Shimwell
-# 06/06/2019 Changed for the use with h5parms -- Alexander Drabent
-
-import matplotlib as mpl
-mpl.use("Agg")
-from losoto.h5parm import h5parm
-from losoto.lib_operations import *
-from pylab import *
-import scipy
-from scipy.optimize import leastsq
-import argparse
-
-def mad(arr):
- """ Median Absolute Deviation: a "Robust" version of standard deviation.
- Indices variabililty of the sample.
- https://en.wikipedia.org/wiki/Median_absolute_deviation
- """
- arr = np.ma.array(arr).compressed() # should be faster to not use masked arrays.
- med = np.median(arr)
- return np.median(np.abs(arr - med))
-
-def model(t, coeffs):
- coeffs[0] = abs(coeffs[0])
- return coeffs[0] + (t**coeffs[2])/coeffs[1]
-def residuals(coeffs, y, t):
- coeffs[0] = abs(coeffs[0])
- return y - model(t, coeffs)
-
-def main(h5parmfile,solset='sol000',soltab='phase000',nr_grid=1,doplot=True,outbasename='ionosphere',output_dir='./',dofilter=True):
- outfile = open(str(output_dir) + '/' + outbasename +'_structure.txt','w')
- data = h5parm(h5parmfile)
- solset = data.getSolset(solset)
- soltab = solset.getSoltab(soltab)
- soltype = soltab.getType()
- if soltype != 'phase':
- logging.warning("Soltab is of type " + soltype + ", but should be phase. Skipping.")
- data.close()
- return(0)
- vals = soltab.val
- station_names = soltab.ant[:]
- stations = solset.getAnt()
- phx=[]
- phy=[]
- allposx=[]
- allposy=[]
- nrtimes = len(soltab.time)
- nrfreqss = len(soltab.freq)
- timestep=int(nrtimes/nr_grid)
- for vals, coord, selection in soltab.getValuesIter(returnAxes=soltab.getAxesNames(), weight=False):
- try:
- vals = reorderAxes( vals, soltab.getAxesNames(), ['pol', 'ant', 'time', 'freq', 'dir'])
- vals = vals[:,:,:,:,0]
- except:
- vals = reorderAxes( vals, soltab.getAxesNames(), ['pol', 'ant', 'time', 'freq'])
-
- valx = vals[0]
- valy = vals[1]
- for i,station_name in enumerate(station_names):
- if not "CS" in station_name:
- continue
- phx.append(valx[i,:,0].flatten())
- allposx.append(stations[station_name])
- phy.append(valy[i,:,0].flatten())
- allposy.append(stations[station_name])
-
- phx=np.array(phx)
- phy=np.array(phy)
- allposx=np.array(allposx)
- allposy=np.array(allposy)
- D=allposx[:,np.newaxis,:]-allposx[np.newaxis,:]
- D2=np.sqrt(np.sum(D**2,axis=-1))
- Dy=allposy[:,np.newaxis,:]-allposy[np.newaxis,:]
- D2y=np.sqrt(np.sum(Dy**2,axis=-1))
- S0s=[]
- betas=[]
- for itime in xrange(nr_grid+1):
- tm=[0,1000000000]
- if itime<nr_grid:
- tm[0]=int(itime*timestep)
- tm[1]=int(tm[0]+timestep)
- dphx=phx[:,np.newaxis,tm[0]:tm[1]]-phx[np.newaxis,:,tm[0]:tm[1]]
- dphy=phy[:,np.newaxis,tm[0]:tm[1]]-phy[np.newaxis,:,tm[0]:tm[1]]
- dphx=np.unwrap(np.remainder(dphx-dphx[:,:,0][:,:,np.newaxis]+np.pi,2*np.pi))
- dvarx=np.var(dphx,axis=-1)
- dphy=np.unwrap(np.remainder(dphy-dphy[:,:,0][:,:,np.newaxis]+np.pi,2*np.pi))
- dvary=np.var(dphy,axis=-1)
- myselect=np.logical_and(D2>0,np.logical_and(np.any(np.logical_and(dvarx>1e-7,dvarx<.1),axis=0)[np.newaxis],np.any(np.logical_and(dvarx>1e-7,dvarx<.1),axis=1)[:,np.newaxis]))
- if dofilter:
- x=D2[myselect]
- y=dvarx[myselect]
- # Seems like real values dont occur above 1.0
- flagselect = np.where(y > 1.0)
- xplot = np.delete(x,flagselect)
- yplot = np.delete(y,flagselect)
-
-
- bins = np.logspace(np.log10(np.min(xplot)),np.log10(np.max(xplot)),10)
- binys = []
- binxs = []
- for i in range(0,len(bins)-1):
- binmin = bins[i]
- binmax = bins[i+1]
- inbin = np.intersect1d(np.where(xplot > binmin),np.where(xplot < binmax))
- binxs.append((binmin+binmax)/2.0)
- binys.append(np.median(yplot[inbin])+10*mad(yplot[inbin]))
- x0 = [0.1,1.0,1.0]
- xfit, flag = scipy.optimize.leastsq(residuals, x0, args=(binys,binxs))
- flagselect = np.where(y > model(x,xfit))
- print xfit,'fitting'
- if doplot:
- x=D2[myselect]
- y=dvarx[myselect]
- subplot(3,1,1)
- scatter(x,y,color='b')
- if dofilter:
- y2 = y[flagselect]
- x2 = x[flagselect]
- scatter(x2,y2,color='r')
- scatter(x,model(x,xfit),color='gray',linestyle='-')#,'gray',linestyle='-',linewidth=3)
- x=D2[np.logical_and(D2>1e3,myselect)]
- y=dvarx[np.logical_and(D2>1e3,myselect)]
- if dofilter:
- flagselect = np.where(y > model(x,xfit))
- x = np.delete(x,flagselect)
- y = np.delete(y,flagselect)
- A=np.ones((2,x.shape[0]),dtype=float)
- A[1,:]=np.log10(x)
- par=np.dot(np.linalg.inv(np.dot(A,A.T)),np.dot(A,np.log10(y)))
- S0=10**(-1*np.array(par)[0]/np.array(par)[1])
- if doplot:
- plot(x,10**(par[0]+par[1]*np.log10(x)),'r-')
- xscale("log")
- yscale("log")
- xlim(30,4000)
- xlabel('baseline length [m]')
- ylabel('XX phase variance [rad$^2$]')
- title('XX diffractive scale: %3.1f km'%(float(S0)/1000.))
- myselect=np.logical_and(D2y>0,np.logical_and(np.any(np.logical_and(dvary>1e-7,dvary<.1),axis=0)[np.newaxis],np.any(np.logical_and(dvary>1e-7,dvary<.1),axis=1)[:,np.newaxis]))
- if dofilter:
- x=D2y[myselect]
- y=dvary[myselect]
- # seems like real values dont occur above 1.0
- flagselect = np.where(y > 1.0)
- xplot = np.delete(x,flagselect)
- yplot = np.delete(y,flagselect)
- bins = np.logspace(np.log10(np.min(xplot)),np.log10(np.max(xplot)),20)
- binys = []
- binxs = []
- for i in range(0,len(bins)-1):
- binmin = bins[i]
- binmax = bins[i+1]
- inbin = np.intersect1d(np.where(xplot > binmin),np.where(xplot < binmax))
- binxs.append((binmin+binmax)/2.0)
- binys.append(np.median(yplot[inbin])+10*mad(yplot[inbin]))
- x0 = [0.1,1.0,1.0]
- xfit, flag = scipy.optimize.leastsq(residuals, x0, args=(binys,binxs))
- flagselect = np.where(y > model(x,xfit))
- print xfit,'fitting'
- if doplot:
- x=D2y[myselect]
- y=dvary[myselect]
- subplot(3,1,3)
- scatter(x,y,color='b')
- if dofilter:
- y2 = y[flagselect]
- x2 = x[flagselect]
- scatter(x2,y2,color='r')
- scatter(x,model(x,xfit),color='gray',linestyle='-')#,'gray',linestyle='-',linewidth=3)
- x=D2y[np.logical_and(D2y>1e3,myselect)]
- y=dvary[np.logical_and(D2y>1e3,myselect)]
- if dofilter:
- flagselect = np.where(y > model(x,xfit))
- x = np.delete(x,flagselect)
- y = np.delete(y,flagselect)
- A=np.ones((2,x.shape[0]),dtype=float)
- A[1,:]=np.log10(x)
- pary=np.dot(np.linalg.inv(np.dot(A,A.T)),np.dot(A,np.log10(y)))
- S0y=10**(-1*np.array(pary)[0]/np.array(pary)[1])
- if doplot:
- plot(x,10**(pary[0]+pary[1]*np.log10(x)),'r-')
- xscale("log")
- yscale("log")
- xlim(30,4000)
- xlabel('baseline length [m]')
- ylabel('YY phase variance [rad$^2$]')
- title('YY diffractive scale: %3.1f km'%(float(S0y)/1000.))
- savefig(output_dir + '/' + outbasename + '_structure.png')
- close()
- cla()
- S0s.append([S0,S0y])
- betas.append([par[1],pary[1]])
- outfile.write('S0s ****%s**** %s beta %s %s\n'%(S0,S0y,par[1],pary[1]))
- data.close()
- return(0)
-
-
-if __name__ == "__main__":
- parser = argparse.ArgumentParser()
- parser.add_argument('h5parmfile',help='Name of input h5parm file')
- parser.add_argument('--solset',help='Name of solution set to use')
- parser.add_argument('--soltab',help='Name of solution table to use')
- parser.add_argument('--outbasename',help='Outfile base name')
- parser.add_argument('--output_dir', help='Output directory', default='./')
- parser.add_argument('-p',dest='doplot',default=True,help='Do plotting')
- parser.add_argument('-f',dest='dofilter',default=True,help='Attempt to filter out odd solutions')
- parser.add_argument('-nr_grid',dest='nr_grid',default=1,help='Time step')
-
- args = parser.parse_args()
- h5parmfile = args.h5parmfile
- solset = args.solset
- soltab = args.soltab
- dofilter = args.dofilter
- doplot = args.doplot
- outbasename = args.outbasename
- output_dir = args.output_dir
- nr_grid = args.nr_grid
-
- main(h5parmfile,solset,soltab,nr_grid,doplot,outbasename,output_dir,dofilter)
diff --git a/Docker/scripts/h5parm_pointingname.py b/Docker/scripts/h5parm_pointingname.py
deleted file mode 100755
index 60d7a879b9541ad280d0ba872abc009f3ec1e630..0000000000000000000000000000000000000000
--- a/Docker/scripts/h5parm_pointingname.py
+++ /dev/null
@@ -1,41 +0,0 @@
-#!/usr/bin/env python
-
-from losoto.h5parm import h5parm
-from losoto.lib_operations import *
-
-import numpy
-########################################################################
-def main(h5parmdb, solsetName='sol000', pointing='POINTING'):
-
-
- data = h5parm(h5parmdb, readonly = False)
- solset = data.getSolset(solsetName)
-
- sources = solset.obj._f_get_child('source')
- direction = list(sources[0][-1])
-
- for i in numpy.arange(len(sources)):
- sources[i] = (pointing, direction)
- pass
-
- data.close()
- return 0
-
-
-########################################################################
-if __name__ == '__main__':
- import argparse
- parser = argparse.ArgumentParser(description='Define the pointing direction of an h5parm.')
-
- parser.add_argument('h5parmdb', type=str,
- help='H5parm of which pointing will be set.')
- parser.add_argument('--solsetName', '--solset', type=str, default='sol000',
- help='Name of the h5parm solution set (default: sol000)')
- parser.add_argument('--pointing', type=str, default='POINTING',
- help='Name of the h5parm solution set (default: POINTING)')
-
- args = parser.parse_args()
-
- h5parmdb = args.h5parmdb
- logging.info("Working on: %s" % (h5parmdb))
- main(h5parmdb, solsetName=args.solsetName, pointing=args.pointing)
diff --git a/Docker/scripts/make_clean_mask.py b/Docker/scripts/make_clean_mask.py
deleted file mode 100755
index b504be4a1ff077cc21f06e527fb24f909c7b0074..0000000000000000000000000000000000000000
--- a/Docker/scripts/make_clean_mask.py
+++ /dev/null
@@ -1,793 +0,0 @@
-#! /usr/bin/env python
-"""
-Script to make a clean mask
-"""
-import argparse
-from argparse import RawTextHelpFormatter
-import casacore.images as pim
-from astropy.io import fits as pyfits
-from astropy.coordinates import Angle
-import pickle
-import numpy as np
-import sys
-import os
-
-# I copy&pasted that into this file
-#from factor.lib.polygon import Polygon
-
-# workaround for a bug / ugly behavior in matplotlib / pybdsf
-# (some installtaions of matplotlib.pyplot fail if there is no X-Server available
-# which is not caught in pybdsf.)
-try:
- import matplotlib
- matplotlib.use('Agg')
-except (RuntimeError, ImportError):
- pass
-try:
- import bdsf
-except ImportError:
- from lofar import bdsm as bdsf
-
-class Polygon:
- """
- Generic polygon class
-
- Polygons are used to define the facet boundaries.
-
- Parameters
- ----------
- x : array
- A sequence of nodal x-coords.
- y : array
- A sequence of nodal y-coords.
-
- """
- def __init__(self, x, y):
- if len(x) != len(y):
- raise IndexError('x and y must be equally sized.')
- self.x = np.asfarray(x)
- self.y = np.asfarray(y)
-
- # Closes the polygon if were open
- x1, y1 = x[0], y[0]
- xn, yn = x[-1], y[-1]
- if x1 != xn or y1 != yn:
- self.x = np.concatenate((self.x, [x1]))
- self.y = np.concatenate((self.y, [y1]))
-
- # Anti-clockwise coordinates
- if _det(self.x, self.y) < 0:
- self.x = self.x[::-1]
- self.y = self.y[::-1]
-
-
- def is_inside(self, xpoint, ypoint, smalld=1e-12):
- """
- Check if point is inside a general polygon.
-
- An improved version of the algorithm of Nordbeck and Rydstedt.
-
- REF: SLOAN, S.W. (1985): A point-in-polygon program. Adv. Eng.
- Software, Vol 7, No. 1, pp 45-47.
-
- Parameters
- ----------
- xpoint : array or float
- The x-coord of the point to be tested.
- ypoint : array or float
- The y-coords of the point to be tested.
- smalld : float
- Tolerance within which point is considered to be on a side.
-
- Returns
- -------
- mindst : array or float
- The distance from the point to the nearest point of the polygon:
-
- If mindst < 0 then point is outside the polygon.
- If mindst = 0 then point in on a side of the polygon.
- If mindst > 0 then point is inside the polygon.
-
- """
- xpoint = np.asfarray(xpoint)
- ypoint = np.asfarray(ypoint)
-
- # Scalar to array
- if xpoint.shape is tuple():
- xpoint = np.array([xpoint], dtype=float)
- ypoint = np.array([ypoint], dtype=float)
- scalar = True
- else:
- scalar = False
- # Check consistency
- if xpoint.shape != ypoint.shape:
- raise IndexError('x and y must be equally sized.')
-
- # If snear = True: Dist to nearest side < nearest vertex
- # If snear = False: Dist to nearest vertex < nearest side
- snear = np.ma.masked_all(xpoint.shape, dtype=bool)
-
- # Initialize arrays
- mindst = np.ones_like(xpoint, dtype=float) * np.inf
- j = np.ma.masked_all(xpoint.shape, dtype=int)
- x = self.x
- y = self.y
- n = len(x) - 1 # Number of sides/vertices defining the polygon
-
- # Loop over each side defining polygon
- for i in range(n):
- d = np.ones_like(xpoint, dtype=float) * np.inf
-
- # Start of side has coords (x1, y1)
- # End of side has coords (x2, y2)
- # Point has coords (xpoint, ypoint)
- x1 = x[i]
- y1 = y[i]
- x21 = x[i + 1] - x1
- y21 = y[i + 1] - y1
- x1p = x1 - xpoint
- y1p = y1 - ypoint
-
- # Points on infinite line defined by
- # x = x1 + t * (x1 - x2)
- # y = y1 + t * (y1 - y2)
- # where
- # t = 0 at (x1, y1)
- # t = 1 at (x2, y2)
- # Find where normal passing through (xpoint, ypoint) intersects
- # infinite line
- t = -(x1p * x21 + y1p * y21) / (x21 ** 2 + y21 ** 2)
- tlt0 = t < 0
- tle1 = (0 <= t) & (t <= 1)
-
- # Normal intersects side
- d[tle1] = ((x1p[tle1] + t[tle1] * x21) ** 2 +
- (y1p[tle1] + t[tle1] * y21) ** 2)
-
- # Normal does not intersects side
- # Point is closest to vertex (x1, y1)
- # Compute square of distance to this vertex
- d[tlt0] = x1p[tlt0] ** 2 + y1p[tlt0] ** 2
-
- # Store distances
- mask = d < mindst
- mindst[mask] = d[mask]
- j[mask] = i
-
- # Point is closer to (x1, y1) than any other vertex or side
- snear[mask & tlt0] = False
-
- # Point is closer to this side than to any other side or vertex
- snear[mask & tle1] = True
-
- if np.ma.count(snear) != snear.size:
- raise IndexError('Error computing distances')
- mindst **= 0.5
-
- # Point is closer to its nearest vertex than its nearest side, check if
- # nearest vertex is concave.
- # If the nearest vertex is concave then point is inside the polygon,
- # else the point is outside the polygon.
- jo = j.copy()
- jo[j == 0] -= 1
- area = _det([x[j + 1], x[j], x[jo - 1]], [y[j + 1], y[j], y[jo - 1]])
- mindst[~snear] = np.copysign(mindst, area)[~snear]
-
- # Point is closer to its nearest side than to its nearest vertex, check
- # if point is to left or right of this side.
- # If point is to left of side it is inside polygon, else point is
- # outside polygon.
- area = _det([x[j], x[j + 1], xpoint], [y[j], y[j + 1], ypoint])
- mindst[snear] = np.copysign(mindst, area)[snear]
-
- # Point is on side of polygon
- mindst[np.fabs(mindst) < smalld] = 0
-
- # If input values were scalar then the output should be too
- if scalar:
- mindst = float(mindst)
- return mindst
-
-
-
-def read_vertices(filename):
- """
- Returns facet vertices stored in input file
- """
- with open(filename, 'r') as f:
- direction_dict = pickle.load(f)
- return direction_dict['vertices']
-
-
-def read_casa_polys(filename, image):
- """
- Reads casa region file and returns polys
-
- Note: only regions of type "poly" are supported
- """
- with open(filename, 'r') as f:
- lines = f.readlines()
-
- polys = []
- for line in lines:
- if line.startswith('poly'):
- poly_str_temp = line.split('[[')[1]
- poly_str = poly_str_temp.split(']]')[0]
- poly_str_list = poly_str.split('], [')
- ra = []
- dec = []
- for pos in poly_str_list:
- RAstr, Decstr = pos.split(',')
- ra.append(Angle(RAstr, unit='hourangle').to('deg').value)
- dec.append(Angle(Decstr.replace('.', ':', 2), unit='deg').to('deg').value)
- poly_vertices = [np.array(ra), np.array(dec)]
-
- # Convert to image-plane polygon
- xvert = []
- yvert = []
- for RAvert, Decvert in zip(np.array(ra), np.array(dec)):
- try:
- pixels = image.topixel([0, 1, Decvert*np.pi/180.0,
- RAvert*np.pi/180.0])
- except:
- pixels = image.topixel([1, 1, Decvert*np.pi/180.0,
- RAvert*np.pi/180.0])
- xvert.append(pixels[2]) # x -> Dec
- yvert.append(pixels[3]) # y -> RA
- polys.append(Polygon(xvert, yvert))
-
- elif line.startswith('ellipse'):
- ell_str_temp = line.split('[[')[1]
- if '], 0.0' not in ell_str_temp and '], 90.0' not in ell_str_temp:
- print('Only position angles of 0.0 and 90.0 are support for CASA '
- 'regions of type "ellipse"')
- sys.exit(1)
- if '], 0.0' in ell_str_temp:
- ell_str = ell_str_temp.split('], 0.0')[0]
- pa = 0
- else:
- ell_str = ell_str_temp.split('], 90.0')[0]
- pa = 90
- ell_str_list = ell_str.split('], [')
-
- # Ellipse center
- RAstr, Decstr = ell_str_list[0].split(',')
- ra_center = Angle(RAstr, unit='hourangle').to('deg').value
- dec_center = Angle(Decstr.replace('.', ':', 2), unit='deg').to('deg').value
- pixels = image.topixel([0, 1, dec_center*np.pi/180.0,
- ra_center*np.pi/180.0])
- x_center = pixels[2] # x -> Dec
- y_center = pixels[3] # y -> RA
-
- # Ellipse semimajor and semiminor axes
- a_str, b_str = ell_str_list[1].split(',')
- a_deg = float(a_str.split('arcsec')[0])/3600.0
- b_deg = float(b_str.split('arcsec')[0])/3600.0
- pixels1 = image.topixel([0, 1, (dec_center-a_deg/2.0)*np.pi/180.0,
- ra_center*np.pi/180.0])
- a_pix1 = pixels1[2]
- pixels2 = image.topixel([0, 1, (dec_center+a_deg/2.0)*np.pi/180.0,
- ra_center*np.pi/180.0])
- a_pix2 = pixels2[2]
- a_pix = abs(a_pix2 - a_pix1)
- ex = []
- ey = []
- for th in range(0, 360, 1):
- if pa == 0:
- # semimajor axis is along x-axis
- ex.append(a_pix * np.cos(th * np.pi / 180.0)
- + x_center) # x -> Dec
- ey.append(a_pix * b_deg / a_deg * np.sin(th * np.pi / 180.0) + y_center) # y -> RA
- elif pa == 90:
- # semimajor axis is along y-axis
- ex.append(a_pix * b_deg / a_deg * np.cos(th * np.pi / 180.0)
- + x_center) # x -> Dec
- ey.append(a_pix * np.sin(th * np.pi / 180.0) + y_center) # y -> RA
- polys.append(Polygon(ex, ey))
-
- elif line.startswith('box'):
- poly_str_temp = line.split('[[')[1]
- poly_str = poly_str_temp.split(']]')[0]
- poly_str_list = poly_str.split('], [')
- ra = []
- dec = []
- for pos in poly_str_list:
- RAstr, Decstr = pos.split(',')
- ra.append(Angle(RAstr, unit='hourangle').to('deg').value)
- dec.append(Angle(Decstr.replace('.', ':', 2), unit='deg').to('deg').value)
- ra.insert(1, ra[0])
- dec.insert(1, dec[1])
- ra.append(ra[2])
- dec.append(dec[0])
- poly_vertices = [np.array(ra), np.array(dec)]
-
- # Convert to image-plane polygon
- xvert = []
- yvert = []
- for RAvert, Decvert in zip(np.array(ra), np.array(dec)):
- try:
- pixels = image.topixel([0, 1, Decvert*np.pi/180.0,
- RAvert*np.pi/180.0])
- except:
- pixels = image.topixel([1, 1, Decvert*np.pi/180.0,
- RAvert*np.pi/180.0])
- xvert.append(pixels[2]) # x -> Dec
- yvert.append(pixels[3]) # y -> RA
- polys.append(Polygon(xvert, yvert))
-
- elif line.startswith('#'):
- pass
-
- else:
- print('Only CASA regions of type "poly", "box", or "ellipse" are supported')
- sys.exit(1)
-
- return polys
-
-
-def make_template_image(image_name, reference_ra_deg, reference_dec_deg,
- imsize=512, cellsize_deg=0.000417):
- """
- Make a blank image and save it to disk
-
- Parameters
- ----------
- image_name : str
- Filename of output image
- reference_ra_deg : float, optional
- RA for center of output mask image
- reference_dec_deg : float, optional
- Dec for center of output mask image
- imsize : int, optional
- Size of output image
- cellsize_deg : float, optional
- Size of a pixel in degrees
-
- """
- shape_out = [1, 1, imsize, imsize]
- hdu = pyfits.PrimaryHDU(np.zeros(shape_out, dtype=np.float32))
- hdulist = pyfits.HDUList([hdu])
- header = hdulist[0].header
-
- # Add WCS info
- header['CRVAL1'] = reference_ra_deg
- header['CDELT1'] = -cellsize_deg
- header['CRPIX1'] = imsize/2.0
- header['CUNIT1'] = 'deg'
- header['CTYPE1'] = 'RA---SIN'
- header['CRVAL2'] = reference_dec_deg
- header['CDELT2'] = cellsize_deg
- header['CRPIX2'] = imsize/2.0
- header['CUNIT2'] = 'deg'
- header['CTYPE2'] = 'DEC--SIN'
-
- # Add STOKES info
- header['CRVAL3'] = 1.0
- header['CDELT3'] = 1.0
- header['CRPIX3'] = 1.0
- header['CUNIT3'] = ''
- header['CTYPE3'] = 'STOKES'
-
- # Add frequency info
- header['RESTFRQ'] = 15036
- header['CRVAL4'] = 150e6
- header['CDELT4'] = 3e8
- header['CRPIX4'] = 1.0
- header['CUNIT4'] = 'HZ'
- header['CTYPE4'] = 'FREQ'
- header['SPECSYS'] = 'TOPOCENT'
-
- # Add equinox
- header['EQUINOX'] = 2000.0
-
- # Add telescope
- header['TELESCOP'] = 'LOFAR'
-
- hdulist[0].header = header
-
- hdulist.writeto(image_name, clobber=True)
- hdulist.close()
-
-
-
-def main(image_name, mask_name, atrous_do=False, threshisl=0.0, threshpix=0.0, rmsbox=None,
- rmsbox_bright=(35, 7), iterate_threshold=False, adaptive_rmsbox=False, img_format='fits',
- threshold_format='float', trim_by=0.0, vertices_file=None, atrous_jmax=6,
- pad_to_size=None, skip_source_detection=False, region_file=None, nsig=1.0,
- reference_ra_deg=None, reference_dec_deg=None, cellsize_deg=0.000417,
- use_adaptive_threshold=False, adaptive_thresh=150.0):
- """
- Make a clean mask and return clean threshold
-
- Parameters
- ----------
- image_name : str
- Filename of input image from which mask will be made. If the image does
- not exist, a template image with center at (reference_ra_deg,
- reference_dec_deg) will be made internally
- mask_name : str
- Filename of output mask image
- atrous_do : bool, optional
- Use wavelet module of PyBDSF?
- threshisl : float, optional
- Value of thresh_isl PyBDSF parameter
- threshpix : float, optional
- Value of thresh_pix PyBDSF parameter
- rmsbox : tuple of floats, optional
- Value of rms_box PyBDSF parameter
- rmsbox_bright : tuple of floats, optional
- Value of rms_box_bright PyBDSF parameter
- iterate_threshold : bool, optional
- If True, threshold will be lower in 20% steps until
- at least one island is found
- adaptive_rmsbox : tuple of floats, optional
- Value of adaptive_rms_box PyBDSF parameter
- img_format : str, optional
- Format of output mask image (one of 'fits' or 'casa')
- threshold_format : str, optional
- Format of output threshold (one of 'float' or 'str_with_units')
- trim_by : float, optional
- Fraction by which the perimeter of the output mask will be
- trimmed (zeroed)
- vertices_file : str, optional
- Filename of file with vertices (must be a pickle file containing
- a dictionary with the vertices in the 'vertices' entry)
- atrous_jmax : int, optional
- Value of atrous_jmax PyBDSF parameter
- pad_to_size : int, optional
- Pad output mask image to a size of pad_to_size x pad_to_size
- skip_source_detection : bool, optional
- If True, source detection is not run on the input image
- region_file : str, optional
- Filename of region file in CASA format. If given, no mask image
- is made (the region file is used as the clean mask)
- nsig : float, optional
- Number of sigma of returned threshold value
- reference_ra_deg : float, optional
- RA for center of output mask image
- reference_dec_deg : float, optional
- Dec for center of output mask image
- cellsize_deg : float, optional
- Size of a pixel in degrees
- use_adaptive_threshold : bool, optional
- If True, use an adaptive threshold estimated from the negative values in
- the image
- adaptive_thresh : float, optional
- If adaptive_rmsbox is True, this value sets the threshold above
- which a source will use the small rms box
-
- Returns
- -------
- result : dict
- Dict with nsig-sigma rms threshold
-
- """
- if rmsbox is not None and type(rmsbox) is str:
- rmsbox = eval(rmsbox)
-
- if type(rmsbox_bright) is str:
- rmsbox_bright = eval(rmsbox_bright)
-
- if pad_to_size is not None and type(pad_to_size) is str:
- pad_to_size = int(pad_to_size)
-
- if type(atrous_do) is str:
- if atrous_do.lower() == 'true':
- atrous_do = True
- threshisl = 4.0 # override user setting to ensure proper source fitting
- else:
- atrous_do = False
-
- if type(iterate_threshold) is str:
- if iterate_threshold.lower() == 'true':
- iterate_threshold = True
- else:
- iterate_threshold = False
-
- if type(adaptive_rmsbox) is str:
- if adaptive_rmsbox.lower() == 'true':
- adaptive_rmsbox = True
- else:
- adaptive_rmsbox = False
-
- if type(skip_source_detection) is str:
- if skip_source_detection.lower() == 'true':
- skip_source_detection = True
- else:
- skip_source_detection = False
-
- if type(use_adaptive_threshold) is str:
- if use_adaptive_threshold.lower() == 'true':
- use_adaptive_threshold = True
- else:
- use_adaptive_threshold = False
-
- if reference_ra_deg is not None and reference_dec_deg is not None:
- reference_ra_deg = float(reference_ra_deg)
- reference_dec_deg = float(reference_dec_deg)
-
- if not os.path.exists(image_name):
- print('Input image not found. Making empty image...')
- if not skip_source_detection:
- print('ERROR: Source detection cannot be done on an empty image')
- sys.exit(1)
- if reference_ra_deg is not None and reference_dec_deg is not None:
- image_name = mask_name + '.tmp'
- make_template_image(image_name, reference_ra_deg, reference_dec_deg,
- cellsize_deg=float(cellsize_deg))
- else:
- print('ERROR: if image not found, a refernce position must be given')
- sys.exit(1)
-
- trim_by = float(trim_by)
- atrous_jmax = int(atrous_jmax)
- threshpix = float(threshpix)
- threshisl = float(threshisl)
- nsig = float(nsig)
- adaptive_thresh = float(adaptive_thresh)
- threshold = 0.0
-
- if not skip_source_detection:
- if vertices_file is not None:
- # Modify the input image to blank the regions outside of the polygon
- temp_img = pim.image(image_name)
- image_name += '.blanked'
- temp_img.saveas(image_name, overwrite=True)
- input_img = pim.image(image_name)
- data = input_img.getdata()
-
- vertices = read_vertices(vertices_file)
- RAverts = vertices[0]
- Decverts = vertices[1]
- xvert = []
- yvert = []
- for RAvert, Decvert in zip(RAverts, Decverts):
- pixels = input_img.topixel([1, 1, Decvert*np.pi/180.0,
- RAvert*np.pi/180.0])
- xvert.append(pixels[2]) # x -> Dec
- yvert.append(pixels[3]) # y -> RA
- poly = Polygon(xvert, yvert)
-
- # Find masked regions
- masked_ind = np.where(data[0, 0])
-
- # Find distance to nearest poly edge and set to NaN those that
- # are outside the facet (dist < 0)
- dist = poly.is_inside(masked_ind[0], masked_ind[1])
- outside_ind = np.where(dist < 0.0)
- if len(outside_ind[0]) > 0:
- data[0, 0, masked_ind[0][outside_ind], masked_ind[1][outside_ind]] = np.nan
-
- # Save changes
- input_img.putdata(data)
-
- if use_adaptive_threshold:
- # Get an estimate of the rms
- img = bdsf.process_image(image_name, mean_map='zero', rms_box=rmsbox,
- thresh_pix=threshpix, thresh_isl=threshisl,
- atrous_do=atrous_do, ini_method='curvature', thresh='hard',
- adaptive_rms_box=adaptive_rmsbox, adaptive_thresh=adaptive_thresh,
- rms_box_bright=rmsbox_bright, rms_map=True, quiet=True,
- atrous_jmax=atrous_jmax, stop_at='isl')
-
- # Find min and max pixels
- max_neg_val = abs(np.min(img.ch0_arr))
- max_neg_pos = np.where(img.ch0_arr == np.min(img.ch0_arr))
- max_pos_val = abs(np.max(img.ch0_arr))
- max_pos_pos = np.where(img.ch0_arr == np.max(img.ch0_arr))
-
- # Estimate new thresh_isl from min pixel value's sigma, but don't let
- # it get higher than 1/2 of the peak's sigma
- threshisl_neg = 2.0 * max_neg_val / img.rms_arr[max_neg_pos][0]
- max_sigma = max_pos_val / img.rms_arr[max_pos_pos][0]
- if threshisl_neg > max_sigma / 2.0:
- threshisl_neg = max_sigma / 2.0
-
- # Use the new threshold only if it is larger than the user-specified one
- if threshisl_neg > threshisl:
- threshisl = threshisl_neg
-
- if iterate_threshold:
- # Start with given threshold and lower it until we get at least one island
- nisl = 0
- while nisl == 0:
- img = bdsf.process_image(image_name, mean_map='zero', rms_box=rmsbox,
- thresh_pix=threshpix, thresh_isl=threshisl,
- atrous_do=atrous_do, ini_method='curvature', thresh='hard',
- adaptive_rms_box=adaptive_rmsbox, adaptive_thresh=adaptive_thresh,
- rms_box_bright=rmsbox_bright, rms_map=True, quiet=True,
- atrous_jmax=atrous_jmax)
- nisl = img.nisl
- threshpix /= 1.2
- threshisl /= 1.2
- if threshpix < 5.0:
- break
- else:
- img = bdsf.process_image(image_name, mean_map='zero', rms_box=rmsbox,
- thresh_pix=threshpix, thresh_isl=threshisl,
- atrous_do=atrous_do, ini_method='curvature', thresh='hard',
- adaptive_rms_box=adaptive_rmsbox, adaptive_thresh=adaptive_thresh,
- rms_box_bright=rmsbox_bright, rms_map=True, quiet=True,
- atrous_jmax=atrous_jmax)
-
- if img.nisl == 0:
- if region_file is None or region_file == '[]':
- print('No islands found. Clean mask cannot be made.')
- sys.exit(1)
- else:
- # Continue on and use user-supplied region file
- skip_source_detection = True
- threshold = nsig * img.clipped_rms
-
- # Check if there are large islands preset (indicating that multi-scale
- # clean is needed)
- has_large_isl = False
- for isl in img.islands:
- if isl.size_active > 100:
- # Assuming normal sampling, a size of 100 pixels would imply
- # a source of ~ 10 beams
- has_large_isl = True
-
- if (region_file is not None and region_file != '[]' and skip_source_detection):
- # Copy region file and return if source detection was not done
- os.system('cp {0} {1}'.format(region_file.strip('[]"'), mask_name))
- if threshold_format == 'float':
- return {'threshold_5sig': threshold}
- elif threshold_format == 'str_with_units':
- # This is done to get around the need for quotes around strings in casapy scripts
- # 'casastr/' is removed by the generic pipeline
- return {'threshold_5sig': 'casastr/{0}Jy'.format(threshold)}
- elif not skip_source_detection:
- img.export_image(img_type='island_mask', mask_dilation=0, outfile=mask_name,
- img_format=img_format, clobber=True)
-
- if (vertices_file is not None or trim_by > 0 or pad_to_size is not None
- or (region_file is not None and region_file != '[]')
- or skip_source_detection):
- # Alter the mask in various ways
- if skip_source_detection:
- # Read the image
- mask_im = pim.image(image_name)
- else:
- # Read the PyBDSF mask
- mask_im = pim.image(mask_name)
- data = mask_im.getdata()
- coordsys = mask_im.coordinates()
- if reference_ra_deg is not None and reference_dec_deg is not None:
- values = coordsys.get_referencevalue()
- values[2][0] = reference_dec_deg/180.0*np.pi
- values[2][1] = reference_ra_deg/180.0*np.pi
- coordsys.set_referencevalue(values)
- imshape = mask_im.shape()
- del(mask_im)
-
- if pad_to_size is not None:
- imsize = pad_to_size
- coordsys['direction'].set_referencepixel([imsize/2, imsize/2])
- pixmin = (imsize - imshape[2]) / 2
- if pixmin < 0:
- print("The padded size must be larger than the original size.")
- sys.exit(1)
- pixmax = pixmin + imshape[2]
- data_pad = np.zeros((1, 1, imsize, imsize), dtype=np.float32)
- data_pad[0, 0, pixmin:pixmax, pixmin:pixmax] = data[0, 0]
- new_mask = pim.image('', shape=(1, 1, imsize, imsize), coordsys=coordsys)
- new_mask.putdata(data_pad)
- else:
- new_mask = pim.image('', shape=imshape, coordsys=coordsys)
- new_mask.putdata(data)
-
- data = new_mask.getdata()
-
- if skip_source_detection:
- # Mask all pixels
- data[:] = 1
-
- if vertices_file is not None:
- # Modify the clean mask to exclude regions outside of the polygon
- vertices = read_vertices(vertices_file)
- RAverts = vertices[0]
- Decverts = vertices[1]
- xvert = []
- yvert = []
- for RAvert, Decvert in zip(RAverts, Decverts):
- try:
- pixels = new_mask.topixel([0, 1, Decvert*np.pi/180.0,
- RAvert*np.pi/180.0])
- except:
- pixels = new_mask.topixel([1, 1, Decvert*np.pi/180.0,
- RAvert*np.pi/180.0])
- xvert.append(pixels[2]) # x -> Dec
- yvert.append(pixels[3]) # y -> RA
- poly = Polygon(xvert, yvert)
-
- # Find masked regions
- masked_ind = np.where(data[0, 0])
-
- # Find distance to nearest poly edge and unmask those that
- # are outside the facet (dist < 0)
- dist = poly.is_inside(masked_ind[0], masked_ind[1])
- outside_ind = np.where(dist < 0.0)
- if len(outside_ind[0]) > 0:
- data[0, 0, masked_ind[0][outside_ind], masked_ind[1][outside_ind]] = 0
-
- if trim_by > 0.0:
- sh = np.shape(data)
- margin = int(sh[2] * trim_by / 2.0 )
- data[0, 0, 0:sh[2], 0:margin] = 0
- data[0, 0, 0:margin, 0:sh[3]] = 0
- data[0, 0, 0:sh[2], sh[3]-margin:sh[3]] = 0
- data[0, 0, sh[2]-margin:sh[2], 0:sh[3]] = 0
-
- if region_file is not None and region_file != '[]':
- # Merge the CASA regions with the mask
- casa_polys = read_casa_polys(region_file.strip('[]"'), new_mask)
- for poly in casa_polys:
- # Find unmasked regions
- unmasked_ind = np.where(data[0, 0] == 0)
-
- # Find distance to nearest poly edge and mask those that
- # are inside the casa region (dist > 0)
- dist = poly.is_inside(unmasked_ind[0], unmasked_ind[1])
- inside_ind = np.where(dist > 0.0)
- if len(inside_ind[0]) > 0:
- data[0, 0, unmasked_ind[0][inside_ind], unmasked_ind[1][inside_ind]] = 1
-
- # Save changes
- new_mask.putdata(data)
- if img_format == 'fits':
- new_mask.tofits(mask_name, overwrite=True)
- elif img_format == 'casa':
- new_mask.saveas(mask_name, overwrite=True)
- else:
- print('Output image format "{}" not understood.'.format(img_format))
- sys.exit(1)
-
- if not skip_source_detection:
- if threshold_format == 'float':
- return {'threshold_5sig': nsig * img.clipped_rms, 'multiscale': has_large_isl}
- elif threshold_format == 'str_with_units':
- # This is done to get around the need for quotes around strings in casapy scripts
- # 'casastr/' is removed by the generic pipeline
- return {'threshold_5sig': 'casastr/{0}Jy'.format(nsig * img.clipped_rms),
- 'multiscale': has_large_isl}
- else:
- return {'threshold_5sig': '0.0'}
-
-
-if __name__ == '__main__':
- descriptiontext = "Make a clean mask.\n"
-
- parser = argparse.ArgumentParser(description=descriptiontext, formatter_class=RawTextHelpFormatter)
- parser.add_argument('image_name', help='Image name')
- parser.add_argument('mask_name', help='Mask name')
- parser.add_argument('-a', '--atrous_do', help='use wavelet fitting', type=bool, default=False)
- parser.add_argument('-i', '--threshisl', help='', type=float, default=3.0)
- parser.add_argument('-p', '--threshpix', help='', type=float, default=5.0)
- parser.add_argument('-r', '--rmsbox', help='rms box width and step (e.g., "(60, 20)")',
- type=str, default='(60, 20)')
- parser.add_argument('--rmsbox_bright', help='rms box for bright sources(?) width and step (e.g., "(60, 20)")',
- type=str, default='(60, 20)')
- parser.add_argument('-t', '--iterate_threshold', help='iteratively decrease threshold until at least '
- 'one island is found', type=bool, default=False)
- parser.add_argument('-o', '--adaptive_rmsbox', help='use an adaptive rms box', type=bool, default=False)
- parser.add_argument('-f', '--img_format', help='format of output mask', type=str, default='casa')
- parser.add_argument('-d', '--threshold_format', help='format of return value', type=str, default='float')
- parser.add_argument('-b', '--trim_by', help='Trim masked region by this number of pixels', type=float, default=0.0)
- parser.add_argument('-v', '--vertices_file', help='file containing facet polygon vertices', type=str, default=None)
- parser.add_argument('--region_file', help='File containing casa regions to be merged with the detected mask', type=str, default=None)
- parser.add_argument('-j', '--atrous_jmax', help='Max wavelet scale', type=int, default=3)
- parser.add_argument('-z', '--pad_to_size', help='pad mask to this size', type=int, default=None)
- parser.add_argument('-s', '--skip_source_detection', help='skip source detection', type=bool, default=False)
-
- args = parser.parse_args()
- erg = main(args.image_name, args.mask_name, atrous_do=args.atrous_do,
- threshisl=args.threshisl, threshpix=args.threshpix, rmsbox=args.rmsbox,
- rmsbox_bright=args.rmsbox_bright,
- iterate_threshold=args.iterate_threshold,
- adaptive_rmsbox=args.adaptive_rmsbox, img_format=args.img_format,
- threshold_format=args.threshold_format, trim_by=args.trim_by,
- vertices_file=args.vertices_file, atrous_jmax=args.atrous_jmax,
- pad_to_size=args.pad_to_size, skip_source_detection=args.skip_source_detection,
- region_file=args.region_file)
- print erg
diff --git a/Docker/scripts/make_source_list.py b/Docker/scripts/make_source_list.py
deleted file mode 100755
index 3ce987507c2497694dad1ac79098e780e0452043..0000000000000000000000000000000000000000
--- a/Docker/scripts/make_source_list.py
+++ /dev/null
@@ -1,136 +0,0 @@
-#! /usr/bin/env python
-"""
-Script to make a source catalog for an image
-"""
-import argparse
-from argparse import RawTextHelpFormatter
-import sys
-import os
-try:
- import matplotlib
- matplotlib.use('Agg')
-except (RuntimeError, ImportError):
- pass
-import bdsf
-import lsmtool
-
-
-def main(image_name, catalog_name, atrous_do=False, threshisl=3.0, threshpix=5.0, rmsbox=(60, 20),
- rmsbox_bright=(35, 7), adaptive_rmsbox=False, atrous_jmax=6, adaptive_thresh=150.0,
- compare_dir=None, format='png'):
- """
- Make a source catalog for an image
-
- Parameters
- ----------
- image_name : str
- Filename of input image from which catalog will be made
- catalog_name : str
- Filename of output source catalog in FITS format
- atrous_do : bool, optional
- Use wavelet module of PyBDSF?
- threshisl : float, optional
- Value of thresh_isl PyBDSF parameter
- threshpix : float, optional
- Value of thresh_pix PyBDSF parameter
- rmsbox : tuple of floats, optional
- Value of rms_box PyBDSF parameter
- rmsbox_bright : tuple of floats, optional
- Value of rms_box_bright PyBDSF parameter
- adaptive_rmsbox : tuple of floats, optional
- Value of adaptive_rms_box PyBDSF parameter
- atrous_jmax : int, optional
- Value of atrous_jmax PyBDSF parameter
- adaptive_thresh : float, optional
- If adaptive_rmsbox is True, this value sets the threshold above
- which a source will use the small rms box
- compare_dir : str, optional
- Directory for output plots/info of comparison of source properties to surveys
- format : str, optional
- Format of output plots: 'png' or 'pdf'
- """
- if type(rmsbox) is str:
- rmsbox = eval(rmsbox)
-
- if type(rmsbox_bright) is str:
- rmsbox_bright = eval(rmsbox_bright)
-
- if type(atrous_do) is str:
- if atrous_do.lower() == 'true':
- atrous_do = True
- threshisl = 4.0 # override user setting to ensure proper source fitting
- else:
- atrous_do = False
-
- if type(adaptive_rmsbox) is str:
- if adaptive_rmsbox.lower() == 'true':
- adaptive_rmsbox = True
- else:
- adaptive_rmsbox = False
-
- if not os.path.exists(image_name):
- print('ERROR: input image not found')
- sys.exit(1)
-
- atrous_jmax = int(atrous_jmax)
- threshpix = float(threshpix)
- threshisl = float(threshisl)
- adaptive_thresh = float(adaptive_thresh)
-
- img = bdsf.process_image(image_name, mean_map='zero', rms_box=rmsbox,
- thresh_pix=threshpix, thresh_isl=threshisl,
- atrous_do=atrous_do, ini_method='curvature', thresh='hard',
- adaptive_rms_box=adaptive_rmsbox, adaptive_thresh=adaptive_thresh,
- rms_box_bright=rmsbox_bright, rms_map=True, quiet=True,
- atrous_jmax=atrous_jmax)
- srcroot = 'prefactor'
- img.write_catalog(outfile=catalog_name, bbs_patches='source', srcroot=srcroot, clobber=True)
-
- # Use LSMTool to make some basic comparisons to surveys
- if compare_dir is not None:
- s = lsmtool.load(catalog_name)
- s.setPatchPositions(method='wmean')
- _, _, refRA, refDec = s._getXY()
- def_dict = s.getDefaultValues()
- if 'ReferenceFrequency' in def_dict:
- ref_freq_hz = def_dict['ReferenceFrequency']
- else:
- ref_freq_hz = 0.0
- if ref_freq_hz > 140e6 and ref_freq_hz < 160e6:
- # use TGSS
- vocat = 'TGSS'
- ignspec = None
- else:
- # too far in freq from TGSS -> use GSM, ignoring sources that lack spectral
- # information (spec_indx = -0.7)
- vocat = 'GSM'
- ignspec = -0.7
- s_vo = lsmtool.load(vocat, VOPosition=[refRA, refDec], VORadius='5 deg')
- s_vo.group('threshold', FWHM='30 arcsec', threshold=0.001)
- s.compare(s_vo, radius='30 arcsec', excludeMultiple=False, outDir=compare_dir,
- name1='LOFAR', name2=vocat, format=format, ignoreSpec=ignspec)
-
-
-if __name__ == '__main__':
- descriptiontext = "Make a source catalog from an image.\n"
-
- parser = argparse.ArgumentParser(description=descriptiontext, formatter_class=RawTextHelpFormatter)
- parser.add_argument('image_name', help='Image name')
- parser.add_argument('catalog_name', help='Output catalog name')
- parser.add_argument('-a', '--atrous_do', help='use wavelet fitting', type=bool, default=False)
- parser.add_argument('-i', '--threshisl', help='', type=float, default=3.0)
- parser.add_argument('-p', '--threshpix', help='', type=float, default=5.0)
- parser.add_argument('-r', '--rmsbox', help='rms box width and step (e.g., "(60, 20)")',
- type=str, default='(60, 20)')
- parser.add_argument('--rmsbox_bright', help='rms box for bright sources(?) width and step (e.g., "(60, 20)")',
- type=str, default='(60, 20)')
- parser.add_argument('-o', '--adaptive_rmsbox', help='use an adaptive rms box', type=bool, default=False)
- parser.add_argument('-j', '--atrous_jmax', help='Max wavelet scale', type=int, default=3)
- parser.add_argument('-c', '--compare_dir', help='', type=str, default=None)
- parser.add_argument('-f', '--format', help='', type=str, default='pdf')
-
- args = parser.parse_args()
- main(args.image_name, args.catalog_name, atrous_do=args.atrous_do,
- threshisl=args.threshisl, threshpix=args.threshpix, rmsbox=args.rmsbox,
- rmsbox_bright=args.rmsbox_bright, adaptive_rmsbox=args.adaptive_rmsbox,
- atrous_jmax=args.atrous_jmax, compare_dir=args.compare_dir, format=args.format)
diff --git a/Docker/scripts/make_summary.py b/Docker/scripts/make_summary.py
deleted file mode 100755
index 70c230b217a9ff27281c643e666c5e60729671f0..0000000000000000000000000000000000000000
--- a/Docker/scripts/make_summary.py
+++ /dev/null
@@ -1,226 +0,0 @@
-#! /usr/bin/env python
-# -*- coding: utf-8 -*-
-
-import os,sys
-import argparse
-import subprocess, platform
-import numpy, re
-import multiprocessing
-
-########################################################################
-def input2strlist_nomapfile(invar):
- """
- from bin/download_IONEX.py
- give the list of MSs from the list provided as a string
- """
- str_list = None
- if type(invar) is str:
- if invar.startswith('[') and invar.endswith(']'):
- str_list = [f.strip(' \'\"') for f in invar.strip('[]').split(',')]
- else:
- str_list = [invar.strip(' \'\"')]
- elif type(invar) is list:
- str_list = [str(f).strip(' \'\"') for f in invar]
- else:
- raise TypeError('input2strlist: Type '+str(type(invar))+' unknown!')
- return str_list
-
-###############################################################################
-def find_flagged_fraction(ms_file):
-
- outputs = os.popen('DPPP msin=' + ms_file + ' msout=. steps=[count] count.type=counter count.warnperc=0.0000001 | grep NOTE').readlines()
- fraction_flagged = { output.split('(')[-1].rstrip(')\n'):output.split('%')[0][-5:].strip() for output in outputs if 'station' in output }
- return fraction_flagged
-
-###############################################################################
-def main(observation_directory = '/data/share/pipeline/Observation', logfile = '/data/share/pipeline/Observation/logs/pipeline.log', h5parmdb = 'solutions.h5', MSfile = '[]'):
- """
- Creates summary of a given prefactor3 run
-
- Parameters
- ----------
- observation_directory: directory where all the processed data and solutions are stored
- logfile: directory where all the log files are stored
- h5parmdb: name of the solutions h5parm database
- MSfile: pattern of the final MS files
-
- """
- ID = os.path.basename(observation_directory)
- summary = os.path.normpath(observation_directory + '/../' + ID + '.log')
- f_summary = open(summary, 'w')
-
- # location of logfile
- print('Summary logfile is written to ' + summary)
-
- if not os.path.exists(logfile):
- f_summary.write(logfile + ' does not exist. Skipping!')
- f_summary.close()
- return(1)
- pass
-
- ## check for the h5parm
- if os.path.exists(observation_directory + '/results/cal_values/' + h5parmdb):
- h5parmdb = observation_directory + '/results/cal_values/' + h5parmdb
- pass
- elif os.path.exists(observation_directory + '/results/cal_values/cal_' + h5parmdb):
- h5parmdb = observation_directory + '/results/cal_values/cal_' + h5parmdb
- else:
- f_summary.write('No h5parm solutions file found in ' + observation_directory + '/results/cal_values/' + h5parmdb)
- f_summary.close()
- return(1)
- pass
-
- ## convert stringlist to list
- mslist = input2strlist_nomapfile(MSfile)
-
- ## get proper solset
- from losoto.h5parm import h5parm
- data = h5parm(h5parmdb, readonly = True)
- if 'target' in data.getSolsetNames():
- solset = 'target'
- pass
- elif 'calibrator' in data.getSolsetNames():
- solset = 'calibrator'
- pass
- else:
- f_summary.write('Neither calibrator or target solset has been found in ' + h5parmdb)
- f_summary.close()
- return(1)
- pass
-
- ## get antenna information
- solset = data.getSolset(solset)
- soltabs = list([soltab.name for soltab in solset.getSoltabs()])
- source = solset.obj._f_get_child('source')[0][0]
- antenna_len = '{:<' + str(max([ len(antenna) for antenna in solset.getAnt()])) + '}'
- soltab_len = '{:^' + str(max([ len(soltab_name) for soltab_name in soltabs ])) + '}'
-
- ## get software versions
- try:
- DPPP_version = subprocess.Popen(['DPPP', '-v'], stdout=subprocess.PIPE).stdout.read()
- losoto_version = 'losoto ' + subprocess.Popen(['losoto', '--version'], stderr=subprocess.PIPE).stderr.read()
- lsmtool_version = subprocess.Popen(['lsmtool', '--version'], stdout=subprocess.PIPE).stdout.read()
- aoflagger_version = subprocess.Popen(['aoflagger', '--version'], stdout=subprocess.PIPE).stdout.read()
- wsclean_version = subprocess.Popen(['wsclean', '--version'], stdout=subprocess.PIPE).stdout.read().split('\n')[1] + '\n'
- python_version = subprocess.Popen(['python', '--version'], stderr=subprocess.PIPE).stderr.read()
-
- import matplotlib, scipy, astropy
- modules_version = 'matplotlib ' + str(matplotlib.__version__) + ', scipy ' + str(scipy.__version__) + ', astropy ' + str(astropy.__version__) + '\n'
- os_version = ' '.join(platform.linux_distribution()) + ' ' + platform.release() + ' ' + platform.version() + '\n'
- f_summary.write('Software versions currently used:\n' + os_version \
- + DPPP_version \
- + aoflagger_version \
- + losoto_version \
- + lsmtool_version \
- + wsclean_version \
- + python_version \
- + modules_version )
- except OSError:
- f_summary.write('Could not find all LOFAR or related software packages. Could not determine the used software versions.\n')
- pass
-
- ## get the list of removed stations:
- baseline_list = list(set([ str(line.split(";")[1:]) for line in open(logfile) if 'baseline:' in line and ';' in line ]))
- bad_antennas = list(set(re.sub("[\[\]\"\ '!*]", "", str(baseline_list)).replace('\\n','').replace('\\','').split(',')))
- if len(bad_antennas) == 1 and bad_antennas[0] == '':
- bad_antennas = 'NONE'
- pass
-
- ## check for A-Team warning:
- Ateam_list = list(set([ str(line.split('source ')[-1]).split(' is')[0] for line in open(logfile) if 'WARNING: The A-Team source' in line ]))
- if len(Ateam_list) == 0:
- Ateam_list = 'NONE'
- pass
-
- f_summary.write('\n')
- f_summary.write('Antennas removed from the data: ' + re.sub("[\[\]']", "", str(bad_antennas)) + '\n')
- f_summary.write('A-Team sources close to the phase reference center: ' + re.sub("[\[\]']", "", str(Ateam_list)) + '\n')
- f_summary.write('\n')
-
- ## diffractive scale
- structure_file = observation_directory + '/results/inspection/' + source + '_structure.txt'
- if os.path.exists(structure_file):
- with open(structure_file, 'r') as infile:
- for line in infile:
- diffractive_scale_XX = float(line.split()[2])
- diffractive_scale_YY = float(line.split()[2])
- f_summary.write('XX diffractive scale: %3.1f km'%(diffractive_scale_XX/1000.) + '\n')
- f_summary.write('YY diffractive scale: %3.1f km'%(diffractive_scale_YY/1000.) + '\n')
-
- ## check whether reference has been used
- if 'bandpass' in soltabs:
- soltab = solset.getSoltab('bandpass')
- history = soltab.getHistory()
- if not history.strip('\n') == '':
- f_summary.write(history + '\n')
- pass
-
- ## get table of flagged solutions
- f_summary.write('\n')
- flagged_solutions = {}
- import losoto.lib_operations as losoto
- for soltab_name in soltabs:
- soltab = solset.getSoltab(soltab_name)
- antennas = soltab.ant
- axes = soltab.getAxesNames()
- axes.insert(0, axes.pop(axes.index('ant'))) ## put antenna table to the second place
- flagged_solutions[soltab_name] = {}
- for vals, weights, coord, selection in soltab.getValuesIter(returnAxes=axes, weight=True):
- vals = losoto.reorderAxes(vals, soltab.getAxesNames(), axes)
- weights = losoto.reorderAxes(weights, soltab.getAxesNames(), axes)
- for i, antenna in enumerate(antennas):
- flagged_fraction = 1. - float(numpy.mean(weights[i]))
- flagged_solutions[soltab_name][antenna] = soltab_len.format('{:.2f}'.format(100 * flagged_fraction) + '%')
- soltab_names = ' '.join([ soltab_len.format(soltab_name) for soltab_name in soltabs ])
- f_summary.write('Amount of flagged solutions per station and solution table:\n')
- f_summary.write(antenna_len.format('Station') + ' ' + soltab_names + '\n')
- for antenna in antennas:
- values_to_print = []
- for soltab_name in soltabs:
- try:
- values_to_print.append(flagged_solutions[soltab_name][antenna])
- except KeyError:
- values_to_print.append(soltab_len.format(' '))
- values_to_print = ' '.join(values_to_print)
- f_summary.write(antenna_len.format(antenna) + ' ' + values_to_print + '\n')
-
-
- ## derive the fraction of flagged data of the entire observation
- pool = multiprocessing.Pool(processes = multiprocessing.cpu_count())
- flagged_fraction_list = pool.map(find_flagged_fraction, mslist)
-
- flagged_fraction_data = {}
- for entry in flagged_fraction_list:
- antennas = entry.keys()
- for antenna in antennas:
- try:
- flagged_fraction_data[antenna].append(float(entry[antenna]))
- except KeyError:
- flagged_fraction_data[antenna] = [float(entry[antenna])]
-
- f_summary.write('\n')
- f_summary.write('Overall amount of flagged data in the final data:\n')
- f_summary.write(antenna_len.format('Station') + '\n')
- for antenna in sorted(flagged_fraction_data.keys()):
- flagged_fraction_data[antenna] = sum(flagged_fraction_data[antenna]) / len(flagged_fraction_data[antenna])
- f_summary.write(antenna_len.format(antenna) + ' ' + '{:>10}'.format('{:.2f}'.format(flagged_fraction_data[antenna]) + '%') + '\n')
-
- print('Summary has been created.')
- return 0
-
-
-if __name__=='__main__':
-
- parser = argparse.ArgumentParser(description='Creates summary of a given prefactor3 run.')
- parser.add_argument('--obsdir', type=str, default='/data/scratch/pipeline/Observation', help='Directory where to find the processed data and solutions.')
- parser.add_argument('--logfile', type=str, default='/data/share/pipeline/Observation/logs/pipeline.log', help='Location of the pipeline logfile.')
- parser.add_argument('--h5parm', '--h5parm', type=str, default='solutions.h5', help='Name of the h5parm solutions file (default: solutions.h5)')
- parser.add_argument('--MSfile', '--MSfile', type=str, default='[]', help='List of MS to be analysed (default: [])')
-
- args = parser.parse_args()
-
- # start running script
- main(args.obsdir, args.logfile, args.h5parm, args.MSfile)
-
- sys.exit(0)
- pass
diff --git a/Docker/scripts/merge_skymodels.py b/Docker/scripts/merge_skymodels.py
deleted file mode 100755
index 20ce30f8b35d3f5c2bedf298b0b4b384d65595ea..0000000000000000000000000000000000000000
--- a/Docker/scripts/merge_skymodels.py
+++ /dev/null
@@ -1,54 +0,0 @@
-#!/usr/bin/env python
-"""
-Script to merge two sky models
-"""
-import argparse
-from argparse import RawTextHelpFormatter
-import os
-import lsmtool
-import sys
-
-
-def main(inmodel1, inmodel2, outmodel, match_by='name', radius=0.0, keep='all'):
- """
- Creates mosaic
-
- Parameters
- ----------
- inmodel1 : str
- Filename of input model 1
- inmodel2 : str
- Filename of input model 2
- outmodel : str
- Filename of output model
- match_by : str, optional
- The match_by parameter of LSMTool
- radius : float, optional
- Radius in degrees for cross match
- keep : str, optional
- The keep parameter of LSMTool
-
- """
- s1 = lsmtool.load(inmodel1)
- s2 = lsmtool.load(inmodel2)
-
- s1.concatenate(s2, matchBy=match_by, radius=float(radius), keep=keep,
- inheritPatches=True)
- s1.group('every')
- s1.write(fileName=outmodel, clobber=True)
-
-
-if __name__ == '__main__':
- descriptiontext = "Merge two sky models.\n"
-
- parser = argparse.ArgumentParser(description=descriptiontext, formatter_class=RawTextHelpFormatter)
- parser.add_argument('inmodel1', help='input model 1 filename')
- parser.add_argument('inmodel2', help='input model 2 filename')
- parser.add_argument('outmodel', help='output model filename')
- parser.add_argument('-m', '--match_by', help='match type', type=str, default='name')
- parser.add_argument('-r', '--radius', help='radius in degrees for matching', type=float, default=0.0)
- parser.add_argument('-k', '--keep', help='', type=str, default='all')
-
- args = parser.parse_args()
- main(args.inmodel1, args.inmodel2, args.outmodel, match_by=args.match_by,
- radius=args.radius, keep=args.keep)
diff --git a/Docker/scripts/pad_image.py b/Docker/scripts/pad_image.py
deleted file mode 100755
index 15f81dce17ba53392bb35e66e1c2427a2d2d78ef..0000000000000000000000000000000000000000
--- a/Docker/scripts/pad_image.py
+++ /dev/null
@@ -1,44 +0,0 @@
-#! /usr/bin/env python
-"""
-Script to pad a FITS image
-"""
-import argparse
-from argparse import RawTextHelpFormatter
-from astropy.io import fits as pyfits
-import numpy as np
-import sys
-import os
-
-
-def main(infile, xypadsize):
- """Pad a fits image with zeros to padsize"""
- pad_xsize, pad_ysize = [int(s) for s in xypadsize.split(' ')]
-
- hdu = pyfits.open(infile)
- imdata = hdu[0].data[0, 0]
- (ysize, xsize) = imdata.shape
- if xsize > pad_xsize or ysize > pad_ysize:
- raise ValueError('pad_image: padded size is smaller than current size!')
- if xsize == pad_xsize and ysize == pad_ysize:
- return
-
- xoffset = (pad_xsize - xsize) / 2
- yoffset = (pad_ysize - ysize) / 2
- newdata=np.zeros((1, 1, pad_ysize, pad_xsize))
-
- newdata[0, 0, yoffset:yoffset+ysize, xoffset:xoffset+xsize] = imdata
- hdu[0].data = newdata
- hdu[0].header['CRPIX1'] += xoffset
- hdu[0].header['CRPIX2'] += yoffset
- hdu.writeto(infile, clobber=True)
-
-
-if __name__ == '__main__':
- descriptiontext = "Pad FITS image.\n"
-
- parser = argparse.ArgumentParser(description=descriptiontext, formatter_class=RawTextHelpFormatter)
- parser.add_argument('infile', help='filename of input image')
- parser.add_argument('xypadsize', help='padded size as "xsize ysize"')
- args = parser.parse_args()
-
- main(args.infile, args.xypadsize)
diff --git a/Docker/scripts/plot_Ateamclipper.py b/Docker/scripts/plot_Ateamclipper.py
deleted file mode 100755
index 56e5d169a2da5e7846dd0e777cd4765bcb15049c..0000000000000000000000000000000000000000
--- a/Docker/scripts/plot_Ateamclipper.py
+++ /dev/null
@@ -1,49 +0,0 @@
-#!/usr/bin/env python
-# -* coding: utf-8 -*-
-
-"""
-Adds phases (=0) and amplitudes (=1) to any missing station if they appear in an h5parm, but not in a particular soltab.
-
-Created on Tue Jul 24 2019
-
-@author: Alexander Drabent
-"""
-
-import argparse
-import matplotlib as mpl
-mpl.use('Agg')
-import matplotlib
-import matplotlib.pyplot as plt
-import numpy
-
-def main(txtfile = 'Ateamclipper.txt', outfile = 'Ateamclipper.png'):
-
- frac_list_xx = []
- frac_list_yy = []
- freq_list = []
- with open(txtfile, 'r') as infile:
- for line in infile:
- freq_list.append(float(line.split()[0]))
- frac_list_xx.append(float(line.split()[1]))
- frac_list_yy.append(float(line.split()[2]))
-
- # Plot the amount of clipped data vs. frequency potentially contaminated by the A-team
- plt.scatter(numpy.array(freq_list) / 1e6, numpy.array(frac_list_xx), marker = '.', s = 10)
- plt.xlabel('frequency [MHz]')
- plt.ylabel('A-team clipping fraction [%]')
- plt.savefig(outfile)
- return(0)
-
-if __name__ == "__main__":
- parser = argparse.ArgumentParser(description='Adds phases and amplitudes to any missing station if they appear in an h5parm, but not in a particular soltab.')
-
- parser.add_argument('txtfile', type=str,
- help='Input text file containing frequency and flag fraction of the XX and YY polarization.')
- parser.add_argument('outfile', type=str,
- help='Input text file containing frequency and flag fraction of the XX and YY polarization.')
-
-
- args = parser.parse_args()
-
- main(txtfile = args.txtfile, outfile = args.outfile)
-
diff --git a/Docker/scripts/plot_image.py b/Docker/scripts/plot_image.py
deleted file mode 100755
index 460865968769d15b33a3de15c420d47bc65682e2..0000000000000000000000000000000000000000
--- a/Docker/scripts/plot_image.py
+++ /dev/null
@@ -1,213 +0,0 @@
-#!/usr/bin/env python
-"""
-plot image
-Original code from Reinout
-"""
-import os
-import matplotlib as mpl
-mpl.use('Agg')
-import aplpy
-import numpy
-import astropy.io.fits
-
-
-def meanclip(indata, clipsig=4.0, maxiter=10, converge_num=0.001, verbose=0):
- """
- Computes an iteratively sigma-clipped mean on a
- data set. Clipping is done about median, but mean
- is returned.
-
- .. note:: MYMEANCLIP routine from ACS library.
-
- :History:
- * 21/10/1998 Written by RSH, RITSS
- * 20/01/1999 Added SUBS, fixed misplaced paren on float call, improved doc. RSH
- * 24/11/2009 Converted to Python. PLL.
-
- Examples
- --------
- >>> mean, sigma = meanclip(indata)
-
- Parameters
- ----------
- indata: array_like
- Input data.
-
- clipsig: float
- Number of sigma at which to clip.
-
- maxiter: int
- Ceiling on number of clipping iterations.
-
- converge_num: float
- If the proportion of rejected pixels is less than
- this fraction, the iterations stop.
-
- verbose: {0, 1}
- Print messages to screen?
-
- Returns
- -------
- mean: float
- N-sigma clipped mean.
-
- sigma: float
- Standard deviation of remaining pixels.
-
- """
- # Flatten array
- skpix = indata.reshape(indata.size,)
-
- ct = indata.size
- iter = 0
- c1 = 1.0
- c2 = 0.0
-
- while (c1 >= c2) and (iter < maxiter):
- lastct = ct
- medval = numpy.median(skpix)
- sig = numpy.std(skpix)
- wsm = numpy.where(abs(skpix-medval) < clipsig*sig)
- ct = len(wsm[0])
- if ct > 0:
- skpix = skpix[wsm]
-
- c1 = abs(ct - lastct)
- c2 = converge_num * lastct
- iter += 1
-
- mean = numpy.mean(skpix)
- sigma = robust_sigma(skpix)
-
- return mean, sigma
-
-
-def find_imagenoise(data):
- mean, rms = meanclip(data)
- return rms
-
-
-def robust_sigma(in_y, zero=0):
- """
- Calculate a resistant estimate of the dispersion of
- a distribution. For an uncontaminated distribution,
- this is identical to the standard deviation.
-
- Use the median absolute deviation as the initial
- estimate, then weight points using Tukey Biweight.
- See, for example, Understanding Robust and
- Exploratory Data Analysis, by Hoaglin, Mosteller
- and Tukey, John Wiley and Sons, 1983.
-
- .. note:: ROBUST_SIGMA routine from IDL ASTROLIB.
-
- Examples
- --------
- >>> result = robust_sigma(in_y, zero=1)
-
- Parameters
- ----------
- in_y : array_like
- Vector of quantity for which the dispersion is
- to be calculated
-
- zero : int
- If set, the dispersion is calculated w.r.t. 0.0
- rather than the central value of the vector. If
- Y is a vector of residuals, this should be set.
-
- Returns
- -------
- out_val : float
- Dispersion value. If failed, returns -1.
-
- """
- # Flatten array
- y = in_y.ravel()
-
- eps = 1.0E-20
- c1 = 0.6745
- c2 = 0.80
- c3 = 6.0
- c4 = 5.0
- c_err = -1.0
- min_points = 3
-
- if zero:
- y0 = 0.0
- else:
- y0 = numpy.median(y)
-
- dy = y - y0
- del_y = abs(dy)
-
- # First, the median absolute deviation MAD about the median:
- mad = numpy.median(del_y) / c1
-
- # If the MAD=0, try the MEAN absolute deviation:
- if mad < eps:
- mad = del_y.mean() / c2
- if mad < eps:
- return 0.0
-
- # Now the biweighted value:
- u = dy / (c3 * mad)
- uu = u * u
- q = numpy.where(uu <= 1.0)
- count = len(q[0])
- if count < min_points:
- print('ROBUST_SIGMA: This distribution is TOO WEIRD! '
- 'Returning {}'.format(c_err))
- return c_err
-
- numerator = numpy.sum((y[q] - y0)**2.0 * (1.0 - uu[q])**4.0)
- n = y.size
- den1 = numpy.sum((1.0 - uu[q]) * (1.0 - c4 * uu[q]))
- siggma = n * numerator / (den1 * (den1 - 1.0))
-
- if siggma > 0:
- out_val = numpy.sqrt(siggma)
- else:
- out_val = 0.0
-
- return out_val
-
-
-def main(fitsimage, outfilename, outdir=None):
-
- # find image noise
- hdulist = astropy.io.fits.open(fitsimage, mode='update')
- data = hdulist[0].data
- imagenoise = find_imagenoise(data)
- dynamicrange = numpy.max(data) / imagenoise
-
- # Store the rms value in the image header, for later reference by metadata.py (see
- # CEP/Pipeline/recipes/sip/helpers/metadata.py in the LOFAR software repository)
- header = hdulist[0].header
- header['MEANRMS'] = imagenoise
- beammaj = header['BMAJ']
- beammin = header['BMIN']
- beampa = header['BPA']
- hdulist.close()
-
- # plot
- f = aplpy.FITSFigure(fitsimage, slices=[0, 0])
- f.show_colorscale(vmax=16*imagenoise, vmin=-4*imagenoise, cmap='bone')
- f.add_grid()
- f.grid.set_color('white')
- f.grid.set_alpha(0.5)
- f.grid.set_linewidth(0.2)
- f.add_colorbar()
- f.colorbar.set_axis_label_text('Flux Density (Jy/beam)')
- f.add_beam()
- f.beam.set_frame(True)
- f.set_title('Mean rms = {0:1.2e} Jy/beam; Dynamic range = {1:1.2e}; \n '
- 'Restoring beam = ({2:3.1f} x {3:3.1f}) arcsec, PA = {4:3.1f} deg'.format(
- imagenoise, dynamicrange, beammaj*3600.0, beammin*3600.0, beampa))
-
- if outdir is not None:
- if not os.path.exists(outdir):
- os.makedirs(outdir)
- outfilename = os.path.join(outdir, os.path.basename(outfilename))
- f.save('{}.png'.format(outfilename), dpi=400)
-
diff --git a/Docker/scripts/plot_subtract_images.py b/Docker/scripts/plot_subtract_images.py
deleted file mode 100755
index 3ffd2bad58f2abbc67d7aec4f7a8f323bde9a947..0000000000000000000000000000000000000000
--- a/Docker/scripts/plot_subtract_images.py
+++ /dev/null
@@ -1,251 +0,0 @@
-#!/usr/bin/env python
-"""
-plot image overlaid with mask contours
-Original code from Reinout
-"""
-import argparse
-from argparse import RawTextHelpFormatter
-#import pyrap.images as pim
-#import pyrap.tables as pt
-#import numpy as np
-import matplotlib as mpl
-mpl.use('Agg')
-import aplpy
-import numpy
-import astropy.io.fits
-
-def meanclip(indata, clipsig=4.0, maxiter=10, converge_num=0.001, verbose=0):
- """
- Computes an iteratively sigma-clipped mean on a
- data set. Clipping is done about median, but mean
- is returned.
-
- .. note:: MYMEANCLIP routine from ACS library.
-
- :History:
- * 21/10/1998 Written by RSH, RITSS
- * 20/01/1999 Added SUBS, fixed misplaced paren on float call, improved doc. RSH
- * 24/11/2009 Converted to Python. PLL.
-
- Examples
- --------
- >>> mean, sigma = meanclip(indata)
-
- Parameters
- ----------
- indata: array_like
- Input data.
-
- clipsig: float
- Number of sigma at which to clip.
-
- maxiter: int
- Ceiling on number of clipping iterations.
-
- converge_num: float
- If the proportion of rejected pixels is less than
- this fraction, the iterations stop.
-
- verbose: {0, 1}
- Print messages to screen?
-
- Returns
- -------
- mean: float
- N-sigma clipped mean.
-
- sigma: float
- Standard deviation of remaining pixels.
-
- """
- # Flatten array
- skpix = indata.reshape( indata.size, )
-
- ct = indata.size
- iter = 0; c1 = 1.0 ; c2 = 0.0
-
- while (c1 >= c2) and (iter < maxiter):
- lastct = ct
- medval = numpy.median(skpix)
- sig = numpy.std(skpix)
- wsm = numpy.where( abs(skpix-medval) < clipsig*sig )
- ct = len(wsm[0])
- if ct > 0:
- skpix = skpix[wsm]
-
- c1 = abs(ct - lastct)
- c2 = converge_num * lastct
- iter += 1
- # End of while loop
-
- mean = numpy.mean( skpix )
- sigma = robust_sigma( skpix )
-
- if verbose:
- prf = 'MEANCLIP:'
- print '%s %.1f-sigma clipped mean' % (prf, clipsig)
- print '%s Mean computed in %i iterations' % (prf, iter)
- print '%s Mean = %.6f, sigma = %.6f' % (prf, mean, sigma)
-
- return mean, sigma
-
-
-def find_imagenoise(data):
- mean, rms = meanclip(data)
- return rms
-
-def robust_sigma(in_y, zero=0):
- """
- Calculate a resistant estimate of the dispersion of
- a distribution. For an uncontaminated distribution,
- this is identical to the standard deviation.
-
- Use the median absolute deviation as the initial
- estimate, then weight points using Tukey Biweight.
- See, for example, Understanding Robust and
- Exploratory Data Analysis, by Hoaglin, Mosteller
- and Tukey, John Wiley and Sons, 1983.
-
- .. note:: ROBUST_SIGMA routine from IDL ASTROLIB.
-
- Examples
- --------
- >>> result = robust_sigma(in_y, zero=1)
-
- Parameters
- ----------
- in_y : array_like
- Vector of quantity for which the dispersion is
- to be calculated
-
- zero : int
- If set, the dispersion is calculated w.r.t. 0.0
- rather than the central value of the vector. If
- Y is a vector of residuals, this should be set.
-
- Returns
- -------
- out_val : float
- Dispersion value. If failed, returns -1.
-
- """
- # Flatten array
- y = in_y.ravel()
-
- eps = 1.0E-20
- c1 = 0.6745
- c2 = 0.80
- c3 = 6.0
- c4 = 5.0
- c_err = -1.0
- min_points = 3
-
- if zero:
- y0 = 0.0
- else:
- y0 = numpy.median(y)
-
- dy = y - y0
- del_y = abs( dy )
-
- # First, the median absolute deviation MAD about the median:
-
- mad = numpy.median( del_y ) / c1
-
- # If the MAD=0, try the MEAN absolute deviation:
- if mad < eps:
- mad = del_y.mean() / c2
- if mad < eps:
- return 0.0
-
- # Now the biweighted value:
- u = dy / (c3 * mad)
- uu = u * u
- q = numpy.where(uu <= 1.0)
- count = len(q[0])
- if count < min_points:
- module_logger.warn('ROBUST_SIGMA: This distribution is TOO WEIRD! '
- 'Returning {}'.format(c_err))
- return c_err
-
- numerator = numpy.sum( (y[q] - y0)**2.0 * (1.0 - uu[q])**4.0 )
- n = y.size
- den1 = numpy.sum( (1.0 - uu[q]) * (1.0 - c4 * uu[q]) )
- siggma = n * numerator / ( den1 * (den1 - 1.0) )
-
- if siggma > 0:
- out_val = numpy.sqrt( siggma )
- else:
- out_val = 0.0
-
- return out_val
-
-
-
-
-def main(fitsimage, maskfits, outfilename):
-
- #imagenoise = 5E-3
-
- #f = aplpy.FITSFigure(fitsimage,slices=[0,0])
- #f.show_colorscale(vmax=10*imagenoise,vmin=-3*imagenoise)
- ##outplotname = fitsimage.replace('.fits','.png')
- #f.add_colorbar()
- #f.save('%s.png'%(outfilename),dpi=750)
- #f.show_contour(maskfits, levels=[0.5], colors='m')
- #f.save('%s_mask.png'%(outfilename),dpi=750)
-
- ##f = aplpy.FITSFigure(maskfits,slices=[0,0])
- ##f.show_colorscale(vmax=0,vmin=1)
- ###f.show_contour(maskfits, levels=[0.5], colors='w')
- ###outplotname = fitsimage.replace('.fits','.png')
- ##f.add_colorbar()
- ##f.save('%s-mask.png'%(outfilename),dpi=500)
-
-
-
- #outplotname = fitsimage.replace('.fits','.png')
-
- # find image noise
- hdulist = astropy.io.fits.open(fitsimage)
- data = hdulist[0].data
- imagenoise = find_imagenoise(data)
- hdulist.close()
- print 'Image noise is: ', imagenoise*1e3, ' mJy'
-
-
- f = aplpy.FITSFigure(fitsimage,slices=[0,0])
- f.show_colorscale(vmax=16*imagenoise,vmin=-4*imagenoise, cmap='bone')
-
-
- #f.add_beam()
- #f.beam.set_color('white')
- #f.beam.set_hatch('.')
-
-
- f.add_grid()
- f.grid.set_color('white')
- f.grid.set_alpha(0.5)
- f.grid.set_linewidth(0.2)
-
- f.add_colorbar()
-
-
- f.show_contour(maskfits, colors='red', levels=[0.0],filled=False, smooth=1,alpha=0.6, linewidths=0.25)
- f.colorbar.set_axis_label_text('Flux (Jy beam$^{-1}$)')
-
- f.save('%s.png'%(outfilename),dpi=400)
-
-
-
-
-if __name__ == '__main__':
- descriptiontext = "Convert a fits image to a CASA image.\n"
-
- parser = argparse.ArgumentParser(description=descriptiontext, formatter_class=RawTextHelpFormatter)
- parser.add_argument('fitsimage', help='name of the image')
- parser.add_argument('maskfits', help='name of the mask')
- parser.add_argument('outname', help='name for the output')
- args = parser.parse_args()
-
- main(args.fitsimage, args.maskfits, args.outname)
diff --git a/Docker/scripts/plot_unflagged_fraction.py b/Docker/scripts/plot_unflagged_fraction.py
deleted file mode 100755
index 49485b265db30a0e84eae39a0e663eaa869322f5..0000000000000000000000000000000000000000
--- a/Docker/scripts/plot_unflagged_fraction.py
+++ /dev/null
@@ -1,57 +0,0 @@
-#!/usr/bin/env python
-import argparse
-import glob
-import signal
-import os
-import sys
-import matplotlib as mpl
-mpl.use('Agg')
-import matplotlib
-import matplotlib.pyplot as plt
-import numpy as np
-import casacore.tables as pt
-from matplotlib import cm
-
-
-def input2strlist_nomapfile(invar):
- """
- from bin/download_IONEX.py
- give the list of MSs from the list provided as a string
- """
- str_list = None
- if type(invar) is str:
- if invar.startswith('[') and invar.endswith(']'):
- str_list = [f.strip(' \'\"') for f in invar.strip('[]').split(',')]
- else:
- str_list = [invar.strip(' \'\"')]
- elif type(invar) is list:
- str_list = [str(f).strip(' \'\"') for f in invar]
- else:
- raise TypeError('input2strlist: Type '+str(type(invar))+' unknown!')
- return str_list
-
-
-def main(ms_list, frac_list, outfile='unflagged_fraction.png'):
-
- ms_list = input2strlist_nomapfile(ms_list)
- frac_list = input2strlist_nomapfile(frac_list)
- frac_list = np.array([float(f) for f in frac_list])
- outdir = os.path.dirname(outfile)
- if not os.path.exists(outdir):
- os.makedirs(outdir)
-
- # Get frequencies
- freq_list = []
- for ms in ms_list:
- # open the main table and print some info about the MS
- t = pt.table(ms, readonly=True, ack=False)
- tfreq = pt.table(t.getkeyword('SPECTRAL_WINDOW'),readonly=True,ack=False)
- ref_freq = tfreq.getcol('REF_FREQUENCY',nrow=1)[0]
- freq_list.append(ref_freq)
- freq_list = np.array(freq_list) / 1e6 # MHz
-
- # Plot the unflagged fraction vs. frequency
- plt.scatter(freq_list, frac_list)
- plt.xlabel('frequency [MHz]')
- plt.ylabel('unflagged fraction')
- plt.savefig(outfile)
diff --git a/Docker/scripts/plot_uvcov.py b/Docker/scripts/plot_uvcov.py
deleted file mode 100755
index f752725b12096d78df85cf08e93cb8457f4d46db..0000000000000000000000000000000000000000
--- a/Docker/scripts/plot_uvcov.py
+++ /dev/null
@@ -1,224 +0,0 @@
-#!/usr/bin/env python
-import argparse
-import glob
-import signal
-import os
-import sys
-import matplotlib as mpl
-mpl.use('Agg')
-import matplotlib
-import matplotlib.pyplot as plt
-import numpy
-import casacore.tables as pt
-from matplotlib import cm
-
-
-def input2strlist_nomapfile(invar):
- """
- from bin/download_IONEX.py
- give the list of MSs from the list provided as a string
- """
- str_list = None
- if type(invar) is str:
- if invar.startswith('[') and invar.endswith(']'):
- str_list = [f.strip(' \'\"') for f in invar.strip('[]').split(',')]
- else:
- str_list = [invar.strip(' \'\"')]
- elif type(invar) is list:
- str_list = [str(f).strip(' \'\"') for f in invar]
- else:
- raise TypeError('input2strlist: Type '+str(type(invar))+' unknown!')
- return str_list
-
-
-def string2bool(instring):
- if isinstance(instring, bool):
- return instring
- if instring.upper() == 'TRUE' or instring == '1':
- return True
- elif instring.upper() == 'FALSE' or instring == '0':
- return False
- else:
- raise ValueError('string2bool: Cannot convert string "'+instring+'" to boolean!')
-
-
-def main(input, output, title='uv coverage', limits=',,,', timeslots='0,10,0', antennas='-1',
- wideband=True):
-
- MSlist = input2strlist_nomapfile(input)
- if len(MSlist) == 0:
- print('Error: You must specify at least one MS name.')
- sys.exit(1)
- plottitle = title
- fileformat = output.split('.')[-1]
- if fileformat not in ['png','pdf','eps','ps']:
- print('Error: Unknown file extension')
- sys.exit(1)
- axlimits = limits.strip().split(',')
- if len(axlimits) == 4:
- xmin,xmax,ymin,ymax = axlimits
- else:
- print('Error: You must specify four axis limits')
- sys.exit(1)
- timeslots = timeslots.split(',')
- if len(timeslots) != 3:
- print('Error: Timeslots format is start,skip,end')
- sys.exit(1)
- for i in range(len(timeslots)):
- timeslots[i] = int(timeslots[i])
- if timeslots[i] < 0:
- print('Error: timeslots values must not be negative')
- sys.exit(1)
- antToPlotSpl = antennas.split(',')
- antToPlot = []
- for i in range(len(antToPlotSpl)):
- tmpspl = antToPlotSpl[i].split('..')
- if len(tmpspl) == 1:
- antToPlot.append(int(antToPlotSpl[i]))
- elif len(tmpspl) == 2:
- for j in range(int(tmpspl[0]),int(tmpspl[1])+1):
- antToPlot.append(j)
- else:
- print('Error: Could not understand antenna list.')
- sys.exit(1)
- wideband = string2bool(wideband)
-
- outdir = os.path.dirname(output)
- if not os.path.exists(outdir):
- os.makedirs(outdir)
-
- badval = 0.0
- xaxisvals = []
- yaxisvals = []
- flagvals = []
- savex = numpy.array([])
- savey = numpy.array([])
- prev_firstTime = None
- prev_lastTime = None
- prev_intTime = None
- for inputMS in MSlist:
- # open the main table and print some info about the MS
- t = pt.table(inputMS, readonly=True, ack=False)
- tfreq = pt.table(t.getkeyword('SPECTRAL_WINDOW'),readonly=True,ack=False)
- ref_freq = tfreq.getcol('REF_FREQUENCY',nrow=1)[0]
- ch_freq = tfreq.getcol('CHAN_FREQ',nrow=1)[0]
- if wideband:
- ref_wavelength = 2.99792458e8/ch_freq
- else:
- ref_wavelength = [2.99792458e8/ref_freq]
-
- firstTime = t.getcell("TIME", 0)
- lastTime = t.getcell("TIME", t.nrows()-1)
- intTime = t.getcell("INTERVAL", 0)
- nTimeslots = (lastTime - firstTime) / intTime
- if timeslots[1] == 0:
- timeskip = 1
- else:
- timeskip = int(timeslots[1])
- if timeslots[2] == 0:
- timeslots[2] = nTimeslots
- # open the antenna subtable
- tant = pt.table(t.getkeyword('ANTENNA'), readonly=True, ack=False)
-
- # Station names
- antList = tant.getcol('NAME')
- if len(antToPlot)==1 and antToPlot[0]==-1:
- antToPlot = list(range(len(antList)))
-
- # select by time from the beginning, and only use specified antennas
- tsel = t.query('TIME >= %f AND TIME <= %f AND ANTENNA1 IN %s AND ANTENNA2 IN %s' % (firstTime+timeslots[0]*intTime,firstTime+timeslots[2]*intTime,str(antToPlot),str(antToPlot)), columns='ANTENNA1,ANTENNA2,UVW,FLAG_ROW')
-
- # Now we loop through the baselines
- i = 0
- nb = (len(antToPlot)*(len(antToPlot)-1))/2
- for tpart in tsel.iter(["ANTENNA1","ANTENNA2"]):
- ant1 = tpart.getcell("ANTENNA1", 0)
- ant2 = tpart.getcell("ANTENNA2", 0)
- if ant1 not in antToPlot or ant2 not in antToPlot: continue
- if ant1 == ant2: continue
- i += 1
- # Get the values to plot (only need to read again if this ms has new times)
- if prev_firstTime != firstTime or prev_lastTime != lastTime or prev_intTime != intTime:
- uvw = tpart.getcol('UVW', rowincr=timeskip)
- savex = numpy.append(savex,[uvw[:,0],-uvw[:,0]])
- savey = numpy.append(savey,[uvw[:,1],-uvw[:,1]])
-
- # Get the flags
- flags = tpart.getcol('FLAG_ROW', rowincr=timeskip)
- for w in ref_wavelength:
- flagvals.extend(flags.tolist())
- flagvals.extend(flags.tolist())
- for w in ref_wavelength:
- xaxisvals.extend((savex/w/1000.).tolist())
- yaxisvals.extend((savey/w/1000.).tolist())
- prev_firstTime = firstTime
- prev_lastTime = lastTime
- prev_intTime = intTime
-
- # Plot the uv coverage
- xaxisvals = numpy.array(xaxisvals)
- yaxisvals = numpy.array(yaxisvals)
- zvals = numpy.zeros(xaxisvals.shape)
- flagvals = numpy.array(flagvals)
- flagged_ind = numpy.where(flagvals)
- zvals[flagged_ind] = 1.0
- uvmax = max(xaxisvals.max(),yaxisvals.max())
- uvmin = min(xaxisvals.min(),yaxisvals.min())
- uvuplim = 0.02*(uvmax-uvmin)+uvmax
- uvlolim = uvmin-0.02*(uvmax-uvmin)
- if xmin == '':
- minx = uvlolim
- else:
- minx = float(xmin)
- if xmax == '':
- maxx = uvuplim
- else:
- maxx = float(xmax)
- if ymin == '':
- miny = uvlolim
- else:
- miny = float(ymin)
- if ymax == '':
- maxy = uvuplim
- else:
- maxy = float(ymax)
- if minx == maxx:
- minx = -1.0
- maxx = 1.0
- if miny == maxy:
- miny = -1.0
- maxy = 1.0
- plt.xlabel(r'u [k$\lambda$]')
- plt.ylabel(r'v [k$\lambda$]')
- plt.xlim([minx,maxx])
- plt.ylim([miny,maxy])
- gridsize = 300
- plt.hexbin(xaxisvals, yaxisvals, C=zvals, gridsize=gridsize, cmap=cm.jet, bins=None)
- plt.title(plottitle)
- plt.axes().set_aspect('equal')
- plt.grid(True)
- cb = plt.colorbar()
- cb.set_label('flagged fraction')
- plt.savefig(output)
-
- # Plot histogram of the flagged fraction vs. uv distance
- uvdist = numpy.sqrt(xaxisvals**2 + yaxisvals**2)
- sorted_ind = numpy.argsort(uvdist)
- uvdist = uvdist[sorted_ind]
- zvals = 1.0-zvals[sorted_ind]
- hist, bins = numpy.histogram(uvdist, 50)
- bins_mean = [0.5 * (bins[i] + bins[i+1]) for i in range(len(bins)-1)]
- inds = []
- for b in bins[:-1]:
- inds.append(numpy.where(uvdist > b)[0][0])
- inds.append(len(bins))
- flagged_mean = [numpy.mean(zvals[inds[i]:inds[i+1]]) for i in range(len(inds)-1)]
- plt.clf()
- plt.bar(bins_mean, flagged_mean, (bins[1]-bins[0])/2.0)
- plt.xlabel(r'uv distance [k$\lambda$]')
- plt.ylabel('unflagged fraction')
- plt.savefig('{0}_uvdist{1}'.format(*os.path.splitext(output)))
-
-
-
-
diff --git a/Docker/scripts/sort_times_into_freqGroups.py b/Docker/scripts/sort_times_into_freqGroups.py
deleted file mode 100755
index b97e2b0ce533cff3f55e0a67bc2814ad7f522bd8..0000000000000000000000000000000000000000
--- a/Docker/scripts/sort_times_into_freqGroups.py
+++ /dev/null
@@ -1,422 +0,0 @@
-#! /usr/bin/env python
-"""
-Script to sort a list of MSs by into frequency groups by time-stamp
-"""
-import pyrap.tables as pt
-import sys, os
-import numpy as np
-from lofarpipe.support.data_map import DataMap
-from lofarpipe.support.data_map import DataProduct
-
-def _calc_edge_chans(inmap, numch, edgeFactor=32):
- """
- Generates a map with strings that can be used as input for NDPPP to flag the edges
- of the input MSs during (or after) concatenation.
-
- inmap - MultiDataMap (not mapfilename!) with the files to be concatenated.
- numch - Number of channels per input file (All files are assumed to have the same number
- of channels.)
- edgeFactor - Divisor to compute how many channels are to be flagged at beginning and end.
- (numch=64 and edgeFactor=32 means "flag two channels at beginning and two at end")
- """
- outmap = DataMap([])
- for group in inmap:
- flaglist = []
- for i in xrange(len(group.file)):
- flaglist.extend(range(i*numch,i*numch+numch/edgeFactor))
- flaglist.extend(range((i+1)*numch-numch/edgeFactor,(i+1)*numch))
- outmap.append(DataProduct(group.host,str(flaglist).replace(' ',''),group.skip))
- print '_calc_edge_chans: flaglist:', str(flaglist).replace(' ','')
- return outmap
-
-def input2bool(invar):
- if invar == None:
- return None
- if isinstance(invar, bool):
- return invar
- elif isinstance(invar, str):
- if invar.upper() == 'TRUE' or invar == '1':
- return True
- elif invar.upper() == 'FALSE' or invar == '0':
- return False
- else:
- raise ValueError('input2bool: Cannot convert string "'+invar+'" to boolean!')
- elif isinstance(invar, int) or isinstance(invar, float):
- return bool(invar)
- else:
- raise TypeError('input2bool: Unsupported data type:'+str(type(invar)))
-
-def input2int(invar):
- if invar == None:
- return None
- if isinstance(invar, int):
- return invar
- elif isinstance(invar, float):
- return int(invar)
- elif isinstance(invar, str):
- if invar.strip().upper() == 'NONE' or invar == 'FALSE':
- return None
- else:
- return int(invar)
- else:
- raise TypeError('input2int: Unsupported data type:'+str(type(invar)))
-
-
-
-def main(ms_input, filename=None, mapfile_dir=None, numSB=-1, hosts=None, NDPPPfill=True, target_path=None, stepname=None,
- mergeLastGroup=False, truncateLastSBs=True, firstSB=None):
- """
- Check a list of MS files for missing frequencies
-
- Parameters
- ----------
- ms_input : list or str
- List of MS filenames, or string with list, or path to a mapfile
- filename: str
- Name of output mapfile
- mapfile_dir : str
- Directory for output mapfile
- numSB : int, optional
- How many files should go into one frequency group. Values <= 0 mean put
- all files of the same time-step into one group.
- default = -1
- hosts : list or str
- List of hostnames or string with list of hostnames
- NDPPPfill : bool, optional
- Add dummy file-names for missing frequencies, so that NDPPP can
- fill the data with flagged dummy data.
- default = True
- target_path : str, optional
- Change the path of the "groups" files to this. (I.e. write output files
- into this directory with the subsequent NDPPP call.)
- default = keep path of input files
- stepname : str, optional
- Add this step-name into the file-names of the output files.
- mergeLastGroup, truncateLastSBs : bool, optional
- mergeLastGroup = True, truncateLastSBs = True:
- not allowed
- mergeLastGroup = True, truncateLastSBs = False:
- put the files from the last group that doesn't have SBperGroup subbands
- into the second last group (which will then have more than SBperGroup entries).
- mergeLastGroup = False, truncateLastSBs = True:
- ignore last files, that don't make for a full group (not all files are used).
- mergeLastGroup = False, truncateLastSBs = False:
- keep inclomplete last group, or - with NDPPPfill=True - fill
- last group with dummies.
- firstSB : int, optional
- If set, then reference the grouping of files to this station-subband. As if a file
- with this station-subband would be included in the input files.
- (For HBA-low, i.e. 0 -> 100MHz, 55 -> 110.74MHz, 512 -> 200MHz)
-
- Returns
- -------
- result : dict
- Dict with the name of the generated mapfile
-
- """
-
- NDPPPfill = input2bool(NDPPPfill)
- mergeLastGroup = input2bool(mergeLastGroup)
- truncateLastSBs = input2bool(truncateLastSBs)
- firstSB = input2int(firstSB)
- numSB = int(numSB)
-
- if not filename or not mapfile_dir:
- raise ValueError('sort_times_into_freqGroups: filename and mapfile_dir are needed!')
- if mergeLastGroup and truncateLastSBs:
- raise ValueError('sort_times_into_freqGroups: Can either merge the last partial group or truncate at last full group, not both!')
-# if mergeLastGroup:
-# raise ValueError('sort_times_into_freqGroups: mergeLastGroup is not (yet) implemented!')
- if type(ms_input) is str:
- if ms_input.startswith('[') and ms_input.endswith(']'):
- ms_list = [f.strip(' \'\"') for f in ms_input.strip('[]').split(',')]
- else:
- map_in = DataMap.load(ms_input)
- map_in.iterator = DataMap.SkipIterator
- ms_list = []
- for fname in map_in:
- if fname.startswith('[') and fname.endswith(']'):
- for f in fname.strip('[]').split(','):
- ms_list.append(f.strip(' \'\"'))
- else:
- ms_list.append(fname.strip(' \'\"'))
- elif type(ms_input) is list:
- ms_list = [str(f).strip(' \'\"') for f in ms_input]
- else:
- raise TypeError('sort_times_into_freqGroups: type of "ms_input" unknown!')
-
- if type(hosts) is str:
- hosts = [h.strip(' \'\"') for h in hosts.strip('[]').split(',')]
- if not hosts:
- hosts = ['localhost']
- numhosts = len(hosts)
- print "sort_times_into_freqGroups: Working on",len(ms_list),"files (including flagged files)."
-
- time_groups = {}
- # sort by time
- for i, ms in enumerate(ms_list):
- # work only on files selected by a previous step
- if ms.lower() != 'none':
- # use the slower but more reliable way:
- obstable = pt.table(ms, ack=False)
- timestamp = int(round(np.min(obstable.getcol('TIME'))))
- #obstable = pt.table(ms+'::OBSERVATION', ack=False)
- #timestamp = int(round(obstable.col('TIME_RANGE')[0][0]))
- obstable.close()
- if timestamp in time_groups:
- time_groups[timestamp]['files'].append(ms)
- else:
- time_groups[timestamp] = {'files': [ ms ], 'basename' : os.path.splitext(ms)[0] }
- print "sort_times_into_freqGroups: found",len(time_groups),"time-groups"
-
- # sort time-groups by frequency
- timestamps = time_groups.keys()
- timestamps.sort() # not needed now, but later
- first = True
- nchans = 0
- for time in timestamps:
- freqs = []
- for ms in time_groups[time]['files']:
- # Get the frequency info
- sw = pt.table(ms+'::SPECTRAL_WINDOW', ack=False)
- freq = sw.col('REF_FREQUENCY')[0]
- if first:
- file_bandwidth = sw.col('TOTAL_BANDWIDTH')[0]
- nchans = sw.col('CHAN_WIDTH')[0].shape[0]
- chwidth = sw.col('CHAN_WIDTH')[0][0]
- freqset = set([freq])
- first = False
- else:
- assert file_bandwidth == sw.col('TOTAL_BANDWIDTH')[0]
- assert nchans == sw.col('CHAN_WIDTH')[0].shape[0]
- assert chwidth == sw.col('CHAN_WIDTH')[0][0]
- freqset.add(freq)
- freqs.append(freq)
- sw.close()
- time_groups[time]['freq_names'] = zip(freqs,time_groups[time]['files'])
- time_groups[time]['freq_names'].sort(key=lambda pair: pair[0])
- #time_groups[time]['files'] = [name for (freq,name) in freq_names]
- #time_groups[time]['freqs'] = [freq for (freq,name) in freq_names]
- print "sort_times_into_freqGroups: Collected the frequencies for the time-groups"
-
- freqliste = np.array(list(freqset))
- freqliste.sort()
- freq_width = np.min(freqliste[1:]-freqliste[:-1])
- if file_bandwidth > freq_width:
- raise ValueError("Bandwidth of files is larger than minimum frequency step between two files!")
- if file_bandwidth < (freq_width*0.51):
- raise ValueError("Bandwidth of files is smaller than 51% of the minimum frequency step between two files! (More than about half the data is missing.)")
- #the new output map
- filemap = MultiDataMap()
- groupmap = DataMap()
- # add 1% of the SB badwidth in case maxfreq might be "exactly" on a group-border
- maxfreq = np.max(freqliste)+freq_width*0.51
- if firstSB != None:
- if freqliste[0] < 100e6:
- # LBA Data
- minfreq = (float(firstSB)/512.*100e6)-freq_width/2.
- elif freqliste[0] > 100e6 and freqliste[0] < 200e6:
- # HBA-Low
- minfreq = (float(firstSB)/512.*100e6)+100e6-freq_width/2.
- elif freqliste[0] > 200e6 and freqliste[0] < 300e6:
- # HBA-high
- minfreq = (float(firstSB)/512.*100e6)+200e6-freq_width/2.
- else:
- raise ValueError('sort_times_into_freqGroups: Frequency of lowest input data is higher than 300MHz!')
- if np.min(freqliste) < minfreq:
- raise ValueError('sort_times_into_freqGroups: Frequency of lowest input data is lower than reference frequency!')
- else:
- minfreq = np.min(freqliste)-freq_width/2.
- groupBW = freq_width*numSB
- if groupBW < 1e6 and groupBW > 0:
- print 'sort_times_into_freqGroups: ***WARNING***: Bandwidth of concatenated MS is lower than 1 MHz. This may cause conflicts with the concatenated file names!'
- if groupBW < 0:
- # this is the case for concatenating all subbands
- groupBW = maxfreq-minfreq
- truncateLastSBs = input2bool(False)
- NDPPPfill = input2bool(True)
- freqborders = np.arange(minfreq,maxfreq,groupBW)
- if mergeLastGroup:
- freqborders[-1] = maxfreq
- elif truncateLastSBs:
- pass #nothing to do! # left to make the logic more clear!
- elif not truncateLastSBs and NDPPPfill:
- freqborders = np.append(freqborders,(freqborders[-1]+groupBW))
- elif not truncateLastSBs and not NDPPPfill:
- freqborders = np.append(freqborders,maxfreq)
-
- freqborders = freqborders[freqborders>(np.min(freqliste)-groupBW)]
- ngroups = len(freqborders)-1
- if ngroups == 0:
- raise ValueError('sort_times_into_freqGroups: Not enough input subbands to create at least one full (frequency-)group!')
-
- print "sort_times_into_freqGroups: Will create",ngroups,"group(s) with",numSB,"file(s) each."
-
- hostID = 0
- for time in timestamps:
- (freq,fname) = time_groups[time]['freq_names'].pop(0)
- for groupIdx in xrange(ngroups):
- files = []
- skip_this = True
- filefreqs_low = np.arange(freqborders[groupIdx],freqborders[groupIdx+1],freq_width)
- for lower_freq in filefreqs_low:
- if freq > lower_freq and freq < lower_freq+freq_width:
- assert freq!=1e12
- files.append(fname)
- if len(time_groups[time]['freq_names'])>0:
- (freq,fname) = time_groups[time]['freq_names'].pop(0)
- else:
- (freq,fname) = (1e12,'This_shouldn\'t_show_up')
- skip_this = False
- elif NDPPPfill:
- files.append('dummy.ms')
- if not skip_this:
- filemap.append(MultiDataProduct(hosts[hostID%numhosts], files, skip_this))
- freqID = int((freqborders[groupIdx]+freqborders[groupIdx+1])/2e6)
- groupname = time_groups[time]['basename']+'_%Xt_%dMHz.ms'%(time,freqID)
- if type(stepname) is str:
- groupname += stepname
- if type(target_path) is str:
- groupname = os.path.join(target_path,os.path.basename(groupname))
- groupmap.append(DataProduct(hosts[hostID%numhosts],groupname, skip_this))
- orphan_files = len(time_groups[time]['freq_names'])
- if freq < 1e12:
- orphan_files += 1
- if orphan_files > 0:
- print "sort_times_into_freqGroups: Had %d unassigned files in time-group %xt."%(orphan_files, time)
- filemapname = os.path.join(mapfile_dir, filename)
- filemap.save(filemapname)
- groupmapname = os.path.join(mapfile_dir, filename+'_groups')
- groupmap.save(groupmapname)
- # genertate map with edge-channels to flag
- flagmap = _calc_edge_chans(filemap, nchans)
- flagmapname = os.path.join(mapfile_dir, filename+'_flags')
- flagmap.save(flagmapname)
- result = {'mapfile': filemapname, 'groupmapfile': groupmapname, 'flagmapfile': flagmapname}
- return result
-
-class MultiDataProduct(DataProduct):
- """
- Class representing multiple files in a DataProduct.
- """
- def __init__(self, host=None, file=None, skip=True):
- super(MultiDataProduct, self).__init__(host, file, skip)
- if not file:
- self.file = list()
- else:
- self._set_file(file)
-
- def __repr__(self):
- """Represent an instance as a Python dict"""
- return (
- "{{'host': '{0}', 'file': {1}, 'skip': {2}}}".format(self.host, self.file, str(self.skip))
- )
-
- def __str__(self):
- """Print an instance as 'host:[filelist]'"""
- return ':'.join((self.host, str(self.file)))
-
- def _set_file(self, data):
- try:
- # Try parsing as a list
- if isinstance(data, list):
- self.file = data
- if isinstance(data, DataProduct):
- self._from_dataproduct(data)
- if isinstance(data, DataMap):
- self._from_datamap(data)
-
- except TypeError:
- raise DataProduct("No known method to set a filelist from %s" % str(file))
-
- def _from_dataproduct(self, prod):
- print 'setting filelist from DataProduct'
- self.host = prod.host
- self.file = prod.file
- self.skip = prod.skip
-
- def _from_datamap(self, inmap):
- print 'setting filelist from DataMap'
- filelist = {}
- for item in inmap:
- if not item.host in filelist:
- filelist[item.host] = []
- filelist[item.host].append(item.file)
- self.file = filelist['i am']
-
- def append(self, item):
- self.file.append(item)
-
-
-class MultiDataMap(DataMap):
- """
- Class representing a specialization of data-map, a collection of data
- products located on the same node, skippable as a set and individually
- """
- @DataMap.data.setter
- def data(self, data):
- if isinstance(data, DataMap):
- mdpdict = {}
- data.iterator = DataMap.SkipIterator
- for item in data:
- if not item.host in mdpdict:
- mdpdict[item.host] = []
- mdpdict[item.host].append(item.file)
- mdplist = []
- for k, v in mdpdict.iteritems():
- mdplist.append(MultiDataProduct(k, v, False))
- self._set_data(mdplist, dtype=MultiDataProduct)
- elif isinstance(data, MultiDataProduct):
- self._set_data(data, dtype=MultiDataProduct)
- elif not data:
- pass
- else:
- self._set_data(data, dtype=MultiDataProduct)
-
- def split_list(self, number):
- mdplist = []
- for item in self.data:
- for i in xrange(0, len(item.file), number):
- chunk = item.file[i:i+number]
- mdplist.append(MultiDataProduct(item.host, chunk, item.skip))
- self._set_data(mdplist)
-
-
-
-if __name__ == '__main__':
- import optparse
- import glob
- import random
-
- opt = optparse.OptionParser(usage='%prog [options] <MSPattern> \n')
- opt.add_option('-v', '--verbose', help='Go Vebose! (default=False)', action='store_true', default=False)
- opt.add_option('-r', '--randomize', help='Randomize order of the input files. (default=False)', action='store_true', default=False)
- opt.add_option('-d', '--decimate', help='Remove every 10th file (after randomization if that is done). (default=False)', action='store_true', default=False)
- opt.add_option('-n', '--numbands', help='Number of how many files should be grouped together in frequency. (default=all files in one group)', type='int', default=-1)
- opt.add_option('-f', '--filename', help='Name for the mapfiles to write. (default=\"test.mapfile\")', type='string', default='test.mapfile')
-
- (options, args) = opt.parse_args()
-
- # Check options
- if len(args) != 1:
- opt.print_help()
- sys.exit()
-
- # first argument: pattern for measurement-sets
- inMSs = glob.glob(args[0])
- if options.randomize:
- random.shuffle(inMSs)
- if options.decimate:
- for i in range((len(inMSs)-1),-1,-10):
- inMSs.pop(i)
-
- ergdict = main(inMSs, options.filename, '.', numSB=options.numbands, hosts=None, NDPPPfill=True)
-
- groupmap = DataMap.load(ergdict['groupmapfile'])
- filemap = MultiDataMap.load(ergdict['mapfile'])
- print "len(groupmap) : %d , len(filemap) : %d " % (len(groupmap),len(filemap))
- if len(groupmap) != len(filemap):
- print "groupmap and filemap have different length!"
- sys.exit(1)
- for i in xrange(len(groupmap)):
- print "Group \"%s\" has %d entries."%(groupmap[i].file,len(filemap[i].file))
diff --git a/Docker/scripts/transfer_solutions.py b/Docker/scripts/transfer_solutions.py
deleted file mode 100755
index 6b59e3e7e236a12ed7421bf96d270d72d9932943..0000000000000000000000000000000000000000
--- a/Docker/scripts/transfer_solutions.py
+++ /dev/null
@@ -1,178 +0,0 @@
-#!/usr/bin/env python
-# Written for prefactor by Alexander Drabent (alex@tls-tautenburg.de), 19th July 2019.
-
-from losoto.h5parm import h5parm
-from losoto.lib_operations import *
-
-import os
-import re
-import numpy
-import subprocess
-########################################################################
-def main(h5parmdb, refh5parm, insolset='sol000', outsolset='sol000', insoltab='amplitude000', outsoltab='amplitude000', antenna = '[FUSPID].*', trusted_sources = ['3C48', '3C147'], parset = None):
-
-
- ### Open up the h5parm, get an example value
- data = h5parm(h5parmdb, readonly=False)
- refdata = h5parm(refh5parm, readonly=True)
- insolset = refdata.getSolset(insolset)
- outsolset = data.getSolset(outsolset)
- insoltab = insolset.getSoltab(insoltab)
- outsoltab = outsolset.getSoltab(outsoltab, useCache = True)
-
- ### Get the calibrator information from the target solset and check whether it is trusted or not.
- calibrators = []
- sources = outsolset.obj._f_get_child('source')
- for i in numpy.arange(len(sources)):
- calibrators.append(sources[i][0])
- calibrator = numpy.unique(calibrators)
- if len(calibrator) > 1:
- logging.error('There is more than one calibrator used in the target solution set: ' + str(calibrator) + '. No solutions will be transferred.')
- data.close()
- refdata.close()
- return(1)
- if calibrator[0].upper() in trusted_sources:
- logging.info(calibrator[0] + ' is a trusted calibrator source. No solutions from a reference solution set will be transferred!')
- data.close()
- refdata.close()
- return(0)
- else:
- logging.info(calibrator[0] + ' is not a trusted calibrator source. Solutions from a reference solution set will be transferred!')
-
-
- ### check for matching antennas
- station_names = outsoltab.ant
- stations_to_transfer = [ station_name for station_name in station_names if re.match(antenna, station_name) ]
- if len(stations_to_transfer) == 0:
- logging.warning('No stations found matching the regular expression: ' + antenna)
- data.close()
- refdata.close()
- return(0)
-
- ### Properly define axes order
- out_axes = outsoltab.getAxesNames()
- out_axes.insert(0, out_axes.pop(out_axes.index('time'))) ## put time table to the first place
- out_axes.insert(1, out_axes.pop(out_axes.index('freq'))) ## put frequency table to the second place
- out_axes.insert(2, out_axes.pop(out_axes.index('ant'))) ## put antenna table to the third place
- in_axes = insoltab.getAxesNames()
- in_axes.insert(0, in_axes.pop(in_axes.index('freq'))) ## put frequency table to the second place
- in_axes.insert(1, in_axes.pop(in_axes.index('ant'))) ## put antenna table to the third place
-
- if 'time' in in_axes:
- logging.warning('Reference solution table provides a time axis. The bandpass is assumed to be constant in time. Thus, only the first time step will be transferred.')
- in_axes.insert(2, in_axes.pop(in_axes.index('time'))) ## put time table to the third place
-
- ### resort the val and weights array according to the new axes order
- for outvals, outweights, coord, selection in outsoltab.getValuesIter(returnAxes=out_axes, weight=True):
- outvals = reorderAxes( outvals, outsoltab.getAxesNames(), out_axes)
- outweights = reorderAxes( outweights, outsoltab.getAxesNames(), out_axes )
- pass
-
- ### resort the val and weights array according to the new axes order
- for invals, inweights, coord, selection in insoltab.getValuesIter(returnAxes=in_axes, weight=True):
- invals = reorderAxes( invals, insoltab.getAxesNames(), in_axes)
- inweights = reorderAxes( inweights, insoltab.getAxesNames(), in_axes )
- pass
-
- ### check for equistance in frequency
- freq_resolution = numpy.unique(numpy.diff(outsoltab.freq))
- if len(freq_resolution) > 1:
- logging.error('Frequency axis is not equidistant!')
- data.close()
- refdata.close()
- return(1)
-
- ### look for nearest neighbours in frequency and write them into a dictionary
- nearest_freq = {}
- freq_idx = []
- for i, frequency in enumerate(numpy.array(outsoltab.freq)):
- nearest_freq[i] = numpy.abs(frequency - numpy.array(insoltab.freq)).argmin()
- freq_idx.append(nearest_freq[i])
-
- ### do the job!
- for station_to_transfer in stations_to_transfer:
- logging.info('Transferring solutions for antenna: ' + str(station_to_transfer))
- unique_elements, counts_elements = numpy.unique(freq_idx, return_counts = True)
- ant_index = list(outsoltab.ant).index(station_to_transfer)
- try:
- ref_ant_index = list(insoltab.ant).index(station_to_transfer)
- except ValueError:
- logging.warning('Station ' + station_to_transfer + ' does not exist in reference soltab. No solutions will be transferred!')
- continue
- for j, frequency in enumerate(numpy.array(outsoltab.freq)):
- if 'time' in in_axes:
- weight = inweights[nearest_freq[j], ref_ant_index, 0]
- else:
- weight = inweights[nearest_freq[j], ref_ant_index]
- if counts_elements[nearest_freq[j]] > 1:
- if frequency < insoltab.freq[0] or frequency > insoltab.freq[-1]:
- counts_elements[nearest_freq[j]] -= 1
- logging.warning('No entry for frequency ' + str(frequency) + '. Solution will be flagged.')
- weight = 0
- for i, time in enumerate(outsoltab.time):
- outweights[i,j,ant_index] = weight
- if 'time' in in_axes:
- outvals[i,j,ant_index] = invals[nearest_freq[j], ref_ant_index, 0]
- else:
- outvals[i,j,ant_index] = invals[nearest_freq[j], ref_ant_index]
-
-
- ### writing to the soltab
- outsoltab.setValues(outvals, selection)
- outsoltab.setValues(outweights, selection, weight = True)
- outsoltab.addHistory('Transferred solutions from ' + os.path.basename(refh5parm) + ' for the stations ' + str(stations_to_transfer).lstrip('[').rstrip(']') )
- outsoltab.flush()
-
- ### plotting new tables if provided
- if parset:
- if os.path.exists(parset):
- download = subprocess.Popen(['losoto', '-v', h5parmdb, parset], stdout=subprocess.PIPE)
- errorcode = download.wait()
- if errorcode != 0:
- logging.error('An error has occured while plotting.')
- data.close()
- refdata.close()
- return(1)
- else:
- logging.error('Parset file ' + parset + ' has not been found.')
- data.close()
- refdata.close()
- return(1)
-
- data.close()
- refdata.close()
- return(0)
-
-########################################################################
-if __name__ == '__main__':
- import argparse
- parser = argparse.ArgumentParser(description='Transfer solutions from a solution table of a reference h5parm to a solution table of a target h5parm.')
-
- parser.add_argument('h5parm', type=str,
- help='H5parm to which the solutions should be transferred.')
- parser.add_argument('--refh5parm', '--outh5parm', type=str,
- help='Name of the h5parm from which the solutions should be transferred.')
- parser.add_argument('--insolset', '--insolset', type=str, default='sol000',
- help='Name of the input h5parm solution set (default: sol000)')
- parser.add_argument('--outsolset', '--outsolset', type=str, default='sol000',
- help='Name of the output h5parm solution set (default: sol000)')
- parser.add_argument('--insoltab', '--insoltab', type=str, default='amplitude000',
- help='Name of the input h5parm solution set (default: amplitude000)')
- parser.add_argument('--outsoltab', '--outsoltab', type=str, default='amplitude000',
- help='Name of the output h5parm solution set (default: amplitude000)')
- parser.add_argument('--antenna', '--antenna', type=str, default='[FUSPID].*',
- help='Regular expression of antenna solutions to be transferred (default: [FUSPID].*)')
-
- args = parser.parse_args()
-
- format_stream = logging.Formatter("%(asctime)s\033[1m %(levelname)s:\033[0m %(message)s","%Y-%m-%d %H:%M:%S")
- format_file = logging.Formatter("%(asctime)s %(levelname)s: %(message)s","%Y-%m-%d %H:%M:%S")
- logging.root.setLevel(logging.INFO)
-
- log = logging.StreamHandler()
- log.setFormatter(format_stream)
- logging.root.addHandler(log)
-
- logging.info('Transferring solutions from ' + args.refh5parm + ' to ' + args.h5parm + '.')
- logging.info('Solutions will be transferred from soltab ' + str(args.insoltab) + ' to ' + str(args.outsoltab) + '.')
- main(args.h5parm, refh5parm = args.refh5parm, insolset=args.insolset, outsolset=args.outsolset, insoltab=args.insoltab, outsoltab=args.outsoltab, antenna = args.antenna, parset = None)