pAff2sEN

1def make_stars(kintable, fovimage, redshift, name, snthresh=100, velocity_thresh=None, project_by_median=True, project_by_redshift=False, cropbox=None, zoom=False, zoombox=None, wcs=True, nancolor='white', colorbar_scalefactor=0.047, vel_vmin=-500, vel_vmax=500, disp_vmin=0, disp_vmax=300, save=True, file_save_directory="./"):
2    
3    table = fits.getdata(kintable)
4    #hdr = fits.getheader(fovimage)
5    
6    x = table['x_cor']
7    y = table['y_cor']
8    
9    # Get the 2D dimensions into which you'll paint this data
10    fovdata = fits.getdata(fovimage)
11    dim = fovdata.shape
12    
13    # Make empty maps of NaNs
14    velmap = np.full((dim[0], dim[1]), np.nan)
15    dispmap = np.full((dim[0], dim[1]), np.nan)
16    
17    # Threshold & Crop
18    if velocity_thresh is not None:
19        mask = (table['vel_fit'] / table['vel_fit_err'] > snthresh) & (table['vel_fit'] < np.nanmedian(table['vel_fit'] + velocity_thresh)) & (table['vel_fit'] > np.nanmedian(table['vel_fit'] - velocity_thresh))
20    else:
21        mask = table['vel_fit'] / table['vel_fit_err'] > snthresh
22    
23    # YES, y,x, in that order. I know it's confusing. 
24    velmap[y[mask], x[mask]] = table['vel_fit'][mask]
25    dispmap[y[mask], x[mask]] = table['disp_fit'][mask]
26    
27    if cropbox is not None:
28        x1, x2, y1, y2 = cropbox 
29
30        # YES, you should be confused by the below line. 
31        # This is *inverted* for the actual mask, but NOT for the zoom. 
32        # So what the user would naturally expect to be x1 is actually y1
33        keep = (x > x1) & (x < x2) & (y > y1) & (y < y2) # YES
34        crop = np.logical_not(keep)
35    
36        velmap[x[crop], y[crop]] = np.nan
37        dispmap[x[crop], y[crop]] = np.nan 
38
39
40    if project_by_redshift is True:
41        if project_by_median is True:
42            project_by_median = False
43            print("project_by_redshift=True which overrides project_by_median=True")
44        velmap = velmap - redshift * const.c.to(u.km / u.s).value
45    if project_by_median is True:
46        redshift = np.nanmedian(velmap) / const.c.to(u.km / u.s).value
47        velmap = velmap - np.nanmedian(velmap)
48    
49
50    # Create the Figures   
51
52    if wcs is True:
53        wcs = WCS(fits.getheader(fovimage, 1))
54
55        velfig = plt.figure(1, figsize=(10,10))
56        velax = velfig.add_subplot(111, projection=wcs)
57        
58        dispfig = plt.figure(2, figsize=(10,10))
59        dispax = dispfig.add_subplot(111, projection=wcs)
60
61        velax.coords[0].set_axislabel('Right Ascension')
62        velax.coords[1].set_axislabel('Declination')
63        
64        dispax.coords[0].set_axislabel('Right Ascension')
65        dispax.coords[1].set_axislabel('Declination')        
66        
67    elif wcs is False:
68        velfig = plt.figure(1, figsize=(10,10))
69        velax = velfig.add_subplot(111)
70        
71        dispfig = plt.figure(2, figsize=(10,10))
72        dispax = dispfig.add_subplot(111)
73        
74        velax.set_xlabel("X")
75        velax.set_ylabel("Y")
76        
77        dispax.set_xlabel("X")
78        dispax.set_ylabel("Y")
79    
80    velax.grid(False)
81    dispax.grid(False)
82    
83    cmap_vel = cm.RdBu_r
84    cmap_vel.set_bad(nancolor) 
85    
86    cmap_disp = cm.plasma
87    cmap_disp.set_bad(nancolor)
88    
89    if project_by_median is True or project_by_redshift is True:
90        velframe = velax.imshow(velmap, origin='lower', vmin=vel_vmin, vmax=vel_vmax, cmap=cmap_vel, interpolation='nearest')
91        velcbar = velfig.colorbar(velframe, fraction=colorbar_scalefactor, pad=0.01)
92        velcbar.set_label(r"Stellar Velocity (km s$^{{-1}}$) relative to z = {}".format(round(redshift, 4)))
93    else:
94        velframe = velax.imshow(velmap, origin='lower', cmap=cmap_vel, interpolation='nearest')
95        velcbar = velfig.colorbar(velframe, fraction=colorbar_scalefactor, pad=0.01)
96        velcbar.set_label(r"Stellar Velocity (km s$^{-1}$)")
97        vmin=None
98        vmax=None
99
100    dispframe = dispax.imshow(dispmap, origin='lower', vmin=disp_vmin, vmax=disp_vmax, cmap=cmap_disp, interpolation='nearest')
101    dispcbar = dispfig.colorbar(dispframe, fraction=colorbar_scalefactor, pad=0.01)
102    dispcbar.set_label(r"Stellar Velocity Dispersion (km s$^{-1}$)")
103        
104    if zoom is True:
105        if zoombox is not None:
106            x1, x2, y1, y2 = zoombox
107            print("Zooming to {}".format(zoombox))
108        elif zoombox is None and cropbox is not None:
109            print("Using cropbox as zoombox")
110        elif zoombox is None and cropbox is None:
111            raise Exception("Zoom is TRUE but you don't have a crop or zoombox. You must specify at least one!")
112        velax.set_xlim(x1, x2)
113        velax.set_ylim(y1, y2)
114
115        dispax.set_xlim(x1, x2)
116        dispax.set_ylim(y1, y2)
117
118    velax.grid(False)
119    dispax.grid(False)
120
121    # Save Everything
122    
123    if save is True:
124        # Check that the file save directory exists. If not, create it. 
125        if not os.path.exists(file_save_directory):
126            os.makedirs(file_save_directory)
127            print("Creating file save directory: {}".format(file_save_directory))
128        else:
129            print("Found file save directory: {}".format(file_save_directory))
130        
131        # Save the PDF figure
132        velfig_pdf_file = "{}_stellar_velocity.pdf".format(name.replace(" ", ""))
133        velfig.savefig(file_save_directory + velfig_pdf_file, dpi=300, bbox_inches='tight')
134        print("Saved Stellar Velocity Figure to {}".format(velfig_pdf_file))
135        
136        dispfig_pdf_file = "{}_stellar_dispersion.pdf".format(name.replace(" ", ""))
137        dispfig.savefig(file_save_directory + dispfig_pdf_file, dpi=300, bbox_inches='tight')
138        print("Saved Stellar Velocity Dispersion Figure to {}".format(dispfig_pdf_file))
139   
140        # Save the FITS file
141        vel_fits_file = "{}_stellar_velocity.fits".format(name.replace(" ", ""))
142        hdr = WCS(fits.getheader(fovimage, 1)).to_header()
143        hdu = fits.PrimaryHDU(velmap, header=hdr)
144        hdulist = fits.HDUList([hdu])
145        hdulist.writeto(file_save_directory + vel_fits_file, overwrite=True, output_verify='silentfix')
146        print("Saved Stellar Velocity Map FITS image to {}".format(vel_fits_file))
147
148        
149        # Save the FITS figure, along with a WCS
150        
151        disp_fits_file = "{}_stellar_dispersion.fits".format(name.replace(" ", ""))
152        
153        disp_fits_file = "{}_stellar_dispersion.fits".format(name.replace(" ", ""))   
154        hdu = fits.PrimaryHDU(dispmap, header=hdr)
155        hdulist = fits.HDUList([hdu])
156        hdulist.writeto(file_save_directory + disp_fits_file, overwrite=True, output_verify='silentfix')
157        print("Saved Stellar Velocity Dispersion Map FITS image to {}".format(disp_fits_file))
158          
159        # Make Kinemetry Table
160        
161        kinemetry_table_filename = "{}_kinemetry_table.dat".format(name.replace(" ", ""))
162
163        bin_number =np.arange(1, len(table["x_cor"][mask]) + 1)
164        kinemetry_table = Table([bin_number, 
165                                table["x_cor"][mask], 
166                                table["y_cor"][mask], 
167                                table["vel_fit"][mask], 
168                                table["vel_fit_err"][mask], 
169                                table["disp_fit"][mask], 
170                                table["disp_fit_err"][mask]], 
171                                names=["#", "XBIN", "YBIN", "VEL", "ER_VEL", "SIG", "ER_SIG"])  
172        print("Saved Kinemetry Table to {}".format(kinemetry_table_filename)) 
173        ascii.write(kinemetry_table, file_save_directory + kinemetry_table_filename, overwrite=True)