example.param 0000664 0021702 0000550 00000000436 11221215351 013224 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
#dniters = 5000 # for initial clean
#mniters = 500 # for masked clean
image.py 0000664 0021702 0000550 00000030175 11417671026 012222 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]
#
#
# 6. moments
# -> *.m0
# -> *.m1
# -> *.m2
# -> *.m-2
#
# 7. jvm
# -> *.spect.txt
# -> *.spect.ps
#
# 8. cd (column density)
# -> *.m0.cd
#
# 8. cleanup
# deletes unnecessary files and makes a tarball
#
# 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 = ""
dniters = 0
mniters = 0
# 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" or value == "moments" \
or value == "jvm" or value == "cleanup" or value == "recover" \
or value == "cd" or value == "spect":
restart = value
else:
print value + " is not a valid restart point."
print "Please choose: readin clean smooth mask reclean moments " + \
"jvm cleanup recover cd spect"
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 == "dniters":
if value:
dniters = int(value)
elif param == "mniters":
if value:
mniters = int(value)
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.hanning"], \
"mask": [galaxy + ".region",\
galaxy + ".mask"], \
"reclean": [galaxy + ".mclean", \
galaxy + ".mrestor", \
galaxy + ".resid.mrestor",\
galaxy + ".mrestor.linmos", \
galaxy + ".cc.mrestor"],
"moments": [galaxy + ".m0", \
galaxy + ".m1", \
galaxy + ".m2", \
galaxy + ".m-2"],\
"jvm": [galaxy + ".spect.ps",\
galaxy + ".spect.txt",\
galaxy + ".mrestor.jvm"],
"cd": [galaxy + ".m0.cd"],
"spect": [galaxy + ".pyspect.txt", \
galaxy + ".pyspect.eps"],
"keep": [galaxy + ".mask",\
galaxy + ".mrestor",\
galaxy + ".mrestor.linmos",\
galaxy + ".resid.mrestor",\
galaxy + ".cc.mrestor",\
galaxy + ".image.log",\
galaxy + ".m0",\
galaxy + ".m1",\
galaxy + ".m2",\
galaxy + ".m-2",\
galaxy + ".spect.ps",\
galaxy + ".spect.txt",\
galaxy + ".mrestor.jvm",\
galaxy + ".m0.cd", \
galaxy + ".pyspect.eps", \
galaxy + ".pyspect.txt"]
}
# run cleanup first if that is the restart value.
if restart == "cleanup":
if not(immodules.doStep(createdFiles["keep"])):
immodules.cleanup(galaxy, createdFiles, pfile)
print "Done!"
sys.exit(0)
else:
print "Not all files exist to clean up. Run all script steps first."
sys.exit(0)
if restart == "recover":
immodules.recover(galaxy, createdFiles["keep"])
sys.exit(0)
# run column density part - AMS added dec. 15 2009
#if restart == "cd":
# if os.path.exists(galaxy + ".beam") and os.path.exists(galaxy + ".m0"):
# if os.path.exists(galaxy + ".m0.cd"):
# os.system("rm -r " + galaxy + ".m0.cd")
# immodules.coldens(galaxy)
# sys.exit(0)
# else:
# print "Either beam or m0 file does not exist. Cannot run restart=cd."
# sys.exit(0)
# remove old files.
save = restart
if restart == "readin":
immodules.remove(createdFiles[restart])
restart = "clean"
if restart == "clean":
immodules.remove(createdFiles[restart])
restart = "smooth"
if restart == "smooth":
immodules.remove(createdFiles[restart])
restart = "mask"
if restart == "mask":
immodules.remove(createdFiles[restart])
restart = "reclean"
if restart == "reclean":
immodules.remove(createdFiles[restart])
restart = "jvm"
if restart == "jvm":
immodules.remove(createdFiles[restart])
restart = "moments"
if restart == "moments":
immodules.remove(createdFiles[restart])
restart = "cd"
if restart == "cd":
immodules.remove(createdFiles[restart])
restart = "spect"
if restart == "spect":
immodules.remove(createdFiles[restart])
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)
# append to parameter file for faster reading next time.
chanstr = ""
for chan in lfchans:
chanstr = chanstr + str(chan) + ","
chanstr = chanstr[0:len(chanstr)-1]
f = open(pfile, mode="a")
f.write("lfchans = " + chanstr + "\n")
f.close()
print "Line Free Channels are: " + str(lfchans)
if not(noise):
noise = immodules.getNoise(galaxy, lfchans=lfchans)
# append to parameter file for faster reading next time.
f = open(pfile, mode="a")
f.write("noise = " + str(noise) + "\n")
f.close()
print "Noise is: " + str(noise) + " mJy."
if immodules.doStep(createdFiles["reclean"]):
immodules.clean(galaxy, 2.5*noise, mask = galaxy + ".mask", \
niters = mniters)
# otherwise, start from the beginning and see if the files exist.
elif not(os.path.exists(galaxy + ".map")\
and os.path.exists(galaxy + ".resid.mrestor") \
and os.path.exists(galaxy + ".cc.mrestor") \
and os.path.exists(galaxy + ".mrestor") \
and os.path.exists(galaxy + ".mask")):
if immodules.doStep(createdFiles["readin"]):
#remove any existing log file
if os.path.exists(galaxy + ".image.log"):
os.system("rm -r " + galaxy + ".image.log")
# check for fits files
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. Add
# them to the log.
if not(lfchans):
lfchans = immodules.getLFChans(galaxy)
# append to parameter file for faster reading next time.
chanstr = ""
for chan in lfchans:
chanstr = chanstr + str(chan) + ","
chanstr = chanstr[0:len(chanstr)-1]
f = open(pfile, mode="a")
f.write("lfchans = " + chanstr + "\n")
f.close()
print "Line Free Channels are: " + str(lfchans)
# check if noise is set. find it and append if not.
if not(noise):
noise = immodules.getNoise(galaxy, lfchans=lfchans)
# append to parameter file for faster reading next time.
f = open(pfile, mode="a")
f.write("noise = " + str(noise) + "\n")
f.close()
print "Noise is: " + str(noise) + " mJy."
if immodules.doStep(createdFiles["clean"]):
immodules.clean(galaxy, 3*noise, niters = dniters)
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", \
niters = mniters)
if os.path.exists(galaxy + ".map")\
and os.path.exists(galaxy + ".resid.mrestor") \
and os.path.exists(galaxy + ".cc.mrestor") \
and os.path.exists(galaxy + ".mrestor") \
and os.path.exists(galaxy + ".mask"):
if immodules.doStep(createdFiles["jvm"]):
immodules.jvm(galaxy)
if immodules.doStep(createdFiles["moments"]):
immodules.moments(galaxy, noise)
if immodules.doStep(createdFiles["cd"]):
immodules.coldens(galaxy)
if immodules.doStep(createdFiles["spect"]):
immodules.spect(galaxy)
print "Done!"
immodules.py 0000664 0021702 0000550 00000072336 11423633161 013136 0 ustar adrienne astro import os
import sys
import math
import numpy
import pylab
import csv
import pyfits
import matplotlib.pyplot as plt
import matplotlib
from matplotlib.ticker import MaxNLocator
# 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:
# list of the form [bmaj, bmin] (fwhm of major axis and minor axis in arcsec)
def getdV(galaxy):
list = ['itemize']
list.append('in=' + galaxy + '.map')
command = " ".join(list)
# open a pipe to read the miriad output and grab the major axis
pipe = os.popen(command + ' | grep cdelt3')
output = pipe.read()
output = output.split()
dv = math.fabs(float(output[2]))
return dv
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)
pipe = os.popen(command + " | grep bmin")
output = pipe.read()
output = output.split()
bmin = math.ceil(float(output[2]) * 206265)
return [bmaj,bmin]
# 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])
lfchans[1] = int(chans[1])
else:
lfchans.append(int(chans[0]))
lfchans.append(int(chans[1]))
if tmp[0] == "UVLSF CHANSEL( 2)":
chans = tmp[1].split('/')[0].split(',')
if len(lfchans) == 2:
lfchans.append(int(chans[0]))
lfchans.append(int(chans[1]))
elif len(lfchans) == 4:
lfchans[2] = int(chans[0])
lfchans[3] = int(chans[1])
if tmp[0] == "UVLSF CHANSEL( 3)":
chans = tmp[1].split('/')[0].split(',')
if len(lfchans) == 4:
lfchans.append(int(chans[0]))
lfchans.append(int(chans[1]))
elif len(lfchans) == 6:
lfchans[4] = int(chans[0])
lfchans[5] = int(chans[1])
if tmp[0] == "UVLSF CHANSEL( 4)":
chans = tmp[1].split('/')[0].split(',')
if len(lfchans) == 6:
lfchans.append(int(chans[0]))
lfchans.append(int(chans[1]))
elif len(lfchans) == 8:
lfchans[6] = int(chans[0])
lfchans[7] = int(chans[1])
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]) + "\)")
list.append("region="+region)
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=50000")
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)[0]
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")
os.system("rm -r " + galaxy + ".cc.restor.convol")
# 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")
def moments(galaxy, noise):
# write to the log:
f = open(galaxy + ".image.log", mode="a")
f.write("\n")
f.write("image.py MOMENTS:\n")
f.close()
# mask the restored cube - should already be done
#list = ["immask"]
#list.append("in=" + galaxy + ".mrestor.jvm")
#list.append("region=mask\(" + galaxy + ".mask\)")
#list.append("logic=and")
#list.append("flag=true")
#run(list)
# run zeroth moment
list = ["moment"]
list.append("in=" + galaxy + ".mrestor.jvm")
list.append("out=" + galaxy + ".m0")
list.append("mom=0")
run(list)
# correct m0 for the primary beam with linmos - already done
#list = ["linmos"]
#list.append("in=" + galaxy + ".m0.tmp")
#list.append("out=" + galaxy + ".m0")
#run(list)
# remove un-corrected m0
#os.system("rm -r " + galaxy + ".m0.tmp")
# first moment
list = ["moment"]
list.append("in=" + galaxy + ".mrestor.jvm")
list.append("out=" + galaxy + ".m1.tmp")
list.append("mom=1")
run(list)
# second moment
list[2] = "out=" + galaxy + ".m2.tmp"
list[3] = "mom=2"
run(list)
# -2 moment - this will be converted to kelvin
list[2] = "out=" + galaxy + ".m-2.tmp"
list[3] = "mom=-2"
run(list)
# convolve the m0 map with 2x the beam size
beam = getBeamsize(galaxy)
bmaj = beam[0]
bmin = beam[1]
# convol
list = ['convol']
list.append('map=' + galaxy + '.m0')
list.append('out=' + galaxy + '.m0.convol')
list.append('fwhm=' + str(2 * bmaj) + ',' + str(2 * bmaj))
list.append('options=final')
run(list)
# convert to column density
list = ["maths"]
list.append('exp=\<' + galaxy + '.m0.convol\>*1.222e6/'+ \
str(4 * bmaj * bmaj) + "/1.420/1.420*1.8224e18")
# 4 * bmaj * bmaj is because we convolved with a larger beam.
list.append('out=' + galaxy + '.m0.convol.cd')
run(list)
# get the units of the file right
list = ["puthd"]
list.append("in=" + galaxy + ".m0.convol.cd/bunit")
list.append("value=cm^-2")
run(list)
# mask mom1 at 5e20 column density
list = ["maths"]
list.append("exp=\<" + galaxy + ".m1.tmp\>")
list.append("mask=\<" + galaxy + ".m0.convol.cd\>.ge.5e20")
list.append("out=" + galaxy + ".m1")
run(list)
# mask mom2 at 5e20 column density
list = ["maths"]
list.append("exp=\<" + galaxy + ".m2.tmp\>")
list.append("mask=\<" + galaxy + ".m0.convol.cd\>.ge.5e20")
list.append("out=" + galaxy + ".m2")
run(list)
# remove .tmp moment maps
os.system("rm -r " + galaxy + ".m?.tmp")
# remove m0.convol files
os.system("rm -r " + galaxy + ".m0.convol*")
# convert m-2 to kelvin units
list = ["maths"]
list.append("exp=\<" + galaxy + ".m-2.tmp\>*1.222e6/1.4/1.4/" + \
str(bmaj) + "/" + str(bmin))
list.append("out=" + galaxy + ".m-2")
run(list)
# put the proper units in the header
list = ["puthd"]
list.append("in=" + galaxy + ".m-2/bunit")
list.append("value=K")
run(list)
os.system("rm -r " + galaxy + ".m-2.tmp")
# delete the mask from the .mrestor file
#os.system("rm " + galaxy + ".mrestor./mask")
# creates spectrum with jvm correction.
def jvm(galaxy):
# eD = C + eR
# e(D - R) = C
# e = C / (D - R)
#
# G = C*D / (D - R)
# write to the log:
f = open(galaxy + ".image.log", mode="a")
f.write("\n")
f.write("image.py JVM:\n")
f.close()
# remove linmos files
if os.path.exists(galaxy + ".map.linmos"):
os.system("rm -r " + galaxy + ".map.linmos")
if os.path.exists(galaxy + ".resid.mrestor.linmos"):
os.system("rm -r " + galaxy + ".resid.mrestor.linmos")
if os.path.exists(galaxy + ".cc.mrestor.linmos"):
os.system("rm -r " + galaxy + ".cc.mrestor.linmos")
# correct all for the primary beam
list = ["linmos"]
list.append("in=" + galaxy + ".map")
list.append("out=" + galaxy + ".map.linmos")
run(list)
list[1] = "in=" + galaxy + ".cc.mrestor"
list[2] = "out=" + galaxy + ".cc.mrestor.linmos"
run(list)
list[1] = "in=" + galaxy + ".resid.mrestor"
list[2] = "out=" + galaxy + ".resid.mrestor.linmos"
run(list)
# mask the resid, cc, dirty cubes.
list = ["immask"]
list.append("in=" + galaxy + ".resid.mrestor.linmos")
list.append("region=mask\(" + galaxy + ".mask\)")
list.append("logic=and")
list.append("flag=true")
run(list)
list[1] = "in=" + galaxy + ".map.linmos"
run(list)
list[1] = "in=" + galaxy + ".cc.mrestor.linmos"
run(list)
# run these through imstat.
list = ["imstat"]
list.append("in=" + galaxy +".map.linmos")
list.append("options=guaranteespaces,noheader")
list.append("log=imstat.dirty.tmp.txt")
run(list)
list = ["imstat"]
list.append("in=" + galaxy +".cc.mrestor.linmos")
list.append("options=guaranteespaces,noheader")
list.append("log=imstat.cc.tmp.txt")
run(list)
list = ["imstat"]
list.append("in=" + galaxy +".resid.mrestor.linmos")
list.append("options=guaranteespaces,noheader")
list.append("log=imstat.resid.tmp.txt")
run(list)
# read in imstat output
npix = []
felocity = []
flux = []
epsilon = []
c = []
d = []
r = []
plane = []
reader = csv.reader(open("imstat.resid.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:
r.append(float(tmp[3]))
# npix and felocity should be the same for all
# maps at this channel.
npix.append(int(tmp[7]))
felocity.append(float(tmp[1]))
plane.append(float(tmp[0]))
reader = csv.reader(open("imstat.cc.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:
c.append(float(tmp[3]))
reader = csv.reader(open("imstat.dirty.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:
d.append(float(tmp[3]))
# figure out constants
beam = getBeamsize(galaxy)
bmaj = beam[0]
bmin = beam[1]
npixPerBeam = 1.133 * bmaj * bmin / 1.5 / 1.5 # cellsize is 1.5"
# get dv - channel width
dv = getdV(galaxy)
# loop through and figure out total flux and integrated flux
fint = 0
for i in range(len(d)):
d[i] = d[i] * npix[i] / npixPerBeam
c[i] = c[i] * npix[i] / npixPerBeam
r[i] = r[i] * npix[i] / npixPerBeam
if math.fabs(d[i] - r[i]) < 0.00001:
epsilon.append(0)
flux.append(0)
else:
epsilon.append(c[i] / (d[i] - r[i]))
flux.append(epsilon[i] * d[i])
fint = fint + flux[i] * dv
pylab.plot(felocity, flux)
pylab.show()
# write these out to a file
f = open(galaxy + ".spect.txt","w")
f.write("# Columns: plane, felocity, flux, epsilon, npix, d, c, r\n")
for i in range(len(d)):
f.write(str(plane[i]) + "\t" + str(felocity[i]) + "\t" + \
str(round(flux[i]*1000000)/1000000) + "\t" + \
str(round(epsilon[i]*1000000)/1000000) + "\t" + \
str(round(npix[i]*1000000)/1000000) + "\t" + \
str(round(d[i]*1000000)/1000000) + "\t" + \
str(round(c[i]*1000000)/1000000) + "\t" + \
str(round(r[i]*1000000)/1000000) + "\n")
f.write("# Total Integrated Flux, corrected: " + \
str(math.ceil(fint*1000)/1000) + " Jy km/s\n")
f.close()
os.system("cp " + galaxy + ".spect.txt data.dat")
os.system("$red/../SCRIPTS/jvm.sh data.dat")
os.system("mv galaxy.ps " + galaxy + ".spect.ps")
os.system("rm data.dat")
# clean up imstat files
os.system("rm imstat.*.tmp.txt")
# write out files to fits and run idl jvm correction program
list = ["fits"]
list.append("in=" + galaxy + ".cc.mrestor.linmos")
list.append("op=xyout")
list.append("out=" + galaxy + ".tmp.cc.fits")
run(list)
list = ["fits"]
list.append("in=" + galaxy + ".resid.mrestor.linmos")
list.append("op=xyout")
list.append("out=" + galaxy + ".tmp.resid.fits")
run(list)
#idl jvm call
os.system("idl<*1.222e6/'+ \
str(bmaj * bmin) + "/1.420/1.420*1.8224e18")
list.append('out=' + galaxy + '.m0.cd')
run(list)
# get the units of the file right
list = ["puthd"]
list.append("in=" + galaxy + ".m0.cd/bunit")
list.append("value=cm^-2")
run(list)
def spect(galaxy):
# get rid of pyspect files if they exist
remove([galaxy + ".pyspect.txt", galaxy + ".pyspect.ps"])
# make sure the .jvm.fits file exists
if not os.path.exists(galaxy + ".mrestor.jvm.fits"):
if os.path.exists(galaxy + ".mrestor.jvm"):
list = ["fits"]
list.append("in=" + galaxy + ".mrestor.jvm")
list.append("out=" + galaxy + ".mrestor.jvm.fits")
list.append("op=xyout")
run(list)
# read it in to python
cube = pyfits.getdata(galaxy + ".mrestor.jvm.fits")
header = pyfits.getheader(galaxy + ".mrestor.jvm.fits")
# deal with nans - just set them to 0.
cube[numpy.isnan(cube).nonzero()] = 0
# calculate the spectra - don't need npix and the mean, just the sum:
# * npix * beam / pix = sum * beam / pix
bmaj, bmin = getBeamsize(galaxy)
spect = numpy.sum(numpy.sum(cube, axis=1), axis=1)
spect = spect * (1.5**2) / (1.133 * bmaj * bmin)
# get the velocity axis
velocity = numpy.arange(header['naxis3'])
slope = header['cdelt3']
intercept = header['crval3'] - slope * (header['crpix3']-1)
velocity = slope * velocity + intercept
velocity = velocity / 1.0e3
ind = (spect > 0).nonzero()
#matplotlib.rc("xtick", labelsize=12)
#matplotlib.rc("ytick", labelsize=12)
matplotlib.rc('font', family="serif")
matplotlib.rc('font', size=24)
matplotlib.rc('axes', linewidth=2)
matplotlib.rc("xtick.major", size = 10)
matplotlib.rc("xtick.major", pad = 10)
matplotlib.rc("ytick.major", size = 12)
matplotlib.rc("ytick.major", pad = 10)
# now plot it
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(velocity[ind], spect[ind], 'k', linewidth=2)
ax.xaxis.set_major_locator(MaxNLocator(5))
ax.yaxis.set_major_locator(MaxNLocator(5, steps=[0.1]))
plt.text(0.9, 0.85, ((galaxy.partition("-"))[0]).upper(), \
transform=ax.transAxes, horizontalalignment="right")
#plt.title(((galaxy.partition("-"))[0]).upper())
plt.subplots_adjust(bottom=0.15, left=0.17)
ax2 = plt.gca()
for line in ax.xaxis.get_ticklines():
line.set_markeredgewidth(2)
for line in ax.yaxis.get_ticklines():
line.set_markeredgewidth(2)
plt.xlabel("Velocity (km/s)")
plt.ylabel("Flux Density (Jy)")
plt.savefig(galaxy + ".pyspect.eps")
# write the spectrum out to a file
file = open(galaxy + ".pyspect.txt", 'w')
file.write("# channel \t velocity \t fluxdensity\n")
file.write("# \t (km/s) \t (Jy)\n")
for chan, vel, flux in zip(numpy.arange((spect.shape)[0]), \
velocity,
spect):
file.write(str(chan) + "\t" + \
str(round(vel * 100) / 100.0) + "\t" + \
str(flux) + "\n")
file.close()
# deletes unnecessary miriad files, writes others out to fits, and makes a tarball.
def cleanup(galaxy, createdFiles, paramFile):
# check that tar file doesn't exist.
if os.path.exists(galaxy + ".image.tar"):
print "File " + galaxy + ".image.tar exists. Please delete " + \
"before continuing."
return
#write out "keep" files to fits
for file in createdFiles["keep"]:
if file == galaxy + ".spect.ps" or file == galaxy + ".spect.txt" or \
file == galaxy + ".image.log":
pass
else:
if os.path.exists(file + ".fits"):
os.system("rm " + file + ".fits")
list = ["fits"]
list.append("in=" + file)
list.append("out=" + file + ".fits")
list.append("op=xyout")
run(list)
# tar up the fits files
command = "tar --create --file=" + galaxy + ".image.tar " + paramFile
for file in createdFiles["keep"]:
if file == galaxy + ".spect.ps" or file == galaxy + ".spect.txt" or \
file == galaxy + ".image.log":
command = command + " " + file
else:
command = command + " " + file + ".fits"
command = command + " " + galaxy + ".isubim.1.fits " + galaxy + \
".bsubim.1.fits"
os.system(command)
# remove all the miriad files (including the "keep" ones because we have
# them in fits format)
remove(createdFiles["readin"])
remove(createdFiles["clean"])
remove(createdFiles["smooth"])
remove(createdFiles["mask"])
remove(createdFiles["reclean"])
remove(createdFiles["moments"])
def recover(galaxy, keepFiles):
if os.path.exists(galaxy + ".image.tar"):
os.system("tar -xvf " + galaxy + ".image.tar")
else:
print "File " + galaxy + ".image.tar does not exist."
return
if os.path.exists(galaxy + ".isubim.1.fits") and \
os.path.exists(galaxy + ".bsubim.1.fits"):
readin(galaxy)
else:
print "Files " + galaxy + ".isubim.1.fits and " + galaxy + \
".bsubim.1.fits are needed to recover setup."
return
for file in keepFiles:
if file == galaxy + ".spect.ps" or file == galaxy + ".spect.txt" \
or file == galaxy + ".image.log":
pass
else:
list = ["fits"]
list.append("in=" + file + ".fits")
list.append("out=" + file)
list.append("op=xyin")
run(list)
# 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:
i = 0
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
elif restart == "moments":
try:
i = lines.index("image.py MOMENTS:\n")
except ValueError:
# it's not found
i = 0
print '"image.py MOMENTS" not found in log file.'
return
elif restart == "jvm":
try:
i = lines.index("image.py JVM:\n")
except ValueError:
# it's not found
i = 0
print '"image.py JVM" 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()
jvm.pro 0000664 0021702 0000550 00000001450 11417150545 012074 0 ustar adrienne astro pro jvm, galaxy
readcol, galaxy + ".spect.txt", plane, felocity, flux, eps, /silent
cc = mrdfits(galaxy + ".tmp.cc.fits",0,header,/silent)
resid = mrdfits(galaxy + ".tmp.resid.fits",/silent)
; get dimensions of cc and resid arrays
dim = size(cc,/dimensions)
; set up jvm-corrected array dimensions
jvm = fltarr(dim[0],dim[1],dim[2])
; idl (arrays start at 0) is one off from miriad (arrays start at 1)
for i=1,dim[2] do begin
ind = where(i eq plane, count)
if count gt 0 then begin
jvm[*,*,i-1] = (fltarr(dim[0],dim[1]) + (eps[ind])[0]) $
* resid[*,*,i-1] + cc[*,*,i-1]
;stop
endif
endfor
; unallocate arrays that don't matter anymore
cc=0
resid=0
mwrfits, reform(jvm, dim[0], dim[1], dim[2], 1), $
galaxy + ".mrestor.jvm.fits", header, /create, /silent
;stop
end
jvm.sh 0000775 0021702 0000550 00000000403 11221226663 011704 0 ustar adrienne astro #!/bin/sh
file=$1
# run SM
sm -s -q << EOF
device postencap galaxy.ps
data $file
read {x 2 y 3}
erase limits x y
expand 1.1
box
expand 1
ptype 8 3
points x y
ltype 0
connect x y
expand 1.5
xlabel Felocity [km/s]
ylabel Flux [Jy]
expand 1
hardcopy
EOF
miriad2fits.py 0000664 0021702 0000550 00000001517 11221215351 013337 0 ustar adrienne astro # This converts miriad files to fits files. It will not delete any existing
# fits files, so beware.
#
# If this is a feature people want, I can change it to delete the fits files
# before making them again.
import os
# Check if miriad tasks can be 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)
# get the files that have a /image subdirectory. These should be miriad files.
pipe = os.popen("ls */image")
files = pipe.read().split()
for file in files:
try:
# take off the /image suffix
i = file.index('/')
file = file[0:i]
# write them out to fits
command = "fits in=" + file + " out=" + file + ".fits op=xyout"
print command
os.system(command)
except ValueError:
pass readme.txt 0000664 0021702 0000550 00000011761 11423634133 012557 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, moments, jvm, cleanup
- 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.
- jvm.sh: shell script that runs sm to make a ps file of the galaxy
spectrum (jvm-corrected of course)
You *must* put this in the $red/../SCRIPTS/ directory (or
make a symlink)!!
########################################################################
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. The parameter file must be:
.param (ie ngc3109-bcd50.param)
in order for the script to automatically save the parameter file
in the final tarball.
########################################################################
CLEANUP OPTION
The cleanup option is slightly different from the regular script. This
will check to see if all the files created by the script exist. If they
do, it will write out the files we will keep to fits format and tar them
up. Then it will delete *all* the miriad files. (We still have the ones
we want to save in fits format.) The script will then exit.
########################################################################
Features to add:
- ask to continue pre-smooth/masking?
########################################################################
Updates:
v1.0: original version. [April 1, 2009]
v1.1: don't remember.. changed cleaning method, I think. [April 7, 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.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.4: fixed typo in createdFiles dictionary. [April 30, 2009]
v1.5: fixed bug in getNoise. [May 4, 2009]
v1.6: added support for 4 line-free channel ranges. getLFChans no
longer adds/subtracts 5 in case of CVEL. It also displays the
noise as a function of channel number so you can make sure it's
working how you want it to. Just close the box to continue
running the script. [May 5, 2009]
v2.0: added capability to make moment maps. This will probably be
edited as it clips on a fixed column density, not any
threshold determined by the noise. [May 5, 2009]
v2.1: Added jvm capability.
Added miriad2fits.py script, which writes out all miriad (image)
files to fits format.
Fixed minor bug in image.py script that does not redo moment
maps if restarting at an earlier step.
Added linmos correction on m0 map.
Added "cleanup" option for restart. This deletes old files and
makes a tarball of the ones generated by this script to keep.
[May 14, 2009]
v2.2: Added check in JVM script to prevent denominator of 0.
[May 18, 2009]
v2.3: Added restart=recover option. Changed number of initial clean
iterations to 50000 (from 5000).
v2.4: Fixed typo in getBeamsize (from 2062065 to 206265). [June 9, 2009]
v2.5: Lots of updates [June 23, 2009]:
- Changed image.py so jvm and moments can be run with minimal
files
- Changed order of jvm and moments
- Changed jvm to create corrected cube in idl (jvm.pro)
- Changed moments to make moment maps from jvm-corrected cube
- Added capability to specify number of iterations in param file
(dniters, mniters for dirty cube and masked cube)
v2.6,2.7: Minor bug fixes
v2.8: - update to jvm.sh - change columns that SM is plotting
- some minor script running problems fixed
- added miriad2fits.py standalone program
v2.9: - added capability to create column density map from m0 (restart=cd)
v2.10: - fixed moment map/column density map bug
(bmaj * bmaj to bmaj * bmin)
v2.11: - ...??
v2.12: - fixed spectra - new CORRECT spectra files are .pyspect.eps and .pyspect.txt