Skip to content

Commit 11bc56c

Browse files
committed
Merge branch 'develop'
2 parents a90f233 + 8bb638a commit 11bc56c

File tree

3 files changed

+27
-29
lines changed

3 files changed

+27
-29
lines changed

src/make_images.py

+15-21
Original file line numberDiff line numberDiff line change
@@ -4,9 +4,7 @@
44
from astropy import constants as const
55
from astropy.coordinates import SkyCoord
66
from astropy.io import fits
7-
from astropy.time import Time
87
from astropy import units as u
9-
import astropy.wcs
108
from astropy.wcs import WCS
119
from matplotlib import colors
1210
from matplotlib.patches import Ellipse
@@ -117,8 +115,7 @@ def make_overlay_usr(source, src_basename, cube_params, patch, opt, base_contour
117115
# Plot negative contours
118116
if np.nanmin(hdulist_hi[0].data) < -base_contour and np.isfinite(base_contour):
119117
ax1.contour(hdulist_hi[0].data, cmap='BuPu_r', linewidths=1.2, linestyles='dashed',
120-
levels=-base_contour * 2 ** np.arange(10, -1, -1),
121-
transform=ax1.get_transform(cubew))
118+
levels=-base_contour * 2 ** np.arange(10, -1, -1), transform=ax1.get_transform(cubew))
122119
ax1.text(0.5, 0.05, nhi_labels, ha='center', va='center', transform=ax1.transAxes,
123120
color='white', fontsize=22)
124121
ax1.add_patch(Ellipse((0.92, 0.9), height=patch['height'], width=patch['width'], angle=cube_params['bpa'],
@@ -190,8 +187,6 @@ def make_overlay(source, src_basename, cube_params, patch, opt, base_contour, sp
190187
ax1 = fig.add_subplot(111, projection=owcs)
191188
if survey == 'hst':
192189
plot_labels(source, ax1, cube_params['default_beam'], x_color='w')
193-
# ax1.imshow(opt[0].data, origin='lower', cmap='twilight', norm=LogNorm(vmax=5))
194-
# ax1.imshow(opt[0].data, origin='lower', cmap='Greys', norm=LogNorm(vmin=-0.003, vmax=30))
195190
ax1.imshow(opt[0].data, origin='lower', cmap='Greys',
196191
norm=PowerNorm(gamma=0.4, vmin=np.percentile(opt[0].data, 20),
197192
vmax=np.percentile(opt[0].data, 99.5)))
@@ -295,10 +290,10 @@ def make_mom0(source, src_basename, cube_params, patch, opt_head, base_contour,
295290
ax1.text(0.5, 0.05, nhi_labels, ha='center', va='center', transform=ax1.transAxes, fontsize=22)
296291
ax1.add_patch(Ellipse((0.92, 0.9), height=patch['height'], width=patch['width'], angle=cube_params['bpa'],
297292
transform=ax1.transAxes, facecolor='darkorange', edgecolor='black', linewidth=1))
298-
cb_ax = fig.add_axes([0.91, 0.11, 0.02, 0.76])
293+
cb_ax = fig.add_axes([0.915, 0.11, 0.02, 0.76])
299294
cbar = fig.colorbar(im, cax=cb_ax)
300295
cbar.set_label("Total Intensity [{}]".format(hdulist_hi[0].header['bunit']), fontsize=22)
301-
cbar.ax.tick_params(labelsize=22)
296+
cbar.ax.tick_params(labelsize=22, length=6, width=2)
302297
# Prevent cbar label from getting too wide e.g. when in Jy/beam Hz units
303298
cbar.formatter.set_powerlimits((0, 3))
304299
ot = cbar.ax.yaxis.get_offset_text()
@@ -375,7 +370,6 @@ def make_snr(source, src_basename, cube_params, patch, opt_head, base_contour, s
375370
norm = colors.BoundaryNorm(boundaries, wa_cmap.N, clip=True)
376371
fig = plt.figure(figsize=(8, 8))
377372
ax1 = fig.add_subplot(111, projection=owcs)
378-
# ax1 = fig.add_subplot(111, projection=hiwcs)
379373
plot_labels(source, ax1, cube_params['default_beam'])
380374
ax1.set(facecolor="white") # Doesn't work with the color im
381375
im = ax1.imshow(np.abs(snr), cmap=wa_cmap, origin='lower', norm=norm, transform=ax1.get_transform(cubew))
@@ -384,10 +378,10 @@ def make_snr(source, src_basename, cube_params, patch, opt_head, base_contour, s
384378
ax1.text(0.5, 0.05, nhi_label, ha='center', va='center', transform=ax1.transAxes, fontsize=22)
385379
ax1.add_patch(Ellipse((0.92, 0.9), height=patch['height'], width=patch['width'], angle=cube_params['bpa'],
386380
transform=ax1.transAxes, facecolor='gold', edgecolor='indigo', linewidth=1))
387-
cb_ax = fig.add_axes([0.91, 0.11, 0.02, 0.76])
381+
cb_ax = fig.add_axes([0.915, 0.11, 0.02, 0.76])
388382
cbar = fig.colorbar(im, cax=cb_ax)
389383
cbar.set_label("Pixel SNR", fontsize=22)
390-
cbar.ax.tick_params(labelsize=22)
384+
cbar.ax.tick_params(labelsize=22, length=6, width=2)
391385

392386
ax1.set_xlim(0, opt_head['NAXIS1'])
393387
ax1.set_ylim(0, opt_head['NAXIS2'])
@@ -538,11 +532,11 @@ def make_mom1(source, src_basename, cube_params, patch, opt_head, opt_view, base
538532
levels = [v_sys-3*vunit, v_sys-2*vunit, v_sys-1*vunit, v_sys, v_sys+1*vunit, v_sys+2*vunit, v_sys+3*vunit]
539533
clevels = ['white', 'lightgray', 'dimgrey', 'black', 'dimgrey', 'lightgray', 'white']
540534
if not singlechansource:
541-
cf = ax1.contour(mom1_d, colors=clevels, levels=levels, linewidths=0.6, transform=ax1.get_transform(cubew))
535+
cf = ax1.contour(mom1_d, colors=clevels, levels=levels, linewidths=1.0, transform=ax1.get_transform(cubew))
542536

543537
v_sys_label = "$v_{{center}}$ = {} km/s".format(int(v_sys))
544538
if source['id'] != 0:
545-
v_sys_label = "$v_{{sys}}$ = {} km/s".format(int(v_sys))
539+
v_sys_label = "$cz_{{sys}}$ = {} km/s".format(int(v_sys))
546540

547541
if source['id'] != 0:
548542
# Plot kin_pa from HI center of galaxy; calculate end points of line
@@ -581,12 +575,12 @@ def make_mom1(source, src_basename, cube_params, patch, opt_head, opt_view, base
581575
transform=ax1.transAxes, color='black', fontsize=22, rotation=90)
582576
ax1.add_patch(Ellipse((0.92, 0.9), height=patch['height'], width=patch['width'], angle=cube_params['bpa'],
583577
transform=ax1.transAxes, edgecolor='darkred', linewidth=1))
584-
cb_ax = fig.add_axes([0.91, 0.11, 0.02, 0.76])
578+
cb_ax = fig.add_axes([0.915, 0.11, 0.02, 0.76])
585579
cbar = fig.colorbar(im, cax=cb_ax)
586580
if not singlechansource:
587581
cbar.add_lines(cf)
588582
cbar.set_label(cbar_label, fontsize=22)
589-
cbar.ax.tick_params(labelsize=22)
583+
cbar.ax.tick_params(labelsize=22, length=6, width=2)
590584

591585
ax1.set_xlim(0, opt_head['NAXIS1'])
592586
ax1.set_ylim(0, opt_head['NAXIS2'])
@@ -684,18 +678,18 @@ def make_mom2(source, src_basename, cube_params, patch, opt_head, base_contour,
684678
levels = np.arange(vunit,vel_max,vunit)
685679
# clevels = ['white', 'lightgray', 'dimgrey', 'black', 'dimgrey', 'lightgray', 'white']
686680
if not singlechansource:
687-
cf = ax1.contour(mom2_d, levels=levels, colors=['k', ], linewidths=0.6, transform=ax1.get_transform(cubew))
681+
cf = ax1.contour(mom2_d, levels=levels, colors=['k', ], linewidths=1.0, transform=ax1.get_transform(cubew))
688682
v_disp_label = r"$\Delta \sigma_{{contours}}$ = {} km/s".format(int(vunit))
689683

690684
ax1.text(0.5, 0.05, v_disp_label, ha='center', va='center', transform=ax1.transAxes, color='black', fontsize=22)
691685
ax1.add_patch(Ellipse((0.92, 0.9), height=patch['height'], width=patch['width'], angle=cube_params['bpa'],
692686
transform=ax1.transAxes, facecolor='#4199B5', edgecolor='#D8424D', linewidth=1))
693-
cb_ax = fig.add_axes([0.91, 0.11, 0.02, 0.76])
687+
cb_ax = fig.add_axes([0.915, 0.11, 0.02, 0.76])
694688
cbar = fig.colorbar(im, cax=cb_ax)
695689
if not singlechansource:
696690
cbar.add_lines(cf)
697691
cbar.set_label(cbar_label, fontsize=22)
698-
cbar.ax.tick_params(labelsize=22)
692+
cbar.ax.tick_params(labelsize=22, length=6, width=2)
699693

700694
ax1.set_xlim(0, opt_head['NAXIS1'])
701695
ax1.set_ylim(0, opt_head['NAXIS2'])
@@ -856,7 +850,6 @@ def make_pv(source, src_basename, cube_params, opt_view=6*u.arcmin, spec_line=No
856850
all_colors = np.vstack((colors_noise, colors_galaxy))
857851
pvd_map = colors.LinearSegmentedColormap.from_list('pvd_map', all_colors)
858852
divnorm = colors.TwoSlopeNorm(vmin=-3*pvd_rms, vcenter=+3*pvd_rms, vmax=15*pvd_rms)
859-
# ax1.imshow(pvd, cmap='gray', aspect='auto', vmin=-3*pvd_rms, vmax=+3*pvd_rms)
860853
ax1.imshow(pvd, cmap=pvd_map, aspect='auto', norm=divnorm)
861854

862855
if np.all(np.isnan(pv[0].data)):
@@ -890,7 +883,7 @@ def make_pv(source, src_basename, cube_params, opt_view=6*u.arcmin, spec_line=No
890883
ax1.plot([0.0, 0.0], [freq1, freq2], c='orange', linestyle='--', linewidth=1.0,
891884
transform=ax1.get_transform('world'))
892885
ax1.set_title(source['name'], fontsize=24)
893-
ax1.tick_params(axis='both', which='major', labelsize=22)
886+
ax1.tick_params(axis='both', which='major', labelsize=22, length=6, width=2)
894887
ax1.coords[0].set_format_unit(u.arcmin)
895888
ax1.coords[0].set_major_formatter('x.xx') # Prevent crowding on x-axis
896889
ax1.set_xlabel('Angular Offset [arcmin]', fontsize=22)
@@ -920,7 +913,7 @@ def make_pv(source, src_basename, cube_params, opt_view=6*u.arcmin, spec_line=No
920913
vel2 = (const.c * (source['freq'] - freq2)/source['freq']).to(u.km / u.s).value
921914
ax2.set_ylim(vel1, vel2)
922915
ax2.set_ylabel('Rest Frame Velocity [km/s]', fontsize=22)
923-
ax2.tick_params(labelsize=22)
916+
ax2.tick_params(labelsize=22, length=6, width=2)
924917
else:
925918
print("\tWARNING: Input cube is in velocity units--no correction to source rest frame velocity has been applied!")
926919
if ('v_rad' in source.colnames) or (cube_params['spec_axis'] == 'VRAD'):
@@ -1090,6 +1083,7 @@ def main(source, src_basename, opt_view=6*u.arcmin, suffix='png', sofia=2, beam=
10901083
# For CHILES: plot HI contours on HST image if desired.
10911084
if ('hst' in surveys) | ('HST' in surveys):
10921085
hst_opt_view = np.array([40,]) * u.arcsec
1086+
if opt_view < hst_opt_view: hst_opt_view = opt_view.to(u.arcsec)
10931087
if np.any(Xsize > hst_opt_view.to(u.arcmin).value / 2) | np.any(Ysize > hst_opt_view.to(u.arcmin).value / 2):
10941088
hst_opt_view = (np.max([Xsize, Ysize]) * 2 * 1.05 * u.arcmin).to(u.arcsec)
10951089
hst_opt = get_hst_cosmos(source, opt_view=hst_opt_view)

src/make_spectra.py

+11-7
Original file line numberDiff line numberDiff line change
@@ -143,7 +143,7 @@ def make_specfull(source, src_basename, cube_params, original, spec_line=None, s
143143
equivalencies=line['convention']).value
144144
maskmax = (spec['velo'][spec['chan'] == source['z_max']] * u.m / u.s).to(u.km / u.s,
145145
equivalencies=line['convention']).value
146-
v_sys_label = "$v_{{sys}}$ = {} $W_{{50}}$ = {}".format(int(v_sys), int(w50_vel))
146+
v_sys_label = "$cz_{{sys}}$ = {} $W_{{50}}$ = {}".format(int(v_sys), int(w50_vel))
147147
if original or len(spec) >= long_format:
148148
v_sys_label += " $W_{{20}}$ = {} km/s".format(int(w20_vel))
149149
else:
@@ -203,7 +203,7 @@ def make_specfull(source, src_basename, cube_params, original, spec_line=None, s
203203
else:
204204
ax2_spec.set_xlabel("{} {} Recessional Velocity [km/s]".format(cube_params['spec_sys'].capitalize(),
205205
line['rad_opt']), fontsize=17)
206-
ax2_spec.tick_params(axis='both', which='major', labelsize=16)
206+
ax2_spec.tick_params(axis='both', which='major', labelsize=16, length=5, width=1.8)
207207
ax2_spec.autoscale(False)
208208
if not original or len(spec) < long_format:
209209
ax2_spec.xaxis.set_major_locator(plt.MaxNLocator(7))
@@ -217,7 +217,7 @@ def make_specfull(source, src_basename, cube_params, original, spec_line=None, s
217217
freq2 = freq2.to(u.GHz)
218218
ax2b_spec.set_xlabel('Frequency [GHz]', fontsize=17)
219219
ax2b_spec.set_xlim(freq1.value, freq2.value)
220-
ax2b_spec.tick_params(labelsize=16)
220+
ax2b_spec.tick_params(labelsize=16, length=5, width=1.8)
221221
ax2b_spec.ticklabel_format(style='plain', useOffset=False)
222222
if not original or len(spec) < long_format:
223223
ax2b_spec.xaxis.set_major_locator(plt.MaxNLocator(7))
@@ -233,6 +233,8 @@ def make_specfull(source, src_basename, cube_params, original, spec_line=None, s
233233
ax2_spec.plot([maskmin, maskmin], [0.95*ymin, 0.95*ymax], ':', color='gray')
234234
ax2_spec.plot([maskmax, maskmax], [0.95*ymin, 0.95*ymax], ':', color='gray')
235235

236+
ax2_spec.plot([v_sys, v_sys], np.array([-0.05, 0.05])*(ymax-ymin), color='gray')
237+
236238
# # Condition from Apertif experience that if the RFI is *really* bad, plot based on strength of HI profile
237239
# ax_margin_percent = 0.15
238240
# if (ymax > 5.*galspec_max) | (ymin < np.nanmin([-2.*np.abs(galspec_min), -5.*np.nanstd(spectrumJy).value])):
@@ -297,10 +299,10 @@ def make_spec(source, src_basename, cube_params, spec_line=None, suffix='png'):
297299
names=['chan', 'velo', 'f_sum', 'n_pix'])
298300
optical_velocity = (spec['velo'] * u.m / u.s).to(u.km / u.s).value
299301
if 'snr' in source.colnames:
300-
v_sys_label = "$v_{{sys}}$ = {} $W_{{50}}$ = {} km/s, SNR = {:.1f}".format(int(v_sys), int(w50_vel),
302+
v_sys_label = "$cz_{{sys}}$ = {} $W_{{50}}$ = {} km/s, SNR = {:.1f}".format(int(v_sys), int(w50_vel),
301303
source['snr'])
302304
else:
303-
v_sys_label = "$v_{{sys}}$ = {} $W_{{50}}$ = {} km/s".format(int(v_sys), int(w50_vel))
305+
v_sys_label = "$cz_{{sys}}$ = {} $W_{{50}}$ = {} km/s".format(int(v_sys), int(w50_vel))
304306

305307
except FileNotFoundError:
306308
print("\tNo *_spec.txt file. Perhaps you ran SoFiA without generating moments?")
@@ -335,7 +337,7 @@ def make_spec(source, src_basename, cube_params, spec_line=None, suffix='png'):
335337
else:
336338
ax1_spec.set_xlabel("{} {} Recessional Velocity [km/s]".format(cube_params['spec_sys'].capitalize(),
337339
line['rad_opt']), fontsize=17)
338-
ax1_spec.tick_params(axis='both', which='major', labelsize=16)
340+
ax1_spec.tick_params(axis='both', which='major', labelsize=16, length=5, width=1.8)
339341
ax1_spec.autoscale(False)
340342
ax1_spec.xaxis.set_major_locator(plt.MaxNLocator(7))
341343
if 'freq' in source.colnames:
@@ -349,8 +351,10 @@ def make_spec(source, src_basename, cube_params, spec_line=None, suffix='png'):
349351
ax1b_spec.set_xlabel('Frequency [GHz]', fontsize=17)
350352
ax1b_spec.set_xlim(freq1.value, freq2.value)
351353
ax1b_spec.ticklabel_format(style='plain', useOffset=False)
352-
ax1b_spec.tick_params(labelsize=16)
354+
ax1b_spec.tick_params(labelsize=16, length=5, width=1.8)
353355
ax1b_spec.xaxis.set_major_locator(plt.MaxNLocator(7))
356+
ymin, ymax = ax1_spec.get_ylim()
357+
ax1_spec.plot([v_sys, v_sys], np.array([-0.05, 0.05])*(ymax-ymin), color='gray')
354358
else:
355359
print('\t{} already exists. Will not overwrite.'.format(outfile1))
356360
fig1, ax1_spec, outfile1 = None, None, None

src/modules/functions.py

+1-1
Original file line numberDiff line numberDiff line change
@@ -465,7 +465,7 @@ def plot_labels(source, ax, default_beam, x_color='k'):
465465
ax.scatter(source['pos_x'], source['pos_y'], marker='x', c=x_color, linewidth=0.75,
466466
transform=ax.get_transform('world'))
467467
ax.set_title(source['name'], fontsize=24)
468-
ax.tick_params(axis='both', which='major', labelsize=22)
468+
ax.tick_params(axis='both', which='major', labelsize=22, length=6, width=2)
469469
ax.coords[x_coord].set_axislabel(x_label, fontsize=22)
470470
ax.coords[y_coord].set_axislabel(y_label, fontsize=22)
471471
if default_beam:

0 commit comments

Comments
 (0)