Added convert_brightness_units to branch

resolve/ubik_tools/convert_brightness_units.py

+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)
+