From b7b01b512efa844af0744231a556d2fa3717b7df Mon Sep 17 00:00:00 2001 From: Pei-Ying Hsieh <95929782+pyhsiehATalma@users.noreply.github.com> Date: Fri, 12 May 2023 22:17:55 +0800 Subject: [PATCH 1/2] script to combine 7m and TP cubes This script does feather and mosaic cubes of 7m and TP data. --- aces/imaging/TP-ACA-feather.py | 200 +++++++++++++++++++++++++++++++++ 1 file changed, 200 insertions(+) create mode 100644 aces/imaging/TP-ACA-feather.py diff --git a/aces/imaging/TP-ACA-feather.py b/aces/imaging/TP-ACA-feather.py new file mode 100644 index 00000000..f78eb9c6 --- /dev/null +++ b/aces/imaging/TP-ACA-feather.py @@ -0,0 +1,200 @@ +# This script is to feather 7m and TP cubes. + + +# A template CMZ cube for mosaic +#templateIM='/almalustre/home/pyhsieh/ESO/spw24/template/CS2-1-temlate-whole-band.image' # whole band +#templateIM='/almalustre/home/pyhsieh/ESO/spw24/template/CS2-1-temlate-narrow-band.image' # only 0~1500 channel +templateIM='CS2-1-temlate-whole-band.image' + +import glob +linename=glob.glob('*.image.pbcor') +linewt=glob.glob('*.weight') + + +fieldnum=38 + +# update rest frequency of 12m cubes (line and weight data) with CS(2-1) line + +for i in range(len(linename)): + imhead( imagename=linename[i], + mode='put', + hdkey='restfreq',hdvalue='97.980953GHz',verbose=False) + +for i in range(len(linewt)): + imhead( imagename=linewt[i], + mode='put', + hdkey='restfreq',hdvalue='97.980953GHz',verbose=False) + + + +# reading spectral line and weight cubes of 7m and TP from dict.txt +# in dict-aca.txt, adding the pairs of 7m and TP mous id + +daca={} +dtp={} +with open("dict-aca.txt") as f: + for line in f: + (key,valaca,valtp)=line.split() + daca[int(key)]=valaca + dtp[int(key)]=valtp + + +imageListACA=[] +imageListTP=[] +imageListACAwt=[] + + +for i in range(fieldnum): + imageListACAwt.append(glob.glob("*"+daca[i]+"*.weight")) + +for i in range(fieldnum): + imageListACA.append(glob.glob("*"+daca[i]+"*.image.pbcor")) + +for i in range(fieldnum): + imageListTP.append(glob.glob("*"+dtp[i]+"*.sd.fits.image")) + + +# preparing a template CMZ cube +## using the mosaic template from globus under /mosaics +importfits(fitsimage = '7m_CS_2-1_m0_mosaic.fits', + imagename = '7m_CS_2-1_m0_mosaic.image') + +imhead(imagename='7m_CS_2-1_m0_mosaic.image', + mode='put', + hdkey='restfreq',hdvalue='97.980953GHz',verbose=False) + +## regrid the mosaic template to Celestial coordiate (it appears that feather is not working well in Galactic coordinate) +imregrid(imagename = '7m_CS_2-1_m0_mosaic.image', + template='J2000', + output='7m_CS_2-1_m0_mosaic.radec.image') + +## regrid (any) 7m line cube to the CMZ template (only X and Y) +imregrid(imagename = 'uid___A001_X1590_X30ac.s38_0.Sgr_A_star_sci.spw24.cube.I.iter1.image.pbcor', + template='7m_CS_2-1_m0_mosaic.radec.image', + output='CS2-1-temlate-whole-band.image', axes=[0, 1]) + +# Imtrans TP map, here needs to swap velocity and stokes axes to match the 7m data header + +for i in range(fieldnum): + imtrans(imagename=imageListTP[i][0], + outfile=imageListTP[i][0]+'.imtrans', + order="0132") + +# regrid TP data to match 7m data +for i in range(fieldnum): + imregrid(imagename=imageListTP[i][0]+'.imtrans', + template=imageListACA[i][0], + output=imageListTP[i][0]+'.imtrans.regrid', + overwrite=False, + axes=[-1]) + +# feather 7m and TP data +for i in range(fieldnum): + feather(imagename=imageListACA[i][0]+'.feather', + highres=imageListACA[i][0], + lowres=imageListTP[i][0]+'.imtrans.regrid', + sdfactor=1.0,effdishdiam=-1.0,lowpassfiltersd=False) + +# smooth feather cubes to the same beam size (Here I check the maximum beam by visual inspection of 7m cube). +# To speed up this process, I only smooth data from channel 0 to 1550 +for i in range(fieldnum): + imsmooth(imagename=imageListACA[i][0]+'.feather', + kernel="gauss", + targetres=True, + major="20arcsec", + minor= "20arcsec", + pa= "0deg", + chans="0~1550", # updated this number depending on your need, or commend out this line if you want to use whole band + outfile=imageListACA[i][0]+'.feather.imsmooth') + +# create a narrower band to speed up mosaic, you don't need this step if you want to use whole band +for i in range(fieldnum): + imsubimage(imagename=imageListACA[i][0]+'.feather.imsmooth', + chans="369~935", + outfile=imageListACA[i][0]+'.feather.imsmooth.imsubimage') + +for i in range(fieldnum): + imsubimage(imagename=imageListACAwt[i][0], + chans="369~935", + outfile=imageListACAwt[i][0]+'.imsubimage') + +## regrid feather cubes to match the template for mosaic + +for i in range(fieldnum): + imregrid(imagename=imageListACA[i][0]+'.feather.imsmooth.imsubimage', + template=templateIM, + output=imageListACA[i][0]+'.feather.imsmooth.imsubimage.regrid', + overwrite=False, + axes=[-1]) + +# deleting masks created during regriding, this step is necessary. +for i in range(fieldnum): + ia.open(imageListACA[i][0]+'.feather.imsmooth.imsubimage.regrid') + ia.maskhandler('delete',['mask0']) + ia.close() + + +for i in range(fieldnum): + imregrid(imagename=imageListACAwt[i][0]+'.imsubimage', + template=templateIM, + output=imageListACAwt[i][0]+'.regrid', + overwrite=False, + axes=[-1]) + +# deleting masks created during regriding, this step is necessary. +for i in range(fieldnum): + ia.open(imageListACAwt[i][0]+'.regrid') + ia.maskhandler('delete',['mask0']) + ia.close() + +# mosacing maps with weighted averaged +# expor = (IM0*IM1+IM2*IM3+...)/(IM1+IM3+...) +# IM0, 2, 4,... are line cubes, IM1, 3, 5, ... are weight cubes +imageList=[] +num = '(' +den = '(' +for i in range(fieldnum): + imageList.append(imageListACA[i][0]+'.feather.imsmooth.imsubimage.regrid') + imageList.append(imageListACAwt[i][0]+'.regrid') + tempNum = 'IM'+str(2*i)+'*IM'+str(2*i+1) + tempDen = 'IM'+str(2*i+1) + if i < len(imageListACA) -1: + num = num + tempNum + '+' + den = den + tempDen + '+' + else: + num = num + tempNum + ')' + den = den + tempDen + ')' +expr = num + '/' + den + +immath(imagename=imageList, + outfile='CS21-feather-Mosaic.image', + expr=expr) + +# convert coordinate from Celestial to Galactic coordinate. +imregrid(imagename='CS21-feather-Mosaic.image', + template='Galactic', + output='CS21-feather-Mosaic.gal.image') + +# trim the blank area to reduce cube size +imsubimage(imagename='CS21-feather-Mosaic.gal.image', + outfile='CS21-feather-Mosaic.gal.imsub.image', + box='0,555,1595,1241') + +# make mom8 map of individual field +for i in range(fieldnum): + immoments(imagename=imageListACA[i][0]+'.feather', + includepix=[20e-2,100], + moments=[8], + outfile=imageListACA[i][0]+'.feather.mom8', + axis='spectral') + +# convert coordinate from Celestial to Galactic coordinate of individual field +for i in range(fieldnum): + imregrid(imagename=imageListACA[i][0]+'.feather.mom8', + template='Galactic', + output=imageListACA[i][0]+'.gal.feather.mom8') + +for i in range(fieldnum): + exportfits(imagename=imageListACA[i][0]+'.gal.feather.mom8', + fitsimage=imageListACA[i][0]+'.gal.feather.mom8.fits') + From f59544daf52032aa38870fc685e4cc98bedc892a Mon Sep 17 00:00:00 2001 From: Pei-Ying Hsieh <95929782+pyhsiehATalma@users.noreply.github.com> Date: Tue, 20 Jun 2023 16:49:38 +0800 Subject: [PATCH 2/2] Update TP-ACA-feather.py changing line name from linename=glob.glob('*.image.pbcor') to linename=glob.glob('*spw33*.image.pbcor') --- aces/imaging/TP-ACA-feather.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/aces/imaging/TP-ACA-feather.py b/aces/imaging/TP-ACA-feather.py index f78eb9c6..28115fd9 100644 --- a/aces/imaging/TP-ACA-feather.py +++ b/aces/imaging/TP-ACA-feather.py @@ -7,7 +7,7 @@ templateIM='CS2-1-temlate-whole-band.image' import glob -linename=glob.glob('*.image.pbcor') +linename=glob.glob('*spw33*.image.pbcor') linewt=glob.glob('*.weight')