From 62ab2d13561fa989bb8433e7591ccba5575f59cd Mon Sep 17 00:00:00 2001 From: Ashley Barnes Date: Mon, 7 Aug 2023 14:41:52 +0200 Subject: [PATCH] scripts to feather ALMA 12m + GBT inital commit --- aces/joint_deconvolution_cont/crop.py | 106 ++++++++++++++++++ aces/joint_deconvolution_cont/feather.py | 23 ++++ aces/joint_deconvolution_cont/finalise.py | 51 +++++++++ aces/joint_deconvolution_cont/process_alma.py | 28 +++++ .../process_mustang.py | 38 +++++++ 5 files changed, 246 insertions(+) create mode 100644 aces/joint_deconvolution_cont/crop.py create mode 100644 aces/joint_deconvolution_cont/feather.py create mode 100644 aces/joint_deconvolution_cont/finalise.py create mode 100644 aces/joint_deconvolution_cont/process_alma.py create mode 100644 aces/joint_deconvolution_cont/process_mustang.py diff --git a/aces/joint_deconvolution_cont/crop.py b/aces/joint_deconvolution_cont/crop.py new file mode 100644 index 00000000..0e47ae7f --- /dev/null +++ b/aces/joint_deconvolution_cont/crop.py @@ -0,0 +1,106 @@ +from glob import glob +from astropy.io import fits +from reproject import reproject_interp +from reproject.mosaicking import find_optimal_celestial_wcs +import numpy as np +from tqdm.auto import tqdm +from astropy.wcs import WCS +from astropy.coordinates import SkyCoord +from astropy.nddata import Cutout2D +import astropy.units as u +import gc +import os + +def get_croppeddata(hdu, region, square=False): + + hdu_crop = hdu.copy() # Create a copy of the input HDU object + # del hdu # Delete the original HDU object to free up memory + # _ = gc.collect() # Perform garbage collection + + wcs = WCS(hdu_crop) # Create a WCS object from the HDU header + + try: + if square: + radius = max([region['width'], region['height']]) # Calculate the radius for the square cutout + cutout = Cutout2D(hdu_crop.data, region['position'], radius, wcs=wcs) # Create a square cutout + else: + cutout = Cutout2D(hdu_crop.data, region['position'], [region['width'], region['height']], wcs=wcs) # Create a rectangular cutout + except: + print('[INFO] NO cross-over') + return(None) + + hdu_crop.data = cutout.data # Update the data in the cropped HDU object + hdu_crop.header.update(cutout.wcs.to_header()) # Update the header of the cropped HDU object with the cutout's WCS information + + header = hdu.header + header_out = hdu_crop.header + header_keywords = np.array(header.cards)[:,0] + keys = ['BMAJ', 'BMIN', 'BPA', 'BUNIT'] + for key in keys: + if key in header_keywords: + header_out[key] = header[key] + + if 'BUNIT' not in header_keywords: + header_out['BUNIT'] = 'Jy/beam' + + hdu_crop.header = header_out + + del cutout # Delete the cutout to free up memory + _ = gc.collect() # Perform garbage collection + + return hdu_crop # Return the cropped HDU object + +def get_region(l, b, w, h, frame='galactic'): + region = {'position': SkyCoord(l=l *u.deg, b=b *u.deg, frame=frame), + 'width': w *u.deg, + 'height': h *u.deg} + return(region) + +# Get a list of all .fits files in the 'data' directory +files = ['../data/processed/cont_tp.fits', '../data/processed/cont_12m.fits'] + +# Define the galactic coordinates and the dimensions of the desired region in the sky +l = 0.4697169 # galactic longitude +b = -0.0125704 # galactic latitude +width = 1.1808514 +height = 1.1487904 + +# Get the defined region +region = get_region(l, b, height, width) + +# Loop over each file. We're using tqdm to create a progress bar. +# tqdm automatically determines the total number of iterations from the length of the 'files' list. +for file in tqdm(files, desc="Processing files", unit="file"): + + # Print the name of the current file being processed + print('[INFO] infile - %s' % file) + + # Open the FITS file + hdu = fits.open(file)[0] + + # Crop the data in the FITS file to the defined region + hdu_crop = get_croppeddata(hdu, region) + + # Create a WCS (World Coordinate System) object from the cropped HDU + wcs = WCS(hdu_crop) + + # Calculate the middle pixel coordinates of the cropped image + y_mid, x_mid = hdu_crop.shape[0]/2, hdu_crop.shape[1]/2 + + # Convert the middle pixel coordinates to galactic coordinates + l_mid = wcs.pixel_to_world(x_mid, y_mid).l + b_mid = wcs.pixel_to_world(x_mid, y_mid).b + + # Update the header information of the cropped HDU + hdu_crop.header['BUNIT'] = 'Jy/beam' # unit of the data + hdu_crop.header['CRPIX1'] = x_mid # x-coordinate of the reference pixel + hdu_crop.header['CRPIX2'] = y_mid # y-coordinate of the reference pixel + hdu_crop.header['CRVAL1'] = l_mid.value # galactic longitude of the reference pixel + hdu_crop.header['CRVAL2'] = b_mid.value # galactic latitude of the reference pixel + + # If the cropping was successful, save the cropped data to a new FITS file + if hdu_crop is not None: + # Construct the name of the output file by replacing '.fits' with '_cropped.fits' in the input file name + output_file = file.replace('.fits', '_cropped.fits') + # Write the cropped HDU to the output file, overwriting it if it already exists + hdu_crop.writeto(output_file, overwrite=True) diff --git a/aces/joint_deconvolution_cont/feather.py b/aces/joint_deconvolution_cont/feather.py new file mode 100644 index 00000000..0d6da97c --- /dev/null +++ b/aces/joint_deconvolution_cont/feather.py @@ -0,0 +1,23 @@ +import os +from glob import glob +from casatasks import importfits, feather, exportfits + +os.system('rm -rf *.last') +os.system('rm -rf *.log') +os.system('rm -rf *.img') + +inputfile = './../data/processed/cont_tp_cropped.fits' +print(inputfile) +importfits(inputfile, inputfile.replace('fits', 'img'), overwrite=True) + +inputfile = './../data/processed/cont_12m_cropped.fits' +print(inputfile) +importfits(inputfile, inputfile.replace('fits', 'img'), overwrite=True) + +feather(imagename='./../data/feathered/cont_12mtp.img', + highres='./../data/processed/cont_12m_cropped.img', + lowres='./../data/processed/cont_tp_cropped.img') + +exportfits(imagename='./../data/feathered/cont_12mtp.img', + fitsimage='./../data/feathered/cont_12mtp.fits', + overwrite=True) \ No newline at end of file diff --git a/aces/joint_deconvolution_cont/finalise.py b/aces/joint_deconvolution_cont/finalise.py new file mode 100644 index 00000000..fcbac565 --- /dev/null +++ b/aces/joint_deconvolution_cont/finalise.py @@ -0,0 +1,51 @@ +# importing required libraries +from glob import glob +from astropy.io import fits +from reproject import reproject_interp +import numpy as np +from tqdm.auto import tqdm +from astropy.wcs import WCS +from astropy.coordinates import SkyCoord +from astropy.nddata import Cutout2D +import astropy.units as u +import gc +import os + +# defining the path of the FITS file to be opened and processed +file = './../data/feathered/cont_12mtp.fits' + +# opening the FITS file and extracting the first Header Data Unit (HDU) +hdu_12mtp = fits.open(file)[0] + +# removing the singleton dimensions of the data +hdu_12mtp.data = np.squeeze(hdu_12mtp.data) + +# removing some specific header entries +del hdu_12mtp.header['*3*'] +del hdu_12mtp.header['*4*'] + +# writing the modified HDU back to a FITS file, with the specified name +hdu_12mtp.writeto('./../data/feathered/cont_12mtp_final.fits') + +# defining the path of another FITS file to be opened +file = './../data/processed/cont_tp_cropped.fits' + +# opening the second FITS file and extracting the first HDU +hdu_tp = fits.open(file)[0] + +# reprojecting the second HDU to match the coordinate system of the first HDU +# and ignoring the footprint (the second output of the reproject_interp function) +data, _ = reproject_interp(hdu_tp, hdu_12mtp.header) + +# calculating the ratio of the beam sizes (based on the major and minor axes) +# of the two HDUs +beam_ratio = (hdu_12mtp.header['BMAJ']*hdu_12mtp.header['BMIN'])/(hdu_tp.header['BMAJ']*hdu_tp.header['BMIN']) + +# multiplying the reprojected data by the beam ratio +data = data * beam_ratio + +# replacing NaN values in the first HDU's data with corresponding values from the reprojected data +hdu_12mtp.data[np.isnan(hdu_12mtp.data)] = data[np.isnan(hdu_12mtp.data)] + +# writing the final modified HDU back to a new FITS file +hdu_12mtp.writeto('./../data/feathered/cont_12mtp_final_filled.fits') diff --git a/aces/joint_deconvolution_cont/process_alma.py b/aces/joint_deconvolution_cont/process_alma.py new file mode 100644 index 00000000..ac4ae3d8 --- /dev/null +++ b/aces/joint_deconvolution_cont/process_alma.py @@ -0,0 +1,28 @@ +# importing required libraries +from glob import glob +from astropy.io import fits +from reproject import reproject_interp +from reproject.mosaicking import find_optimal_celestial_wcs +import numpy as np +from tqdm.auto import tqdm +from astropy.wcs import WCS +from astropy.coordinates import SkyCoord +from astropy.nddata import Cutout2D +import astropy.units as u +import gc +import os + +# defining the path of the FITS file to be opened and processed +file = '../data/raw/12m_continuum_commonbeam_circular_reimaged_mosaic.fits' + +# opening the FITS file and extracting the first Header Data Unit (HDU) +# HDU is the highest level component of the FITS file structure +hdu = fits.open(file)[0] + +# setting the BUNIT (brightness unit) to 'Jy/beam' +# This indicates that the data is in units of Janskys per beam +hdu.header['BUNIT'] = 'Jy/beam' + +# writing the modified HDU back to a FITS file, with the specified name, +# in the specified location, overwriting the file if it already exists. +hdu.writeto('../data/processed/cont_12m.fits', overwrite=True) diff --git a/aces/joint_deconvolution_cont/process_mustang.py b/aces/joint_deconvolution_cont/process_mustang.py new file mode 100644 index 00000000..94d34701 --- /dev/null +++ b/aces/joint_deconvolution_cont/process_mustang.py @@ -0,0 +1,38 @@ +# importing required libraries +from glob import glob +from astropy.io import fits +from reproject import reproject_interp +from reproject.mosaicking import find_optimal_celestial_wcs +import numpy as np +from tqdm.auto import tqdm +from astropy.wcs import WCS +from astropy.coordinates import SkyCoord +from astropy.nddata import Cutout2D +import astropy.units as u +import gc +import os + +# defining the path of the FITS file to be opened and processed +file = '../data/raw/SgrB2_5pass_1_.0.2_10mJy_10mJy_w_session5_final_smooth4_PlanckCombined_10feb2020.fits' + +# opening the FITS file and extracting the first Header Data Unit (HDU) +# HDU is the highest level component of the FITS file structure +hdu = fits.open(file)[0] + +# modifying the BMAJ and BMIN header values +# BMAJ and BMIN typically represent the major and minor axes of the beam, +# in degrees. Here they are set to 9 arcseconds. +hdu.header['BMAJ'] = 9/3600 +hdu.header['BMIN'] = 9/3600 + +# setting the BPA (Beam Position Angle) header value to 0 +# This angle is usually measured in degrees from north through east +hdu.header['BPA'] = 0 + +# setting the BUNIT (brightness unit) to 'Jy/beam' +# This indicates that the data is in units of Janskys per beam +hdu.header['BUNIT'] = 'Jy/beam' + +# writing the modified HDU back to a FITS file, with the specified name, +# in the specified location, overwriting the file if it already exists. +hdu.writeto('../data/processed/cont_tp.fits', overwrite=True)