image.py 0000664 0021702 0000550 00000017012 11176425101 012206 0 ustar adrienne astro #
# Calling sytax:
#
# python image.py [--restart=]
#
# optional keyword "restart" can have the following values:
# - readin
# - clean
# - smooth
# - mask
# - reclean
#
# ie: python image.py --restart=clean
#
# The program will output various miriad directories from each task
# as well as a log (galaxy.image.log) that saves the miriad output.
#
# Steps run by the program:
#
# 1. readin
# -> *.map [keep]
# -> *.beam [keep]
#
# 2. clean
# -> *.clean
# -> *.restor
# -> *.cc.restor
# -> *.resid.restor
#
# 3. smooth
# -> *.restor.convol
# -> *.restor.convol.hanning
# -> *.blank
#
# 4. mask
# -> *.region [keep?]
# -> *.mask [keep]
#
# 5. reclean
# -> *.mclean [keep]
# -> *.mrestor [keep]
# -> *.cc.mrestor [keep]
# -> #.resid.mrestor [keep]
#
#
#
# To Add:
#
# 1. Figure out how to write to a log and have the input scroll
# (quickly) in the terminal as well, instead of in big chunks.
#
import os
import sys
import math
import immodules
# check if miriad is able to run:
if os.path.expandvars("$MIRBIN") == "$MIRBIN":
print ""
print "MIRIAD tasks cannot be called. " + \
"Source MIRRC or MIRRC.sh in the MIRIAD directory."
print ""
sys.exit(0)
# intialize variables for this scope
pfile = ""
restart = ""
noise = 0.0
lfchans = []
galaxy = ""
log = ""
# check calling sequence and set variables
if len(sys.argv) == 2:
pfile = sys.argv[1]
elif len(sys.argv) == 3:
pfile = sys.argv[1]
tmp = sys.argv[2]
tmp = tmp.split('=')
value = tmp[1]
value = value.lower()
if value == "readin" or value == "clean" or value == "smooth" or \
value == "mask" or value == "reclean":
restart = value
else:
print value + " is not a valid restart point."
print "Please choose: readin clean smooth mask reclean"
sys.exit(0)
else:
print "Calling sequence: python " + sys.argv[0] + \
" [--restart=]"
sys.exit(0)
# check that parameter file actually exists.
if not os.path.exists(pfile):
print "Parameter file " + pfile + " does not exist."
sys.exit(0)
f = open(pfile,mode="r")
# read file and set parameters
for line in f:
if line[0] == "#":
# ignore this line
continue
else:
this = line.split()
# search for = in list
try:
i = this.index("=")
except ValueError:
# this line does not have an =
continue
if len(this) != 3:
continue
param = this[0].lower()
value = this[2].lower()
# stupid switch statement
if param == "galaxy":
galaxy = value
elif param == "noise":
if value:
noise = float(value) # in mJy!
elif param == "lfchans":
tmp = value.split(',')
if len(tmp) >= 2:
lfchans = [int(tmp[0]), int(tmp[1])]
if len(tmp) >= 4:
lfchans.append(int(tmp[2]))
lfchans.append(int(tmp[3]))
if len(tmp) >= 6:
lfchans.append(int(tmp[4]))
lfchans.append(int(tmp[5]))
if len(tmp) >= 8:
lfchans.append(int(tmp[6]))
lfchans.append(int(tmp[7]))
# if len(tmp) == 2:
# lfchans = [int(tmp[0]), int(tmp[1])]
# if len(tmp) == 4:
# lfchans = [int(tmp[0]), int(tmp[1]), int(tmp[2]), int(tmp[3])]
# if len(tmp) ==
else:
print "Paramter " + param + " is not valid in " + pfile + "."
sys.exit(0)
if not(galaxy):
print "No galaxy parameter given."
print "Add this to the parameter file to continue:"
print "galaxy = "
sys.exit(0)
# set up dictionary to keep track of files that will be created.
createdFiles = {"readin": [galaxy + ".map", \
galaxy + ".beam"],
"clean": [galaxy + ".clean",\
galaxy + ".restor", \
galaxy + ".resid.restor", \
galaxy + ".restor.linmos", \
galaxy + ".cc.restor"], \
"smooth": [galaxy + ".cc.restor.convol", \
galaxy + ".cc.restor.convol.hanning"], \
"mask": [galaxy + ".region",\
galaxy + ".mask"], \
"reclean": [galaxy + ".mclean", \
galaxy + ".mrestor", \
galaxy + ".resid.mrestor",\
galaxy + ".restor.linmos", \
galaxy + ".cc.mrestor"] }
# remove old files.
save = restart
if restart == "readin":
for file in createdFiles[restart]:
if os.path.exists(file):
os.system("rm -r " + file)
restart = "clean"
if restart == "clean":
for file in createdFiles[restart]:
if os.path.exists(file):
os.system("rm -r " + file)
restart = "smooth"
if restart == "smooth":
for file in createdFiles[restart]:
if os.path.exists(file):
os.system("rm -r " + file)
restart = "mask"
if restart == "mask":
for file in createdFiles[restart]:
if os.path.exists(file):
os.system("rm -r " + file)
restart = "reclean"
if restart == "reclean":
for file in createdFiles[restart]:
if os.path.exists(file):
os.system("rm -r " + file)
restart = save
# set up the log.
if restart:
immodules.resetLog(galaxy, restart)
# now figure out where to start.
if os.path.exists(galaxy + ".map") and os.path.exists(galaxy + ".beam") and \
os.path.exists(galaxy + ".mask"):
# no masking needed - it's already done!
if not(lfchans):
lfchans = immodules.getLFChans(galaxy)
print "Line Free Channels are: " + str(lfchans)
if not(noise):
noise = immodules.getNoise(galaxy, lfchans=lfchans)
print "Noise is: " + str(noise) + " mJy."
if immodules.doStep(createdFiles["reclean"]):
immodules.clean(galaxy, 2.5*noise, mask = galaxy + ".mask")
# otherwise, start from the beginning and see if the files exist.
else:
if immodules.doStep(createdFiles["readin"]):
#remove any existing log file
if os.path.exists(galaxy + ".image.log"):
os.system("rm -r " + galaxy + ".image.log")
if os.path.exists(galaxy + ".isubim.1.fits") and \
os.path.exists(galaxy + ".bsubim.1.fits"):
immodules.readIn(galaxy)
else:
if not(os.path.exists(galaxy + ".isubim.1.fits")):
print "File "+ galaxy + \
".isubim.1.fits does not exist."
else:
print "File "+ galaxy + \
".bsubim.1.fits does not exist."
sys.exit(0)
# once the map is read in, set lfchans and noise if necessary.
if not(lfchans):
lfchans = immodules.getLFChans(galaxy)
print "Line Free Channels are: " + str(lfchans)
if not(noise):
noise = immodules.getNoise(galaxy, lfchans=lfchans)
print "Noise is: " + str(noise) + " mJy."
if immodules.doStep(createdFiles["clean"]):
immodules.clean(galaxy, 3*noise)
if immodules.doStep(createdFiles["smooth"]):
immodules.smooth(galaxy)
if immodules.doStep(createdFiles["mask"]):
immodules.mask(galaxy, noise)
if immodules.doStep(createdFiles["reclean"]):
immodules.clean(galaxy, 2.5*noise, mask = galaxy + ".mask")
print "Done!"
immodules.py 0000664 0021702 0000550 00000027626 11176425361 013146 0 ustar adrienne astro
import os
import sys
import math
import numpy
import pylab
import csv
# runs a miriad command
# inputs:
# list: a list of strings to be joined
# log [optional]: an output log (not yet implemented)
def run(list, log=""):
command = " ".join(list)
# this way doesn't work, but it's here in case I want to mess
# with it later.
#lines = []
#for line in os.popen(command, 'r', 1):
# lines.append(line)
# print line.split('\n')[0]
# print len(lines)
#if os.path.exists(log):
# file = open(log, 'a')
# file.writelines(lines)
# file.close()
# print the output to the screen and append to a file.
if log:
# uncomment this if you want to use tee instead of >>
# make sure to comment the second command if you do so.
#command = command + " | tee -a " + log
command = command + " >> " + log
os.system(command)
return
# removes a group of files if they exist
def remove(files):
for file in files:
if os.path.exists(file):
os.system("rm -r " + file)
# decides whether or not to run a reduction step
# inputs:
# files: a list of files created by the task in question
#
# output:
# boolean value: 1 if all the files exist and 0 if not all of them exist
def doStep(files):
skipStep = bool(1)
for file in files:
skipStep = skipStep & os.path.exists(file)
if not(skipStep):
for file in files:
remove(files)
return not(skipStep)
# find the size of the beam given a *.beam file.
# input:
# galaxy: base galaxy name
# output:
# bmaj: beam major axis fwhm
def getBeamsize(galaxy):
list = ['itemize']
list.append('in=' + galaxy + '.beam')
command = " ".join(list)
# open a pipe to read the miriad output and grab the major axis
pipe = os.popen(command + ' | grep bmaj')
output = pipe.read()
output = output.split()
bmaj = math.ceil(float(output[2]) * 206265)
return bmaj
# finds the line free channels from the history
#
# you can have up to 2 pairs of line free channels (ie beginning and end)
#
def getLFChans(galaxy):
# run prthis to find lf channel approximations from UVLSF inputs.
list = ["prthis"]
list.append("in=" + galaxy + ".map")
list.append("log=hist.tmp.txt")
list.append("task=uvlsf")
run(list)
pipe = os.popen("grep CHANSEL hist.tmp.txt")
output = pipe.read()
output = output.split("\n")
lfchans = [0,0]
# assumes that the last UVLSF entry is the one to use. Who knows.
for line in output:
tmp = line.split('=')
if tmp[0] == "UVLSF CHANSEL( 1)":
chans = tmp[1].split('/')[0].split(',')
if len(lfchans):
lfchans[0] = int(chans[0]) + 5
lfchans[1] = int(chans[1]) - 5
else:
lfchans.append(int(chans[0]) + 5)
lfchans.append(int(chans[1]) - 5)
if tmp[0] == "UVLSF CHANSEL( 2)":
chans = tmp[1].split('/')[0].split(',')
if len(lfchans) == 2:
lfchans.append(int(chans[0]) + 5)
lfchans.append(int(chans[1]) - 5)
elif len(lfchans) == 4:
lfchans[2] = int(chans[0]) + 5
lfchans[3] = int(chans[1]) - 5
os.system("rm hist.tmp.txt")
return lfchans
# determines the noise in a map.
def getNoise(galaxy, lfchans=[]):
# units are Jy/beam, colums are:
# 1: plane
# 2: Felocity
# 3: sum
# 4: mean
# 5: rms
# 6: max
# 7: min
# 8: npoints
# run imstat
list = ["imstat"]
list.append("in=" + galaxy +".map")
list.append("options=guaranteespaces,noheader")
list.append("log=imstat.tmp.txt")
if lfchans:
pass
else:
lfchans = getLFChans(galaxy)
region = ""
for i,chan in enumerate(lfchans):
if i % 2:
# it's odd
region = region + str(chan) + "\),"
else:
region = region + "image\(" + str(chan) + ","
# get rid of trailing comma
region = region[0:len(region)-1]
# if len(lfchans) == 2:
# list.append("region=image\(" + str(lfchans[0]) + "," + \
# str(lfchans[1]) + "\)")
# elif len(lfchans) == 4:
# list.append("region=image\(" + str(lfchans[0]) + "," + \
# str(lfchans[1]) + "\)," + \
# "image\(" + str(lfchans[2]) + "," + str(lfchans[3]) + "\)")
run(list)
# process imstat output
rms = []
chan = []
reader = csv.reader(open("imstat.tmp.txt","rb"))
reader_list = []
reader_list.extend(reader)
for line in reader_list:
tmp = line[0]
tmp = tmp.split()
if tmp[0] == "All":
pass
else:
rms.append(float(tmp[4]))
chan.append(int(tmp[0]))
pylab.plot(chan, rms)
pylab.show()
# find the median:
rms.sort()
n = len(rms)
if n & 1:
noise = rms[n // 2]
else:
noise = (rms[n // 2 - 1] + rms[n // 2]) / 2.0
os.system("rm imstat.tmp.txt")
# return the noise in mJy
return noise*1000
# read in the map and beam fits files and run velsw on the map to switch
# to the correct velocity units
def readIn(galaxy):
# write to the log:
f = open(galaxy + ".image.log", mode="a")
f.write("\n")
f.write("image.py READIN:\n")
f.close()
print "Running fits on the map!"
# run fits on the map
list = ['fits']
list.append('in=' + galaxy + '.isubim.1.fits')
list.append('op=xyin')
list.append('out=' + galaxy + '.map')
run(list, log=galaxy + ".image.log")
print "Running fits on the beam!"
# run fits on the beam
list = ['fits']
list.append('in=' + galaxy + '.bsubim.1.fits')
list.append('op=xyin')
list.append('out=' + galaxy + '.beam')
run(list, log=galaxy + ".image.log")
# run velsw on the map
print "Running velsw on the map!"
list = ['velsw']
list.append('in=' + galaxy + '.map')
list.append('axis=opt,bary')
run(list, log=galaxy + ".image.log")
# cleans to some depth with an optional mask. this task also restores
# the clean components/residuals/c+r map and corrects for the primary
# beam.
# Note: depth is in mJy
def clean(galaxy, depth, mask="", niters=0):
# for clarity later.
if mask:
m = "m"
# write to the log:
f = open(galaxy + ".image.log", mode="a")
f.write("\n")
f.write("image.py RECLEAN:\n")
f.close()
else:
m = ""
# write to the log:
f = open(galaxy + ".image.log", mode="a")
f.write("\n")
f.write("image.py CLEAN:\n")
f.close()
if m:
print "Running masked clean!"
else:
print "Running initial clean!"
# run clean
list = ['clean']
list.append('map=' + galaxy + '.map')
list.append('beam=' + galaxy + '.beam')
if niters:
list.append("niters=" + str(niters))
else:
list.append("niters=5000")
list.append('cutoff=' + str(depth * 0.001))
list.append("out=" + galaxy + "." + m + "clean")
if mask:
list.append('region=mask\(' + galaxy + '.mask\)')
run(list, log=galaxy + ".image.log")
print "Running restor for c+r map!"
# run restor
list = ["restor"]
list.append("map=" + galaxy + ".map")
list.append("beam=" + galaxy + ".beam")
list.append("model=" + galaxy + "." + m + "clean")
list.append("out=" + galaxy + "." + m + "restor")
run(list, log=galaxy + ".image.log")
# residual map
list.append("mode=resid")
list[4] = "out=" + galaxy + ".resid." + m + "restor"
print "Running restor for residual map!"
run(list, log=galaxy + ".image.log")
print "Running maths to find cc map!"
# subtract to get clean components...
list=["maths"]
list.append("exp=\<" + galaxy + "." + m + \
"restor\>-\<" + galaxy + ".resid." + m + "restor\>")
list.append("out=" + galaxy + ".cc." + m + "restor")
run(list)
# correct for the primary beam with linmos
list = ["linmos"]
list.append("in=" + galaxy + "." + m + "restor")
list.append("out=" + galaxy + "." + m + "restor.linmos")
run(list)
# smooth finds the beam size, convolves the .restor map with
# twice the beam, hanning smooths.
def smooth(galaxy):
# write to the log:
f = open(galaxy + ".image.log", mode="a")
f.write("\n")
f.write("image.py SMOOTH:\n")
f.close()
# get the beam size
bmaj = getBeamsize(galaxy)
print "Running convol!"
# convol
list = ['convol']
list.append('map=' + galaxy + '.cc.restor')
list.append('out=' + galaxy + '.cc.restor.convol')
list.append('fwhm=' + str(2 * bmaj) + ',' + str(2 * bmaj))
list.append('options=final')
run(list, log=galaxy + ".image.log")
print "Running hanning!"
# hanning
list = ['hanning']
list.append('in=' + galaxy + '.cc.restor.convol')
list.append('out=' + galaxy + '.cc.restor.convol.hanning')
list.append('width=3')
run(list, log=galaxy + ".image.log")
# mask runs cgcurs and then immask.
def mask(galaxy, noise):
# write to the log:
f = open(galaxy + ".image.log", mode="a")
f.write("\n")
f.write("image.py MASK:\n")
f.close()
# blank
list = ['maths']
list.append('exp=\<' + galaxy + '.cc.restor.convol.hanning\>')
list.append('mask=\<' + galaxy + '.cc.restor.convol.hanning\>.' +
'gt.'+ str(noise*0.001))
list.append('out=' + galaxy + '.mask')
print "Blanking cube!"
run(list, log=galaxy + ".image.log")
# cgcurs
list = ['cgcurs']
list.append('in=' + galaxy + '.mask')
list.append('type=pix')
list.append('device=/xwin')
list.append('range=0,0.00001')
list.append('nxy=1')
list.append("options=3pixel,3value,region")
print "Running cgcurs!"
run(list)
os.system('mv cgcurs.region ' + galaxy + '.region')
print "Creating the mask!"
# mask
list = ["immask"]
list.append("in=" + galaxy + ".mask")
list.append("region=@" + galaxy + ".region")
list.append("logic=and")
list.append("flag=true")
run(list, log=galaxy + ".image.log")
# sets up the log.
# this includes erasing lines for steps that will be repeated.
def resetLog(galaxy, restart=""):
if os.path.exists(galaxy + ".image.log"):
f = open(galaxy + ".image.log", mode='r')
lines = f.readlines()
f.close()
else:
return -1
if restart:
if restart == "readin":
# start over.
i = 0
elif restart == "clean":
try:
i = lines.index("image.py CLEAN:\n")
except ValueError:
# it's not found.
i=0
print '"image.py CLEAN" not found in log file.'
return
elif restart == "smooth":
try:
i = lines.index("image.py SMOOTH:\n")
except ValueError:
# it's not found.
i=0
print '"image.py SMOOTH" not found in log file.'
return
elif restart == "mask":
try:
i = lines.index("image.py MASK:\n")
except ValueError:
# it's not found.
i=0
print '"image.py MASK" not found in log file.'
return
elif restart == "reclean":
try:
i = lines.index("image.py RECLEAN:\n")
except ValueError:
# it's not found.
i = 0
print '"image.py RECLEAN" not found in log file.'
return
# remove the old file and write the new one.
os.system("rm " + galaxy + ".image.log")
if i:
f = open(galaxy + ".image.log", "w")
f.writelines(lines[:i])
f.close()
example.param 0000664 0021702 0000550 00000000325 11164742543 013237 0 ustar adrienne astro # Example parameter file for image.py
#
# Anything beginning with # will be skipped.
#
# noise is in mJy
# lfchans can have one or two pairs
#
#
galaxy = n4190.bcd05
#noise = 2.0
#lfchans = 1,2,3,4
#lfchans = 1,2
readme.txt 0000664 0021702 0000550 00000004127 11176425307 012563 0 ustar adrienne astro - image.py: main script that controls the imaging pipeline. Run using:
python image.py [--restart=value]
Value can be one of the following: readin, clean, smooth, mask,
reclean
- immodules.py: contains modules that actually image the uv data.
- example.param: an example parameter file. Comments are started with
#. All available parameters are commented out.
########################################################################
In order to initially run the script, you must first set up the
example.param file and make sure that your fits files have the correct
name (which they should if you output them from AIPS with the AIPS
Imaging scripts).
The .fits files should be named following this convention:
.isubim.1.fits
.bsubim.1.fits
The only way to change this (currently) is by editing
immodules.py. The output of the AIPS scripts will give you files named
(e.g. for ngc3109):
ngc3109-bcd50.isubim.1.fits (robust 5 image)
ngc3109-bcd50.bsubim.1.fits (robust 5 beam)
ngc3109-bcd05.isubim.1.fits (robust 0.5 image)
ngc3109-bcd05.bsubim.1.fits (robust 0.5 beam)
Next, set up your parameter file. Change the line:
galaxy =
to be:
galaxy = ngc3109-bcd50
for example, for ngc3109, robust 5. If I were you I would
change the name of the parameter file to be:
.param
but that is optional.
########################################################################
Updates:
v1.5: added support for 4 line-free channel ranges. [not yet released]
v1.4: fixed typo in createdFiles dictionary. [April 30, 2009]
v1.3: fixed bug where the linmos file was not deleted if
restarting. Changed output name of linmos file to be
.restor.linmos instead of .restor.pbcor. Added
this readme file to the tarball! [April 29, 2009]
v1.2: Added linmos to correct for primary beam in "clean" step. Added
checks to see if the .fits files actually exist. [April 27, 2009]
v1.1: don't remember.. changed cleaning method, I think. [April 7, 2009]
v1.0: original version. [April 1, 2009]