forked from rohitk-10/make_cutouts
-
Notifications
You must be signed in to change notification settings - Fork 0
/
cont_to_reg.py
167 lines (119 loc) · 5.77 KB
/
cont_to_reg.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
import numpy as np
from matplotlib import pyplot as plt
from astropy.io import fits
from astropy.wcs import wcs, WCS
from os import remove
import time
import os
from astropy.table import Table
import glob
############################################
def cont_to_ds9reg(lofar_ra, lofar_dec, regfile_name, half_region_size, nlevels, im_data, im_wcs, rms_imdata):
"""
Generate contours using pyplot and then parse them as ds9 polygon regions and write to file
Inputs:
lofar_ra: RA (deg) of the LOFAR source
lofar_dec: DEC (deg) of the LOFAR source
regfile_name: Output filename of the region file
half_region_size: Half the size of the region on which
nlevels : Number of contour levels
Outputs:
ds9 .reg file written to file
"""
# Convert world to pixel coordinates
lof_x, lof_y = im_wcs.all_world2pix(lofar_ra, lofar_dec, 0)
lof_x = int(lof_x)
lof_y = int(lof_y)
# Crop around this region
cropped_region = im_data[lof_y - half_region_size:lof_y + half_region_size, lof_x - half_region_size: lof_x + half_region_size]
# Generate a cropped region around the source in the RMS map and take the median to get the RMS value for contour levels
rms_crop_region = rms_imdata[lof_y - 100:lof_y + 100, lof_x - 100: lof_x + 100]
lofar_rms = np.nanmedian(rms_crop_region)
# plt.imshow(cropped_region, cmap='rainbow', vmin=0.00008, vmax=0.0009)
# Scale factor to generate contour levels from RMS
sf = 3
# Generate levels for this LOFAR source based on the RMS values
cont_levels = np.sqrt(2)**(np.arange(nlevels))
cont_levels = cont_levels * lofar_rms * sf
cont_levels[0] = -1 # Set the first contour level to show "holes"
# Generate the contour colours for each level
cont_colours = ['green' for i in range(len(cont_levels))]
cont_colours[0] = 'red' # Green for the -1 level. The rest are coloured white
# Get the size of the image to get x and y contour
naxis1 = np.size(cropped_region, axis=0)
naxis2 = np.size(cropped_region, axis=1)
# Generate the contours
cont = plt.contour(np.arange(naxis1), np.arange(naxis2), cropped_region, levels=cont_levels, linewidths=2.)
# Get the contour segments
cont_segs = cont.allsegs
# Remove the file if it already exists
try:
remove(regfile_name)
except OSError:
pass
# Generate the region file of the polygons for each level
with open(regfile_name, 'w') as regout:
regout.write("# Region file format: DS9 version 4.1\n")
regout.write('global color=green dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 highlite=1 dash=0 fixed=0 edit=1 move=1 delete=1 include=1 source=1\n')
regout.write("fk5\n")
# Now write in a for loop, the polygons for each level
for ii in range(len(cont_segs)):
# Get the segments of the ii th level
level_segs = cont_segs[ii]
# For each polygon array in a given level
for k in range(len(level_segs)):
# Shift the pixel coordinates back to original positions before writing to file
kth_polygon = level_segs[k]
kth_polygon[:,0] = kth_polygon[:,0] + lof_x - half_region_size
kth_polygon[:,1] = kth_polygon[:,1] + lof_y - half_region_size
# Convert to world coordinates
ra_cont,dec_cont = im_wcs.all_pix2world(kth_polygon[:,0],kth_polygon[:,1],0)
regout.write("polygon(")
# For each row in the polygon array
# Write the first row out of the for loop
regout.write(str(ra_cont[0]) + ',' + str(dec_cont[0]))
for aa in range(1,len(ra_cont)):
regout.write(',' + str(ra_cont[aa]) + ',' + str(dec_cont[aa]))
regout.write(") # color=" + cont_colours[ii] + '\n')
# End of writing the region file
return
# Make a function that reads in a FITS image and returns the wcs and image
def read_fits(fits_filename):
"""
Function to read the FITS file and store the image array and WCS
Inputs:
fits_filename: Name of the FITS filename
Outputs:
wcs_info: WCS of the FITS image
image_array: Array of the FITS image
"""
list_of_hdu = fits.open(fits_filename)
# Get the image array
image_array = fits.getdata(fits_filename)
image_array = image_array[0][0]
# Get the WCS of the image
wcs_info = wcs.WCS(list_of_hdu[0].header, list_of_hdu).celestial
return image_array, wcs_info
######################################################################################
LOFAR_IM_PATH = "../data/image_full_ampphase_di_m.NS_shift.int.facetRestored.blanked-crop.fits"
LOFAR_RMS_PATH = "../data/image_full_ampphase_di_m.NS_shift.blanked.scaled.rms.fits"
# Read in the image and the RMS map and store the data
lofar_image_data, lofar_im_wcs = read_fits(LOFAR_IM_PATH)
hdul_rms = fits.open(LOFAR_RMS_PATH)
rms_image_data = hdul_rms[0].data[0, 0, :, :]
# Read in the catalogue
cata = Table.read("workflow.txt", format='ascii')
sname_all = cata["Source_Name"]
final_fname = glob.glob("../data/final_cross_match_catalogue-v0.3.fits")[-1]
mlfin_srl = Table.read(final_fname, character_as_bytes=False)
t1 = time.time()
BASE_OUT = "contour_regions"
if not os.path.exists(BASE_OUT):
os.makedirs(BASE_OUT)
for ii, sname in enumerate(sname_all):
print(sname)
reg_fnames = "{0}/cont_{1}.reg".format(BASE_OUT, sname)
full_cat_ind = np.isin(mlfin_srl["Source_Name"], sname)
cont_to_ds9reg(mlfin_srl["RA"][full_cat_ind],mlfin_srl["DEC"][full_cat_ind], reg_fnames, 400, 9, lofar_image_data, lofar_im_wcs, rms_image_data)
print(time.time() - t1)
# plt.show()