Yuzhu Sun

Blog

I will upload some mini-guides (just based on my experience) of software and packages here, as kinds of working notes.

Useful tools

ALFALFA https://vizier.cds.unistra.fr/viz-bin/VizieR-3?-source=J/ApJ/861/49/table2

RA, DEC convert: https://www.astrouw.edu.pl/~jskowron/ra-dec/

plotting: aplpy

PV plotting

1. Find the path you want to extract the pv-plot using Carta

Choose coordinate: image

You can find the 'start' and the 'end' of the line in pixel

Start:(202.84,218.65)

End:(190.81,182,77)

2. Back to python

from spectral_cube import SpectralCube

from pvextractor import extract_pv_slice, Path

cube = SpectralCube.read('/data/yzsun/13A-028/combining/final_cubes_2/make_moments/line_only_AGC226178_BCDD_nat_sub_freq_v2_NHI.fits')

#from the start to the end

path = Path([(202.84,218.65),(190.81,182,77)])

#you can plot the path on one channel, it did work before, but now it encountered: 'Path' object has no attribute 'show_on_axis'. If this can work, you do not even need to find the path in Carta.

ax = pl.subplot(111, projection=cube.wcs.celestial)

#choose one channel to show

ax.imshow(cube[100].value)

path.show_on_axis(ax, spacing=1, color='r')

ax.set_xlabel(f"Right Ascension [{cube.wcs.wcs.radesys}]")

ax.set_ylabel(f"Declination [{cube.wcs.wcs.radesys}]")

3.Extract pv slice

pvdiagram = extract_pv_slice(cube=sub_cube, path=path, spacing=3)

4. Plotting

fig = plt.figure(figsize=(12,12))

ax = pl.subplot(111, projection=wcs.WCS(pvdiagram.header))

im = ax.imshow(pvdiagram.data)

# adjust the shape of the plot if you are not satisfied with it

ax.set_aspect(0.5)

cb = pl.colorbar(mappable=im)

# we could specify the colorbar units like this:

cb.set_label(cube.unit)

# but the 'BUNIT' keyword is not set for these data, so we don't know the unit. We instead manually specify:

#cb.set_label("Brightness Temperature [K]")

#level = (0.0024,0.0048,0.0072)

#ax.contour(pvdiagram.data, levels=level, colors=['white'])

baygaud-PI installation

updated 13 Nov 2023

# create environment

conda create -n baygaud

# activate environment

source activate baygaud

# install python

conda install python=3.10

# download from GitHub

git clone https://github.com/seheon-oh/baygaud-PI.git

# go to the downloaded file

cd baygaud-PI

# installation

pip install .

# using (in the environment)

python3.10 baygaud.py _baygaud_params.xxxx.yaml

#once done, run classify file to classify the warm, hot, cold components

python3.10 baygaud_classify.py _baygaud_params.xxxx.yaml 1

#check the result

python3.10 baygaud_viewer.py _baygaud_params.xxxx.yaml 1

#there will be a window where you can see the gaussian decomposition result for every pixel.

#and the final result (fits files and npy files) can be checked in the merge folder (generated when run baygaud_classify.py).

parameter setting

please check:

https://sites.google.com/view/baygaud

https://github.com/seheon-oh/baygaud-PI/blob/master/README.md

for details

something should pay attention to

the region to gaussian fitting should be identified in 'naxis', cause it will affect the time the code took

the VELOCITY UNIT of the input data cube MUST be in m/s.

s/r should be adjusted based on the quantity of the data, sometimes the default value is too large.

Overlay HI contour (Ha contour) on the Optical image

----------------------------------------------------------

1. Optical cutout making

----------------------------------------------------------

Cut out using astrocut

from astropy.coordinates import SkyCoord

from astrocut import fits_cut

import aplpy

input_files = ['/data/yzsun/NGVS/NGVS+4-2.l.g.Mg004.fits']

center_coord = SkyCoord('191.6804 10.37', unit='deg')

cutout_size = [1000,1000]

cutout_file = fits_cut(input_files, center_coord, cutout_size)

print(cutout_file)

----------------------------------------------------------

Another method to make the cutout

BUT IT CAN NOT BE USED HERE

NGVS42g = fits.open('/data/yzsun/NGVS/NGVS+4-2.l.g.Mg004.fits')

ngvs42g = NGVS42g[0].data

ngvs_wcs = WCS(NGVS42g[0].header)

position = SkyCoord(191.6804*u.deg, 10.37*u.deg,frame="icrs")

cutout = Cutout2D(ngvs42g, position, (1000, 1000), ngvs_wcs)

norm = ImageNormalize(stretch=LogStretch(a=2000), vmin=0, vmax=1000)

plt.figure(figsize=(15, 15))

plt.imshow(cutout.data, origin='lower', cmap='Greys_r',norm=norm)

path_AGC226178_g = 'AGC226178_g.fits'

fits.writeto(path_AGC226178_g, cutout.data, overwrite=True)

----------------------------------------------------------

2. Overlay contour on the optical image

# read the optical image

gc = aplpy.FITSFigure('./cutout_191.676600_10.366100_500-x-500_astrocut.fits')

gc.show_grayscale(invert=True)

# show Ha contour

gc.show_contour('/data/yzsun/Ha/AGC226178_Ha_NET_flux.fits.gz',levels=(3*0.75,5*0.75), colors='red')

#show HI contour

gc.show_contour('/data/yzsun/13A-028/combining/final_cubes/make_moment/1016_226178_mom_nat_4_mom0.fits',levels = (50*7.4e-4,60*7.4e-4,70*7.4e-4,80*7.4e-4,85*7.4e-4,87*7.4e-4), colors='blue')

gc.show_contour('/data/yzsun/13A-028/combining/final_cubes/make_moment/1016_226178_mom_nat_4_mom0.fits',levels = (3*7.4e-4,100), colors='white')

#add scale bar

gc.add_scalebar(length=10*u.arcsecond)

gc.scalebar.set_label('10 arcsec')

gc.scalebar.set_color('white')

plt.title('AGC226178 u-band')