Skip to content

Commit

Permalink
fits_helpers
Browse files Browse the repository at this point in the history
jurjen93 committed Nov 12, 2023
1 parent 3826dd0 commit 7d72057
Showing 2 changed files with 313 additions and 0 deletions.
44 changes: 44 additions & 0 deletions fits_helpers/cut_fitsfile.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,44 @@
"""With this script you can cut a region from a fits file"""

import pyregion
from astropy.io import fits
from numpy import nan
from argparse import ArgumentParser


def parse_args():
"""
Command line argument parser
:return: parsed arguments
"""
parser = ArgumentParser(description='Cut fits file based with region file')
parser.add_argument('--fits_input', help='fits input file', required=True, type=str)
parser.add_argument('--fits_output', help='fits output file', required=True, type=str)
parser.add_argument('--region', help='region file', required=True, type=str)
return parser.parse_args()


def main():
""" Main function"""
args = parse_args()

fitsfile = args.fits_input
regionfile = args.region
outputfits = args.fits_output

hdu = fits.open(fitsfile)

header = hdu[0].header

r = pyregion.open(regionfile).as_imagecoord(header=header)
mask = r.get_mask(hdu=hdu[0], shape=(header["NAXIS1"], header["NAXIS2"])).astype(int)
imagedata = hdu[0].data.reshape(header["NAXIS1"], header["NAXIS2"])
imagedata *= mask
imagedata[imagedata == 0] = nan

hdu = fits.PrimaryHDU(header=header, data=imagedata)
hdu.writeto(outputfits, overwrite=True)


if __name__ == '__main__':
main()
269 changes: 269 additions & 0 deletions fits_helpers/make_mosaic.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,269 @@
"""
This script is meant to mosaic together facets.
The algorithm uses feathering for overlapping facets.
For this you will need:
- facet fits files
- region files
Make sure that if you sort the facet fits files and the region files on string name, the fits and region correspond with
each other on index.
For example, it is best to name the fits and region files:
facet_1.fits, facet_2.fits, ... facet_10.fits, ...
region_1.reg, region_2.reg, ... region_10.reg, ...
Also update the coordinates and object name below the script if you image another field than ELAIS-N1
"""

from __future__ import print_function
import sys
from argparse import ArgumentParser
from reproj_test import reproject_interp_chunk_2d
from auxcodes import flatten
from shapely import geometry
import numpy as np
from astropy.wcs import WCS
from astropy.io import fits
import matplotlib.pyplot as plt
from matplotlib.colors import SymLogNorm, PowerNorm
import pyregion
from astropy.visualization.wcsaxes import WCSAxes
import astropy.units as u
from astropy.coordinates import SkyCoord

plt.style.use('ggplot')

def make_header(fitsfile, fullpixsize):
"""
Make header
:param fitsfile: fits file input
:param fullpixsize: full wide-field pixel size
:return: header
"""
hdu = fits.open(fitsfile)
himsize = fullpixsize // 2
header = fits.Header()
header['BITPIX'] = -32
header['NAXIS'] = 2
header['WCSAXES'] = 2
header['NAXIS1'] = 2 * himsize
header['NAXIS2'] = 2 * himsize
header['CTYPE1'] = 'RA---SIN'
header['CTYPE2'] = 'DEC--SIN'
header['CUNIT1'] = 'deg'
header['CUNIT2'] = 'deg'
header['CRPIX1'] = himsize
header['CRPIX2'] = himsize
header['CRVAL1'] = CRVAL1
header['CRVAL2'] = CRVAL2
header['CDELT1'] = -hdu[0].header['CDELT2']
header['CDELT2'] = hdu[0].header['CDELT2']
header['LATPOLE'] = header['CRVAL2']
header['BMAJ'] = hdu[0].header['BMAJ']
header['BMIN'] = hdu[0].header['BMIN']
header['BPA'] = hdu[0].header['BPA']
header['TELESCOPE'] = 'LOFAR'
header['OBSERVER'] = 'LoTSS'
header['BUNIT'] = 'JY/BEAM'
header['BSCALE'] = 1.0
header['BZERO'] = 0
header['BTYPE'] = 'Intensity'
header['OBJECT'] = OBJECT_NAME
return header


def get_polygon_center(regionfile):
"""
get polygon center
:param regionfile: region file
:return: polygon center
"""
regionfile = open(regionfile, 'r')
lines = regionfile.readlines()
regionfile.close()
try:
polygon = lines[4]
except IndexError:
polygon = lines[3]
polyp = [float(p) for p in polygon.replace('polygon(', '').replace(')', '').replace('\n', '').split(',')]
poly_geo = geometry.Polygon(tuple(zip(polyp[0::2], polyp[1::2])))
return SkyCoord(f'{poly_geo.centroid.x}deg', f'{poly_geo.centroid.y}deg', frame='icrs')


def get_array_coordinates(pix_array, wcsheader):
"""
Get coordinates from pixel
:param pix_array: array with pixel coordinates
:param wcsheader: wcs header
:return: array with coordinates from pixel array
"""
pixarray = np.argwhere(pix_array)
return wcsheader.pixel_to_world(pixarray[:, 0], pixarray[:, 1], 0, 0)[0]


def get_distance_weights(center, coord_array):
"""
Get weights based on center polygon to coordinates in array
:param center: center polygon
:param coord_array: coordinates field
:return: weights from center polygon
"""
return 1 / center.separation(coord_array).value

def rms(image_data):
"""
from Cyril Tasse/kMS
:param image_data: image data array
:return: rms (noise measure)
"""
from past.utils import old_div

maskSup = 1e-7
m = image_data[np.abs(image_data)>maskSup]
rmsold = np.std(m)
diff = 1e-1
cut = 3.
med = np.median(m)
for _ in range(10):
ind = np.where(np.abs(m - med) < rmsold*cut)[0]
rms = np.std(m[ind])
if np.abs(old_div((rms-rmsold), rmsold)) < diff: break
rmsold = rms
print(f'Noise : {str(round(rms * 1000, 4))} {u.mJy/u.beam}')
return rms

def make_image(image_data=None, hdu=None, save=None, cmap: str = 'CMRmap', header=None):
"""
Image your data with this method.
image_data -> insert your image_data or plot full image
cmap -> choose your preferred cmap
"""

RMS = rms(image_data)
vmin = RMS
vmax = RMS*25

if hdu is None:
wcs = WCS(header, naxis=2)
else:
wcs = WCS(hdu[0].header, naxis=2)


fig = plt.figure(figsize=(7, 10), dpi=200)
plt.subplot(projection=wcs)
WCSAxes(fig, [0.1, 0.1, 0.8, 0.8], wcs=wcs)

im = plt.imshow(image_data, origin='lower', cmap=cmap)

if cmap=='Blues':
im.set_norm(SymLogNorm(linthresh = vmin * 2, vmin=vmin, vmax = vmin * 15))
else:
im.set_norm(PowerNorm(vmin=0, vmax=vmax, gamma=1 / 2))

plt.xlabel('Right Ascension (J2000)', size=14)
plt.ylabel('Declination (J2000)', size=14)
plt.tick_params(axis='both', which='major', labelsize=12)

plt.grid(False)
plt.grid('off')

if save is not None:
plt.savefig(save, dpi=250, bbox_inches='tight')
else:
plt.show()

def parse_arg():
"""
Command line argument parser
"""

parser = ArgumentParser(description='make wide-field by combining facets')
parser.add_argument('--resolution', help='resolution in arcsecond', required=True, type=float)
parser.add_argument('--fits', type=str, nargs='+', help='facets to merge')
parser.add_argument('--regions', type=str, nargs='+', help='regions corresponding to facets')
return parser.parse_args()

def main():
"""
Main function
"""

args = parse_arg()

resolution = args.resolution
facets = sorted(args.fits)
regions = sorted(args.regions)

if resolution == 0.3:
pixelscale = 0.1 # arcsec
elif resolution == 0.6:
pixelscale = 0.2
elif resolution == 1.2:
pixelscale = 0.4
else:
sys.exit('ERROR: only use resolution 0.3 or 1.2')

fullpixsize = int(2.5 * 3600 / pixelscale)

header_new = make_header(facets[0], fullpixsize)
print(header_new)

xsize = header_new['NAXIS1']
ysize = header_new['NAXIS2']

isum = np.zeros([ysize, xsize], dtype="float32")
weights = np.zeros_like(isum, dtype="float32")
fullmask = np.zeros_like(isum, dtype=bool)

for n, facet in enumerate(facets):
print(facet)

hdu = fits.open(facet)
hduflatten = flatten(hdu)
wcsheader = WCS(hdu[0].header)

imagedata, _ = reproject_interp_chunk_2d(hduflatten, header_new, hdu_in=0, parallel=False)

reg = regions[n]
polycenter = get_polygon_center(reg)
r = pyregion.open(reg).as_imagecoord(header=header_new)
mask = r.get_mask(hdu=hdu[0], shape=(header_new["NAXIS1"], header_new["NAXIS2"])).astype(int)

make_image(mask*imagedata, None, facet+'.png', 'CMRmap', header_new)

fullmask |= ~np.isnan(imagedata)
coordinates = get_array_coordinates(imagedata, wcsheader)
facetweight = get_distance_weights(polycenter, coordinates).reshape(imagedata.shape) * mask
facetweight[~np.isfinite(facetweight)] = 0 # so we can add
imagedata *= facetweight
imagedata[~np.isfinite(imagedata)] = 0 # so we can add
isum += imagedata
weights += facetweight

make_image(isum, None, facet+'full.png', 'CMRmap', header_new)

print('FACET COMPLETED ... ADDING CORRECT WEIGHTS ...')

isum /= weights
isum[isum == np.inf] = np.nan
isum[~fullmask] = np.nan

make_image(isum, None, 'finalfaceted.png', 'CMRmap', header_new)

hdu = fits.PrimaryHDU(header=header_new, data=isum)

hdu.writeto('full-mosaic.fits', overwrite=True)


if __name__ == '__main__':
#ELAIS-N1
CRVAL1=-117.25
CRVAL2=54.95
OBJECT_NAME='ELAIS-N1'

main()

0 comments on commit 7d72057

Please sign in to comment.