Source code for admit.at.Ingest_AT

""" .. _Ingest-at-api:

   **Ingest_AT** --- Ingests a (FITS) data cube.
   ---------------------------------------------

   This module defines the Ingest_AT class.
"""
import os, sys

import admit
from admit.AT import AT
from admit.Summary import SummaryEntry
import admit.util.bdp_types as bt
from admit.util.Image import Image
import admit.util.utils as utils
import admit.util.casautil as casautil
import admit.util.ImPlot  as ImPlot
from admit.bdp.SpwCube_BDP import SpwCube_BDP
from admit.bdp.Image_BDP   import Image_BDP
from admit.util.AdmitLogging import AdmitLogging as logging

import numpy as np
import math
try:
    import taskinit
    import casa
    from specsmooth import specsmooth
    from impbcor import impbcor
    from imtrans import imtrans
except:
    print "WARNING: No CASA; Ingest task cannot function."

import random

# @todo 
# - rel path should be to ../basename.fits
# = edge/box are the keywords to cut a cube
#   len(box) = 0       edge allowed
#              2       no edge allowed, they are z1,z2
#              4       edge allowed, they are blc,trc
#              6       no edge allowed
# - vlsr needs to be stored, in a veldef?   For this we need access to ASDM/Source.xml
#   and will need a small tool in util.py to parse the xml and return what we need
# - resolve the confusion between symlink= and copy=
# - allow LineCube instead of the default SpwCube
# - autobox?
# - if not a fits file, there is no smoothing
#   (obviously if it remains a symlink, of course we cannot change the data)
# - check if spectral reference is LSRK (e.g. some ALMA is TOPO, that's bad)
#   although some may not have it, e.g. test6503
#   SPECSYS = 'TOPO' or 'TOPOCENT' (casa 3.3.0)
#             'BARY'  (helio)
#   imhead->['reffreqtype']
#  - smooth and decimate option?
#  - vlsr=0.0 cannot be given technically?


[docs]class Ingest_AT(AT): """Ingest an image (cube) into ADMIT, normally to bootstrap a flow. See also :ref:`Ingest-AT-Design` for the design document. Ingest an image cube (usually FITS, but a CASA image or MIRIAD image are also natively supported by CASA) into CASA. Expect I/O penalties if you use FITS or MIRIAD because CASA images are tiled. A number of selections and corrections to the cube can be made, as specified by the keywords. Notably, a sub-cube can be taken out of the cube by trimming off spatial and spectral edges, a primary beam correction can be made as well. Internally ADMIT will store images as 4D CASA images, with any missing 3rd or 4th axis created redundantly (FREQ as axis 3, and POL as axis 4) **Keywords** **file**: string Input filename. Usually 'basename.fits' where 'basename' can be long, and can involve a directory hierarchy. The admit directory will be 'basename.admit', which will include the whole directory path if one was given. A symbolic link within ADMIT will then be used to resolve this as a local file. An absolute address is advised to be used if your working directory can change in a flow. If a CASA image is given, a symlink is used when no modifications (e.g. box=) are used, but this will cause uncertain integrity of the input BDP, as other clients could be modifying it. A MIRIAD image can also be given, but I/O (like with FITS) will be slower. [no default] **basename**: string New short basename (like an alias) of the output file. This is meant to override a possibly long basename in the input file. The **alias** keyword in the baseclass can still be used to create **basename-alias** type filenames. Warning: if you use a dash in your basename, there is a good risk this will be confused with the alias separator in the basename in a flow, as these are "basename-alias.extension". Default: empty string, basename inherited from the input file name. **pb**: string If given, this is the filename of the Primary Beam by which the input file needs to be multiplied to get a noise flat image for processing the ADMIT flow. The ALMA pipeline product of the Primary Beam should be in 'longname.pb.fits' where the flux flat cube is 'longname.pbcor.fits' Note: PB correction is slow. Default: empty string, no PB correction is done. **usepb**: boolean If True, the PB is actually used in an assumed flux flat input file to create a noise flat input BDP. If False, it is assumed the user has given a noise flat file, and the PB can be used downstream to compute proper PB corrected fluxes. Note that the PB is always stored as an 2D image, not as a cube. Default: True **box**: blc,tlc (a list of 2, 4 or 6 integers) Select a box region from the cube. For example box=[xmin,ymin,xmax,ymax], which takes all channels, or [xmin,ymin,zmin,xmax,ymax,zmax], which also selects a range in channels. You can also select just some channels, with box=[zmin,zmax]. As always, pixels and channels are 0 based in CASA. Arbitrary CASA regions are not implemented here, we only support a box/edge selection. **edge**: Z_start,Z_end (a list of 1 or 2 integers) You can use edge= to remove edge channels, e.g. if box= was not specified, or when only an XY box was given. If box contains 2 or 6 numbers, any edge specification would be ignored. If one number is given, the edge rejection is the same at the upper and lower end. Default: not used. **smooth**: [nx,ny,[nz]] (a list of 2 or 3 integers) You can convolve your cube spatially, and optionally spectrally as well, by supplying the number of pixels by which it is convolved. Spatially the FWHM is used. If nz is 1, a Hanning smooth is applied, else a boxcar of size nz is used. See also :ref:`Smooth-AT-api`, where a common beam can be computed or supplied. A future version should contain a decimation option. By default no smoothing is applied. **mask**: boolean If True, a mask needs to be created where the cube has 0's. This option is automatically bypassed if the input CASA image had a mask. [False] **vlsr**: float (km/s) VLSR of the source (km/s). If not set, the frequency at the center of the band is compared with the rest frequency from the header, which we call VLSRc. If the input file is not FITS, or header items are missing if the input file is CASA or MIRIAD already, unexpected things may happen. This VLSR (or VLSRc) is added to the ADMIT summary, which will be used downstream in the flow by other AT's (e.g. LineID) Default: -999999.99 (not set). **restfreq**: float (GHz) An alternative method providing the source VLSR would be to specify the true restfreq (f0) where the fits header has a 'fake' restfreq (f). This technique is sometimes used by the PI to avoid complex high-z doppler calculations and supply the redshifted line directly. In this case VLSR = c * (1-f/f0), in the radio definition, with z in the optical convention of course. We call this VLSRf. NOTE: clarify/check if the "1+z" velocity scale of the high-z object is correct. Default: -1.0 (method not used). Units must be GHz! **Input BDPs** None. The input is specific via the file= keyword. **Output BDPs** **SpwCube**: count: 1 The output spectral window cube. The convention is that the name of the BDP inherits from the basename of the input fits cube,and adding an extension "im". Each AT has a hidden keyword called alias=, use this keyword if you want to modify the cube name to "alias.im", and hence the BDP to "alias.im.bdp". Note this is an exception from the usual rule, where alias= is used to create a dashed-prefix to an extension, e.g. "x.alias-y". **Image_BDP**: 1 Output PB Map. If the input PB is a cube, the central channel (being representative for the avarage PB across the spectrum window) is used. New extension will be ".pb" Parameters ---------- keyval : dictionary, optional Keyword-value pairs, directly passed to the contructor for ease of assignment. Attributes ---------- _version : string Version ID for some future TBD use. Also should not be documented here, as underscore attributes are for internal usage only. """ #### DEPRECATED KEYWORDS BUT STILL ACTIVE IN CODE BY THEIR DEFAULT """ **symlink** : True/False: If True, A symlink is kept to the input file without any conversion (if that was needed) This is used in those cases where your whole flow can work with the fits file, without need to convert to a CASA image Setting to True, also disabled all other processing (mask/region/pbcor) Use with caution! [False] DEPRECATION """ def __init__(self,**keyval): keys = { 'file' : "", # fitsfile cube or map (or casa/miriad) 'basename': "", # override basename (useful for shorter names) 'pb' : "", # PB cube or map 'usepb' : True, # use PB, or was it just given for downstream 'mask' : True, # define a mask where data==0.0 if no mask present 'box' : [], # [] or z1,z2 or x1,y1,x2,y2 or x1,y1,z1,x2,y2,z2 'edge' : [], # [] or zl,zr - number of edge channels 'smooth' : [], # pixel smoothing size applied to data (can be slow) - see also Smooth_AT 'variflow': False, # requires manual sub-flow management for now 'vlsr' : -999999.99, # force a VLSR (see also LineID) 'restfreq': -1.0, # alternate VLSRf specification # 'symlink' : False, # # 'autobox' : False, # automatically cut away spatial and spectral slices that are masked # 'cbeam' : 0.5, # channel beam variation allowed in terms of pixel size to use median beam } AT.__init__(self,keys,keyval) self._version = "1.1.5" self.set_bdp_in() # no input BDP self.set_bdp_out([(SpwCube_BDP, 1), # one or two output BDPs (Image_BDP, 0), # optional PB if there was an pb= input ])
[docs] def summary(self): """Returns the summary dictionary from the AT, for merging into the ADMIT Summary object. Ingest_AT adds the following to ADMIT summary:: Key type Description -------- ------ ----------- fitsname string Pathless filename of FITS cube casaname string Pathless filename of CASA cube object string Object (or field) name naxis integer Number of axes naxisn integer size of axis n (n=1 to naxis0) crpix1 float Reference pixel axis 1 (n=1 to naxis) crvaln float axis value at CRPIX1 (n=1 to naxis) ctypen string axis type 1 (n=1 to naxis0 cdeltn float axis increment 1 (n=1 to naxis) cunitn string axis unit 1 (n=1 to naxis) equinox string equinox restfreq float rest frequency, Hz bmaj float beam major axis, radians bmin float beam minor axis, radians bpa float beam position angle, deg bunit string units of pixel values telescop string telescope name observer string observer name date-obs string date of observation datamax float maximum data value datamin float minimum data value badpixel float fraction of invalid pixels in the cube (a number between 0 and 1) vlsr float Object line-of-sight velocity (km/s) """ # @todo the master list is in $ADMIT/admit/summary_defs.tab # 1) we duplicate that here.... # 2) should it be with code? or in $ADMIT/etc ? if hasattr(self,"_summary"): return self._summary else: return {}
[docs] def run(self): # self._summary = {} # prepare to make a summary here dt = utils.Dtime("Ingest") # timer for debugging do_cbeam = True # enforce a common beam # pb = self.getkey('pb') do_pb = len(pb) > 0 use_pb = self.getkey("usepb") # create_mask = self.getkey('mask') # create a new mask ? box = self.getkey("box") # corners in Z, XY or XYZ edge = self.getkey("edge") # number of edge channels to remove restfreq = self.getkey("restfreq") # < 0 means not activated ckms = utils.c # 299792.458 km/s # smooth= could become deprecated, and/or include a decimation option to make it useful # again, Smooth_AT() does this also , at the cost of an extra cube to store smooth = self.getkey("smooth") # # vlsr = self.getkey("vlsr") # see also LineID, where this could be given again # first place a fits file in the admit project directory (symlink) # this is a bit involved, depending on if an absolute or relative path was # give to Ingest_AT(file=) fitsfile = self.getkey('file') if fitsfile[0] != os.sep: fitsfile = os.path.abspath(os.getcwd() + os.sep + fitsfile) logging.debug('FILE=%s' % fitsfile) if fitsfile[0] != os.sep: raise Exception,"Bad file=%s, expected absolute name",fitsfile # now determine if it could have been a CASA (or MIRIAD) image already # which we'll assume if it's a directory; this is natively supported by CASA # but there are tools where if you pass it a FITS or MIRIAD # MIRIAD is not recommended for serious work, especially big files, since there # is a performance penalty due to tiling. file_is_casa = casautil.iscasa(fitsfile) loc = fitsfile.rfind(os.sep) # find the '/' ffile0 = fitsfile[loc+1:] # basename.fits basename = self.getkey('basename') # (new) basename allowed (allow no dots?) if len(basename) == 0: basename = ffile0[:ffile0.rfind('.')] # basename logging.info("basename=%s" % basename) target = self.dir(ffile0) if not os.path.exists(target) : cmd = 'ln -s "%s" "%s"' % (fitsfile, target) logging.debug("CMD: %s" % cmd) os.system(cmd) readonly = False if file_is_casa: logging.debug("Assuming input %s is a CASA (or MIRIAD) image" % ffile0) bdpfile = self.mkext(basename,"im") if bdpfile == ffile0: logging.warning("No selections allowed on CASA image, since no alias was given") readonly = True b1 = SpwCube_BDP(bdpfile) self.addoutput(b1) b1.setkey("image", Image(images={bt.CASA:bdpfile})) # @todo b2 and PB? else: # construct the output name and construct the BDP based on the CASA image name # this also takes care of the behind the scenes alias= substitution bdpfile = self.mkext(basename,"im") if bdpfile == basename: raise Exception,"basename and bdpfile are the same, Ingest_AT needs a fix for this" b1 = SpwCube_BDP(bdpfile) self.addoutput(b1) if do_pb: print "doing the PB" bdpfile2 = self.mkext(basename,"pb") b2 = Image_BDP(bdpfile2) self.addoutput(b2) # @todo we should also set readonly=True if no box, no mask etc. and still an alias # that way it will speed up and not make a copy of the image ? # fni and fno are full (abspath) filenames, ready for CASA # fni is the same as fitsfile fni = self.dir(ffile0) fno = self.dir(bdpfile) if do_pb: fno2 = self.dir(bdpfile2) dt.tag("start") ia = taskinit.iatool() rg = taskinit.rgtool() if file_is_casa: ia.open(fni) else: if do_pb and use_pb: # @todo this needs a fix for the path for pb, only works if abs path is given # impbcor(im.fits,pb.fits,out.im,overwrite=True,mode='m') if False: # this may seem like a nice shortcut, to have the fits->casa conversion be done # internally in impbcor, but it's a terrible performance for big cubes. (tiling?) # we keep this code here, perhaps at some future time (mpi?) this performs better # @todo fno2 impbcor(fni,pb,fno,overwrite=True,mode='m') dt.tag("impbcor-1") else: # the better way is to convert FITS->CASA and then call impbcor() # the CPU savings are big, but I/O overhead can still be substantial _pbcor = utils.tmp_file('_pbcor_','.') _pb = utils.tmp_file('_pb_', '.') ia.fromfits(_pbcor,fni,overwrite=True) ia.fromfits(_pb, pb, overwrite=True) dt.tag("impbcor-1f") if False: impbcor(_pbcor,_pb,fno,overwrite=True,mode='m') # @todo fno2 utils.remove(_pbcor) utils.remove(_pb) dt.tag("impbcor-2") else: # immath appears to be even faster (2x in CPU) # https://bugs.nrao.edu/browse/CAS-8299 # @todo this needs to be confirmed that impbcor is now good to go (r36078) casa.immath([_pbcor,_pb],'evalexpr',fno,'IM0*IM1') dt.tag("immath") if True: # use the mean of all channels... faster may be to use the middle plane # barf; edge channels can be with fewer subfields in a mosaic ia.open(_pb) s=ia.summary() #ia1=taskinit.ia.moments(moments=[-1],drop=True,outfile=fno2) # fails for cont maps ia1=ia.collapse(outfile=fno2, function='mean', axes=2) # @todo ensure 2=freq axis ia1.done() ia.close() dt.tag("moments") utils.remove(_pbcor) utils.remove(_pb) dt.tag("impbcor-3") elif do_pb and not use_pb: # cheat case: PB was given, but not meant to be used # not implemented yet print "cheat case dummy PB not implemented yet" else: # no PB given if False: # re-running this was more consistently faster in wall clock time # note that zeroblanks=True will still keep the mask logging.debug("casa::ia.fromfits(%s) -> %s" % (fni,bdpfile)) ia.fromfits(fno,fni,overwrite=True) #taskinit.ia.fromfits(fno,fni,overwrite=True,zeroblanks=True) dt.tag("fromfits") else: # not working to extend 3D yet, but this would solve the impv() 3D problem logging.debug("casa::importfits(%s) -> %s" % (fni,bdpfile)) #casa.importfits(fni,fno,defaultaxes=True,defaultaxesvalues=[None,None,None,'I']) # possible bug: zeroblanks=True has no effect? casa.importfits(fni,fno,zeroblanks=True,overwrite=True) dt.tag("importfits") ia.close() ia.open(fno) if len(smooth) > 0: # smooth here, but Smooth_AT is another option # here we only allow pixel smoothing # spatial: gauss # spectral: boxcar/hanning (check for flux conservation) # is the boxcar wrong, not centered, but edged? # @todo CASA BUG: this will loose the object name (and maybe more?) from header, so VLSR lookup fails fnos = fno + '.smooth' ia.convolve2d(outfile=fnos, overwrite=True, pa='0deg', major='%gpix' % smooth[0], minor='%gpix' % smooth[1], type='gaussian') ia.close() srcname = casa.imhead(fno,mode="get",hdkey="object") # work around CASA bug #@todo use safer ia.rename() here. # https://casa.nrao.edu/docs/CasaRef/image.rename.html utils.rename(fnos,fno) casa.imhead(fno,mode="put",hdkey="object",hdvalue=srcname) # work around CASA bug dt.tag("convolve2d") if len(smooth) > 2 and smooth[2] > 0: if smooth[2] == 1: # @todo only 1 channel option specsmooth(fno,fnos,axis=2,function='hanning',dmethod="") else: # @todo may have the wrong center specsmooth(fno,fnos,axis=2,function='boxcar',dmethod="",width=smooth[2]) #@todo use safer ia.rename() here. # https://casa.nrao.edu/docs/CasaRef/image.rename.html utils.rename(fnos,fno) dt.tag("specsmooth") ia.open(fno) s = ia.summary() if len(s['shape']) != 4: logging.warning("Adding dummy STOKES-I axis") fnot = fno + '_4' ia2=ia.adddegaxes(stokes='I',outfile=fnot) ia.close() ia2.close() #@todo use safer ia.rename() here. # https://casa.nrao.edu/docs/CasaRef/image.rename.html utils.rename(fnot,fno) ia.open(fno) dt.tag("adddegaxes") else: logging.info("SHAPE: %s" % str(s['shape'])) s = ia.summary() dt.tag("summary-0") if s['hasmask'] and create_mask: logging.warning("no extra mask created because input image already had one") create_mask = False # if a box= or edge= was given, only a subset of the cube needs to be ingested # this however complicates PB correction later on if len(box) > 0 or len(edge) > 0: if readonly: raise Exception,"Cannot use box= or edge=, data is read-only, or use an basename/alias" if len(edge) == 1: edge.append(edge[0]) nx = s['shape'][0] ny = s['shape'][1] nz = s['shape'][2] logging.info("box=%s edge=%s processing with SHAPE: %s" % (str(box),str(edge),str(s['shape']))) if len(box) == 2: # select zrange if len(edge)>0: raise Exception,"Cannot use edge= when box=[z1,z2] is used" r1 = rg.box([0,0,box[0]] , [nx-1,ny-1,box[1]]) elif len(box) == 4: if len(edge) == 0: # select just an XY box r1 = rg.box([box[0],box[1]] , [box[2],box[3]]) elif len(edge) == 2: # select an XY box, but remove some edge channels r1 = rg.box([box[0],box[1],edge[0]] , [box[2],box[3],nz-edge[1]-1]) else: raise Exception,"Bad edge= for len(box)=4" elif len(box) == 6: # select an XYZ box r1 = rg.box([box[0],box[1],box[2]] , [box[3],box[4],box[5]]) elif len(edge) == 2: # remove some edge channels, but keep the whole XY box r1 = rg.box([0,0,edge[0]] , [nx-1,ny-1,nz-edge[1]-1]) else: raise Exception,"box=%s illegal" % box logging.debug("BOX/EDGE selection: %s %s" % (str(r1['blc']),str(r1['trc']))) #if taskinit.ia.isopen(): taskinit.ia.close() logging.info("SUBIMAGE") subimage = ia.subimage(region=r1,outfile=fno+'.box',overwrite=True) ia.close() ia.done() subimage.rename(fno,overwrite=True) subimage.close() subimage.done() ia.open(fno) dt.tag("subimage-1") else: # the whole cube is passed onto ADMIT if readonly and create_mask: raise Exception,"Cannot use mask=True, data read-only, or use an alias" if file_is_casa and not readonly: # @todo a miriad file - which should be read only - will also create a useless copy here if no alias used ia.subimage(overwrite=True,outfile=fno) ia.close() ia.open(fno) dt.tag("subimage-0") if create_mask: if readonly: raise Exception,"Cannot create mask, data read-only, or use an alias" # also check out the 'fromfits::zeroblanks = False' # calcmask() will overwrite any previous pixelmask #taskinit.ia.calcmask('mask("%s") && "%s" != 0.0' % (fno,fno)) ia.calcmask('"%s" != 0.0' % fno) dt.tag("mask") s = ia.summary() dt.tag("summary-1") # do a fast statistics (no median or robust) s0 = ia.statistics() dt.tag("statistics") if len(s0['npts']) == 0: raise Exception,"No statistics possible, are there valid data in this cube?" # There may be multiple beams per plane so we can't # rely on the BEAM's 'major', 'minor', 'positionangle' being present. # ia.commonbeam() is guaranteed to return beam parameters # if present if do_cbeam and 'perplanebeams' in s: # report on the beam extremities, need to loop over all, # first and last don't need to be extremes.... n = s['perplanebeams']['nChannels'] ab0 = '*0' bb0 = s['perplanebeams']['beams'][ab0]['*0'] bmaj0 = bb0['major']['value'] bmin0 = bb0['minor']['value'] beamd = 0.0 for i in range(n): ab1 = '*%d' % i bb1 = s['perplanebeams']['beams'][ab1]['*0'] bmaj1 = bb1['major']['value'] bmin1 = bb1['minor']['value'] beamd = max(beamd,abs(bmaj0-bmaj1),abs(bmin0-bmin1)) logging.warning("MAX-BEAMSPREAD %f" % (beamd)) # if True: logging.info("Applying a commonbeam from the median beam accross the band") # imhead is a bit slow; alternatively use ia.summary() at the half point for setrestoringbeam() h = casa.imhead(fno,mode='list') b = h['perplanebeams']['median area beam'] ia.setrestoringbeam(remove=True) ia.setrestoringbeam(beam=b) commonbeam = ia.commonbeam() else: # @todo : this will be VERY slow - code not finished, needs renaming etc. # this is however formally the better solution logging.warning("commmonbeam code not finished") cb = ia.commonbeam() ia.convolve2d(outfile='junk-common.im', major=cb['major'], minor=cb['minor'], pa=cb['pa'], targetres=True, overwrite=True) dt.tag('convolve2d') commonbeam = {} else: try: commonbeam = ia.commonbeam() except: nppb = 4.0 logging.warning("No synthesized beam found, faking one to prevent downstream problems: nppb=%f" % nppb) s = ia.summary() cdelt2 = abs(s['incr'][0]) * 180.0/math.pi*3600.0 bmaj = nppb * cdelt2 # use a nominal 4 points per (round) beam bmin = nppb * cdelt2 bpa = 0.0 ia.setrestoringbeam(major='%farcsec' % bmaj, minor='%farcsec' % bmin, pa='%fdeg' % bpa) commonbeam = {} logging.info("COMMONBEAM[%d] %s" % (len(commonbeam),str(commonbeam))) first_point = ia.getchunk(blc=[0,0,0,0],trc=[0,0,0,0],dropdeg=True) logging.debug("DATA0*: %s" % str(first_point)) ia.close() logging.info('BASICS: [shape] npts min max: %s %d %f %f' % (s['shape'],s0['npts'][0],s0['min'][0],s0['max'][0])) logging.info('S/N (all data): %f' % (s0['max'][0]/s0['rms'][0])) npix = 1 nx = s['shape'][0] ny = s['shape'][1] nz = s['shape'][2] for n in s['shape']: npix = npix * n ngood = int(s0['npts'][0]) fgood = (1.0*ngood)/npix logging.info('GOOD PIXELS: %d/%d (%f%% good or %f%% bad)' % (ngood,npix,100.0*fgood,100.0*(1 - fgood))) if s['hasmask']: logging.warning('MASKS: %s' % (str(s['masks']))) if not file_is_casa: b1.setkey("image", Image(images={bt.CASA:bdpfile})) if do_pb: b2.setkey("image", Image(images={bt.CASA:bdpfile2})) # cube sanity: needs to be either 4D or 2D. But p-p-v cube # alternative: ia.subimage(dropdeg = True) # see also: https://bugs.nrao.edu/browse/CAS-5406 shape = s['shape'] if len(shape)>3: if shape[3]>1: # @todo this happens when you ingest a fits or casa image which is ra-dec-pol-freq if nz > 1: msg = 'Ingest_AT: cannot deal with real 4D cubes yet' logging.critical(msg) raise Exception,msg else: # @todo this is not working yet when the input was a casa image, but ok when fits. go figure. fnot = fno + ".trans" if True: # this works #@todo use safer ia.rename() here. # https://casa.nrao.edu/docs/CasaRef/image.rename.html utils.rename(fno,fnot) imtrans(fnot,fno,"0132") utils.remove(fnot) else: # this does not work, what the heck imtrans(fno,fnot,"0132") #@todo use safer ia.rename() here. # https://casa.nrao.edu/docs/CasaRef/image.rename.html utils.rename(fnot,fno) nz = s['shape'][3] # get a new summary 's' ia.open(fno) s = ia.summary() ia.close() logging.warning("Using imtrans, with nz=%d, to fix axis ordering" % nz) dt.tag("imtrans4") # @todo ensure first two axes are position, followed by frequency elif len(shape)==3: # the current importfits() can do defaultaxes=True,defaultaxesvalues=['', '', '', 'I'] # but then appears to return a ra-dec-pol-freq cube # this branch probably never happens, since ia.fromfits() will # properly convert a 3D cube to 4D now !! # NO: when NAXIS=3 but various AXIS4's are present, that works. But not if it's pure 3D # @todo box= logging.warning("patching up a 3D to 4D cube") raise Exception,"SHOULD NEVER GET HERE" fnot = fno + ".trans" casa.importfits(fni,fnot,defaultaxes=True,defaultaxesvalues=['', '', '', 'I']) utils.remove(fno) # ieck imtrans(fnot,fno,"0132") utils.remove(fnot) dt.tag("imtrans3") logging.regression('CUBE: %g %g %g %d %d %d %f' % (s0['min'],s0['max'],s0['rms'],nx,ny,nz,100.0*(1 - fgood))) # if the cube has only 1 plane (e.g. continuum) , create a visual (png or so) # for 3D cubes, rely on something like CubeSum if nz == 1: implot = ImPlot(pmode=self._plot_mode,ptype=self._plot_type,abspath=self.dir()) implot.plotter(rasterfile=bdpfile,figname=bdpfile) # @todo needs to be registered for the BDP, right now we only have the plot # ia.summary() doesn't have this easily available, so run the more expensive imhead() h = casa.imhead(fno,mode='list') telescope = h['telescope'] # work around CASA's PIPELINE bug/feature? if 'OBJECT' is blank, try 'FIELD' srcname = h['object'] if srcname == ' ': logging.warning('FIELD used for OBJECT') srcname = casa.imhead(fno,mode='get',hdkey='field') if srcname == False: # if no FIELD either, we're doomed. yes, this did happen. srcname = 'Unknown' casa.imhead(fno,mode="put",hdkey="object",hdvalue=srcname) h['object'] = srcname logging.info('TELESCOPE: %s' % telescope) logging.info('OBJECT: %s' % srcname) logging.info('REFFREQTYPE: %s' % h['reffreqtype']) if h['reffreqtype'].find('TOPO')>=0: msg = 'Ingest_AT: cannot deal with cubes with TOPOCENTRIC frequencies yet - winging it' logging.warning(msg) #raise Exception,msg # Ensure beam parameters are available if there are multiple beams # If there is just one beam, then we are just overwriting the header # variables with their identical values. if len(commonbeam) != 0: h['beammajor'] = commonbeam['major'] h['beamminor'] = commonbeam['minor'] h['beampa'] = commonbeam['pa'] # cheat add some things that need to be passed to summary.... h['badpixel'] = 1.0-fgood if vlsr < -999998.0: vlsr = admit.VLSR().vlsr(h['object'].upper()) if vlsr == 0.0: vlsr = -999999.99 h['vlsr'] = vlsr logging.info("VLSR = %f (from source catalog)" % vlsr) taskargs = "file=" + fitsfile if create_mask == True: taskargs = taskargs + " mask=True" if len(box) > 0: taskargs = taskargs + " " + str(box) if len(edge) > 0: taskargs = taskargs + " " + str(edge) r2d = 57.29577951308232 logging.info("RA Axis 1: %f %f %f" % (h['crval1']*r2d,h['cdelt1']*r2d*3600.0,h['crpix1'])) logging.info("DEC Axis 2: %f %f %f" % (h['crval2']*r2d,h['cdelt2']*r2d*3600.0,h['crpix2'])) if nz > 1: # @todo check if this is really a freq axis (for ALMA it is, but...) t3 = h['ctype3'] df = h['cdelt3'] fc = h['crval3'] + (0.5*(float(shape[2])-1)-h['crpix3'])*df # center freq; 0 based pixels if 'restfreq' in h: fr = float(h['restfreq'][0]) # casa cheats, it may put 0 in here if FITS is missing it if fr == 0.0: fr = fc else: if vlsr < -999998.0: vlsr = (1-fc/fr)*ckms h['vlsr'] = vlsr else: fr = fc fw = df*float(shape[2]) print "PJT:",fr/1e9,fc/1e9,fw/1e9 dv = -df/fr*ckms logging.info("Freq Axis 3: %g %g %g" % (h['crval3']/1e9,h['cdelt3']/1e9,h['crpix3'])) logging.info("Cube Axis 3: type=%s velocity increment=%f km/s @ fc=%f fw=%f GHz" % (t3,dv,fc/1e9,fw/1e9)) # @todo sort out this restfreq/vlsr # report 'reffreqtype', 'restfreq' 'telescope' # if the fits file has ALTRVAL/ALTRPIX, this is lost in CASA? # but if you do fits->casa->fits , it's back in fits (with some obvious single precision loss of digits) # @todo ZSOURCE is the proposed VLSR slot in the fits header, but this has frame issues (it's also optical) # # Another method to get the vlsr is to override the restfreq (f0) with an AT keyword # and the 'restfreq' from the header (f) is then used to compute the vlsr: v = c (1 - f/f0) # if shape[2] > 1 and 'restfreq' in h: logging.info("RESTFREQ: %g %g %g" % (fr/1e9,h['restfreq'][0]/1e9,restfreq)) if shape[2] > 1: # v_radio of the center of the window w.r.t. restfreq vlsrc = ckms*(1-fc/fr) # @todo rel frame? vlsrw = dv*float(shape[2]) if restfreq > 0: vlsrf = ckms*(1-fr/restfreq/1e9) h['vlsr'] = vlsrf else: vlsrf = 0.0 logging.info("VLSRc = %f VLSRw = %f VLSRf = %f VLSR = %f" % (vlsrc, vlsrw, vlsrf, vlsr)) if h['vlsr'] == 0.0: # @todo! This fails if vlsr actually is zero. Need another magic number h['vlsr'] = vlsrc logging.warning("Warning: No VLSR found, substituting VLSRc = %f" % vlsrc) else: msg = 'Ingest_AT: missing RESTFREQ' print msg # @todo LINTRN is the ALMA keyword that designates the expected line transition in a spw self._summarize(fitsfile, bdpfile, h, shape, taskargs) dt.tag("done") dt.end()
def _summarize(self, fitsname, casaname, header, shape, taskargs): """Convenience function to populate dictionary for items to add to the ADMIT Summary. The contract function is self.summary(), called by AT() """ self._summary = {} self._summary['fitsname'] = SummaryEntry(fitsname) self._summary['casaname'] = SummaryEntry(casaname) # these are one-to-one match keywords easy = [ 'object' , 'equinox', 'observer', 'date-obs', 'datamax', 'datamin' , 'badpixel', 'vlsr', ] naxis = len(shape) self._summary['naxis'] = SummaryEntry(naxis) for i in range(naxis): j = i+1 jay = str(j) easy.append('crpix'+jay) easy.append('ctype'+jay) easy.append('crval'+jay) easy.append('cdelt'+jay) easy.append('cunit'+jay) self._summary['naxis'+jay] = SummaryEntry(int(shape[i])) # FITS is only 8 chars. if 'telescope' in header: self._summary['telescop'] = SummaryEntry(header['telescope']) if 'imtype' in header: self._summary['bunit'] = SummaryEntry(header['bunit']) for k in easy: if k in header: self._summary[k] = SummaryEntry(header[k]) if 'restfreq' in header: self._summary['restfreq'] = SummaryEntry(header['restfreq'][0]) # These are in imhead returned as dictionaries {'unit','value'} # so we have to munge them # convert beam parameters qa = taskinit.qatool() if 'beampa' in header: self._summary['bpa'] = SummaryEntry(qa.convert(header['beampa'],'deg')['value']) if 'beammajor' in header: self._summary['bmaj'] = SummaryEntry(qa.convert(header['beammajor'],'rad')['value']) if 'beamminor' in header: self._summary['bmin'] = SummaryEntry(qa.convert(header['beamminor'],'rad')['value']) # Now tag all summary items with task name and task ID. for k in self._summary: self._summary[k].setTaskname("Ingest_AT") self._summary[k].setTaskID(self.id(True)) self._summary[k].setTaskArgs(taskargs)