diff --git a/resolve/ubik_tools/convert_brightness_units.py b/resolve/ubik_tools/convert_brightness_units.py new file mode 100644 index 0000000000000000000000000000000000000000..8397b6afcc95ed863a0ad482462f5afa23e51a03 --- /dev/null +++ b/resolve/ubik_tools/convert_brightness_units.py @@ -0,0 +1,64 @@ +from astropy.io import fits +import numpy as np +import nifty8 as ift +import resolve as rve +import sys + +def _convert_brightness_units(inputfile, outputfile, observations=[], mode='to_per_str'): + fits_file = fits.open(inputfile) + fits_file.info() + fits_data = fits_file[0].data + fits_data = fits_file[0].data + + # keys modified already + used_keys = {'BITPIX','BUNIT','CDELT1','CDELT2','CDELT3','CDELT4','CRPIX1','CRPIX2','CRPIX3','CRPIX4','CRVAL1','CRVAL2','CRVAL3','CRVAL4','CTYPE1','CTYPE2','CTYPE3','CTYPE4','CUNIT1','CUNIT2','CUNIT3','CUNIT4','DATE-MAP','NAXIS','NAXIS1','NAXIS2','NAXIS3','NAXIS4','OBSDEC','OBSRA','ORIGIN','SIMPLE'} + + # compute beam arear + theta_maj = fits_file[0].header['BMAJ'] # deg + theta_min = fits_file[0].header['BMIN'] # deg + + Omega = np.pi*theta_maj*theta_min / (4 * np.log(2)) # deg**2 + Omega = ((np.pi/180)**2) * Omega # sr + + if mode == 'to_per_str': + fits_data /= Omega # convert to Jy / sr + elif mode == 'to_per_beam': + fits_data *= Omega + else: + raise ValueError(f"Invalid mode encountered: {mode}") + + delt = np.abs(fits_file[0].header['CDELT1']*np.pi/180) # smallest angular size + + sky_space = ift.RGSpace([fits_data.shape[2], fits_data.shape[3]], distances=delt) + sky_domain = rve.default_sky_domain(sdom=sky_space) + + data = np.zeros(sky_domain.shape) + data[0,0,0,:,:] = fits_data[0,0,:,:].T + + data = ift.makeField(sky_domain, data) + + header_override = fits_file[0].header + for key in used_keys: + header_override.pop(key, None) + + if mode == 'to_per_beam': + header_override['BUNIT'] = 'Jy/beam' + + rve.ubik_tools.field2fits(data, outputfile, observations=observations, header_override=header_override) + + +if __name__ == "__main__": + if len(sys.argv) < 3 or len(sys.argv) > 4: + print('Please use "python3 convert_brightness_units.py inputfile outputfile [mode]". Exitting.') + sys.exit(0) + else: + if len(sys.argv) == 4: + mode = sys.argv[3] + else: + mode = 'to_per_str' + + inputfile = sys.argv[1] + outputfile = sys.argv[2] + + _convert_brightness_units(inputfile, outputfile, mode=mode) +