######################################################## # Started Logging At: 2023-02-25 19:13:25 ######################################################## ######################################################## # # Started Logging At: 2023-02-25 19:13:25 ######################################################## ######################################################## # Started Logging At: 2023-02-25 19:23:15 ######################################################## ######################################################## # # Started Logging At: 2023-02-25 19:23:19 ######################################################## get_ipython().run_line_magic('run', 'ACES_Animations.py') get_ipython().run_line_magic('run', 'ACES_Animation.py') get_ipython().run_line_magic('pwd', '') #[Out]# '/orange/adamginsburg/cmz/DataSetVisualizations' get_ipython().run_line_magic('matplotlib', 'inline') import pylab as pl import numpy as np pl.rcParams['image.interpolation'] = 'nearest' import os from astropy.utils.data import download_file from astropy.io import fits from spectral_cube import SpectralCube from astropy.visualization import simple_norm import shutil import reproject from reproject import reproject_from_healpix, reproject_to_healpix from astropy.wcs import WCS import matplotlib.pyplot as plt from astropy.visualization.wcsaxes.frame import EllipticalFrame import astropy.visualization.wcsaxes from matplotlib.colors import rgb_to_hsv, hsv_to_rgb from astropy import coordinates from astropy.coordinates import SkyCoord from astropy import units as u, constants from astropy.convolution import convolve_fft, Gaussian2DKernel, Gaussian1DKernel import healpy import PIL import pyavm import regions def h_rot(rgb, rot): hsv = rgb_to_hsv(rgb) hsv[:,:,0] += rot # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) return rgb_scaled def make_rgb(rgb, basename, hsv_rotation=0, layernames='RGB', do_layers=True, header=target_header, axlims=None, frame_class=EllipticalFrame): plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(header), frame_class=frame_class) #ax.coords.grid(color='white') try: ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) except Exception: ax.coords['ra'].set_ticklabel(visible=False) ax.coords['dec'].set_ticklabel(visible=False) ax.coords['ra'].set_ticks_visible(False) ax.coords['dec'].set_ticks_visible(False) if do_layers: for ind, color in zip((0, 1, 2), (layernames)): rgbc = np.zeros_like(rgb) rgbc[:, :, ind] = rgb[:, :, ind] if hsv_rotation != 0: rgbc = h_rot(rgbc, hsv_rotation) ax.imshow(rgbc) pl.savefig(f'{basename}_{color}.png', bbox_inches='tight', transparent=True) if hsv_rotation != 0: rgb = h_rot(rgb, hsv_rotation) ax.imshow(rgb) if axlims is not None: ax.axis(axlims) pl.savefig(f'{basename}_RGB.png', bbox_inches='tight', transparent=True) return ax get_ipython().run_line_magic('run', 'ACES_Animation.py') ######################################################## # Started Logging At: 2023-02-25 19:25:02 ######################################################## ######################################################## # # Started Logging At: 2023-02-25 19:25:03 ######################################################## get_ipython().run_line_magic('pwd', '') #[Out]# '/orange/adamginsburg/cmz/DataSetVisualizations' get_ipython().run_line_magic('matplotlib', 'inline') import pylab as pl import numpy as np pl.rcParams['image.interpolation'] = 'nearest' import os from astropy.utils.data import download_file from astropy.io import fits from spectral_cube import SpectralCube from astropy.visualization import simple_norm import shutil import reproject from reproject import reproject_from_healpix, reproject_to_healpix from astropy.wcs import WCS import matplotlib.pyplot as plt from astropy.visualization.wcsaxes.frame import EllipticalFrame import astropy.visualization.wcsaxes from matplotlib.colors import rgb_to_hsv, hsv_to_rgb from astropy import coordinates from astropy.coordinates import SkyCoord from astropy import units as u, constants from astropy.convolution import convolve_fft, Gaussian2DKernel, Gaussian1DKernel import healpy import PIL import pyavm import regions target_header_medres = fits.Header.fromstring(""" NAXIS = 2 NAXIS1 = 960 NAXIS2 = 480 CTYPE1 = 'GLON-MOL' CRPIX1 = 480.5 CRVAL1 = 0.0 CDELT1 = -0.3375 CUNIT1 = 'deg ' CTYPE2 = 'GLAT-MOL' CRPIX2 = 240.5 CRVAL2 = 0.0 CDELT2 = 0.3375 CUNIT2 = 'deg ' COORDSYS= 'icrs ' """, sep='\n') target_header = fits.Header.fromstring(""" NAXIS = 2 NAXIS1 = 1920 NAXIS2 = 960 CTYPE1 = 'GLON-MOL' CRPIX1 = 960.5 CRVAL1 = 0.0 CDELT1 = -0.16875 CUNIT1 = 'deg ' CTYPE2 = 'GLAT-MOL' CRPIX2 = 480.5 CRVAL2 = 0.0 CDELT2 = 0.16875 CUNIT2 = 'deg ' COORDSYS= 'icrs ' """, sep='\n') def h_rot(rgb, rot): hsv = rgb_to_hsv(rgb) hsv[:,:,0] += rot # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) return rgb_scaled def make_rgb(rgb, basename, hsv_rotation=0, layernames='RGB', do_layers=True, header=target_header, axlims=None, frame_class=EllipticalFrame): plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(header), frame_class=frame_class) #ax.coords.grid(color='white') try: ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) except Exception: ax.coords['ra'].set_ticklabel(visible=False) ax.coords['dec'].set_ticklabel(visible=False) ax.coords['ra'].set_ticks_visible(False) ax.coords['dec'].set_ticks_visible(False) if do_layers: for ind, color in zip((0, 1, 2), (layernames)): rgbc = np.zeros_like(rgb) rgbc[:, :, ind] = rgb[:, :, ind] if hsv_rotation != 0: rgbc = h_rot(rgbc, hsv_rotation) ax.imshow(rgbc) pl.savefig(f'{basename}_{color}.png', bbox_inches='tight', transparent=True) if hsv_rotation != 0: rgb = h_rot(rgb, hsv_rotation) ax.imshow(rgb) if axlims is not None: ax.axis(axlims) pl.savefig(f'{basename}_RGB.png', bbox_inches='tight', transparent=True) return ax target_header_orion = fits.Header.fromstring(""" NAXIS = 2 NAXIS1 = 1920 NAXIS2 = 960 CTYPE1 = 'RA---MOL' CRPIX1 = 960.5 CRVAL1 = 83.82 CDELT1 = -0.16875 CUNIT1 = 'deg ' CTYPE2 = 'DEC--MOL' CRPIX2 = 480.5 CRVAL2 = -5.39 CDELT2 = 0.16875 CUNIT2 = 'deg ' COORDSYS= 'icrs ' """, sep='\n') def fix_nans(img): kernel = Gaussian2DKernel(2) sm = convolve_fft(img, kernel) img[np.isnan(img)] = sm[np.isnan(img)] return img def fix_rosat_nans(tb): #kernel = Gaussian1DKernel(3) kernel = Gaussian2DKernel(1) # convolving in 2D convolves adjacent scans # tb is a 2D array, but we want to convolve in 1D #sm = convolve_fft(tb, kernel.array[None, :]) sm = convolve_fft(tb, kernel) tb[np.isnan(tb)] = sm[np.isnan(tb)] return tb get_ipython().run_line_magic('pwd', '') #[Out]# '/orange/adamginsburg/cmz/DataSetVisualizations' fn_rosat5 = "RASS_SXRB_R5.fits" if not os.path.exists(fn_rosat5): rosat5 = download_file(f"https://www.jb.man.ac.uk/research/cosmos/rosat/{fn_rosat5}") shutil.move(rosat5, fn_rosat5) rosat5 = fits.open(fn_rosat5) data_rosat5 = rosat5[1].data['COUNT_RATE'] data_rosat5 = fix_rosat_nans(np.where(data_rosat5==0, np.nan, data_rosat5)) pl.hist(data_rosat5.ravel(), bins=np.logspace(-14, 5), log=True) pl.loglog(); #img_rosat5, footprint = reproject_from_healpix(fn_rosat5, target_header) data_rosat5 = rosat5[1].data['COUNT_RATE'] data_rosat5 = fix_rosat_nans(np.where(data_rosat5<=1e-8, np.nan, data_rosat5)) img_rosat5, footprint = reproject_from_healpix((data_rosat5.ravel(), rosat5[1].header["COORDSYS"]), target_header, nested=rosat5[1].header["ORDERING"].lower() == "nested" ) img_rosat5 = fix_nans(np.where(img_rosat5<=1e-9, np.nan, img_rosat5)) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header), frame_class=EllipticalFrame, ) ax.imshow(img_rosat5, norm=simple_norm(img_rosat5, min_percent=1, max_percent=99.99, stretch='log')) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) fn_halpha = "Halpha_map.fits" if not os.path.exists(fn_halpha): halpha = download_file(f"https://lambda.gsfc.nasa.gov/data/foregrounds/fink_halpha/{fn_halpha}") shutil.move(halpha, fn_halpha) halpha = fits.open(fn_halpha) # img_halpha = halpha[0].data # header = halpha[0].header img_halpha, footprint = reproject.reproject_interp(fn_halpha, target_header) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header), frame_class=EllipticalFrame, ) ax.imshow(img_halpha, norm=simple_norm(img_halpha, min_percent=1, max_percent=99.99, stretch='log')) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) fn_akari65 = "akari_mollweide_65_1_4096.fits" if not os.path.exists(fn_akari65): akari65 = download_file(f"https://lambda.gsfc.nasa.gov/data/foregrounds/akari/images/{fn_akari65}") shutil.move(akari65, fn_akari65) akari65 = fits.open(fn_akari65) #img_akari65, footprint = reproject_from_healpix(fn_akari65, target_header) img_akari65 = fix_nans(np.where(akari65[1].data == 0, np.nan, akari65[1].data)) header = akari65[1].header plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(header), #frame_class=EllipticalFrame, ) ax.imshow(img_akari65, norm=simple_norm(img_akari65, min_percent=1, max_percent=99.99, stretch='log')) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) fn_akari90 = "akari_mollweide_WideS_1_4096.fits" if not os.path.exists(fn_akari90): akari90 = download_file(f"https://lambda.gsfc.nasa.gov/data/foregrounds/akari/images/{fn_akari90}") shutil.move(akari90, fn_akari90) akari90 = fits.open(fn_akari90) #img_akari90, footprint = reproject_from_healpix(fn_akari90, target_header) img_akari90 = fix_nans(np.where(akari90[1].data == 0, np.nan, akari90[1].data)) header = akari90[1].header plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(header), #frame_class=EllipticalFrame, ) ax.imshow(img_akari90, norm=simple_norm(img_akari90, min_percent=1, max_percent=99.99, stretch='log')) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) fn_akari140 = "akari_mollweide_WideL_1_4096.fits" if not os.path.exists(fn_akari140): akari140 = download_file(f"https://lambda.gsfc.nasa.gov/data/foregrounds/akari/images/{fn_akari140}") shutil.move(akari140, fn_akari140) akari140 = fits.open(fn_akari140) #img_akari140, footprint = reproject_from_healpix(fn_akari140, target_header) img_akari140 = fix_nans(np.where(akari140[1].data == 0, np.nan, akari140[1].data)) header = akari140[1].header plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(header), #frame_class=EllipticalFrame, ) ax.imshow(img_akari140, norm=simple_norm(img_akari140, min_percent=1, max_percent=99.99, stretch='log')) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) get_ipython().run_line_magic('run', 'ACES_Animation.py') fn_akari160 = "akari_mollweide_160_1_4096.fits" if not os.path.exists(fn_akari160): akari160 = download_file(f"https://lambda.gsfc.nasa.gov/data/foregrounds/akari/images/{fn_akari160}") shutil.move(akari160, fn_akari160) akari160 = fits.open(fn_akari160) #img_akari160, footprint = reproject_from_healpix(fn_akari160, target_header) img_akari160 = fix_nans(np.where(akari160[1].data == 0, np.nan, akari160[1].data)) header = akari160[1].header plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(header), #frame_class=EllipticalFrame, ) ax.imshow(img_akari160, norm=simple_norm(img_akari160, min_percent=1, max_percent=99.99, stretch='log')) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) #fn_chipass = "lambda_chipass_HPX_r10.fits" fn_chipass = "lambda_chipass_mollweide.fits" if not os.path.exists(fn_chipass): chipass = download_file(f"https://lambda.gsfc.nasa.gov/data/foregrounds/chipass/{fn_chipass}") shutil.move(chipass, fn_chipass) chipass = fits.open(fn_chipass) #img_chipass, footprint = reproject_from_healpix(fn_chipass, target_header) img_chipass = chipass[1].data img_chipass[img_chipass <= 0] = chipass[2].data[img_chipass <= 0] img_chipass[img_chipass <= 0] = np.nan header = chipass[1].header plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(header), #frame_class=EllipticalFrame, ) ax.imshow(img_chipass, norm=simple_norm(img_chipass, min_percent=1, max_percent=99.99, stretch='asinh')) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) fn_sve1420 = "STOCKERT+VILLA-ELISA_1420MHz_1_256.fits" if not os.path.exists(fn_sve1420): sve1420 = download_file(f"https://lambda.gsfc.nasa.gov/data/foregrounds/reich_reich/{fn_sve1420}") shutil.move(sve1420, fn_sve1420) sve1420 = fits.open(fn_sve1420) img_sve1420, footprint = reproject_from_healpix(fn_sve1420, target_header) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header), frame_class=EllipticalFrame) ax.imshow(img_sve1420, norm=simple_norm(img_sve1420, min_percent=1, max_percent=99.99, stretch='log')) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) fn_lfi30 = 'LFI_SkyMap_030_1024_R2.01_full.fits' if not os.path.exists(fn_lfi30): lfi30 = download_file(f"https://irsa.ipac.caltech.edu/data/Planck/release_2/all-sky-maps/maps/{fn_lfi30}") shutil.move(lfi30, fn_lfi30) lfi30 = fits.open(fn_lfi30) # don't use this any more m30 = healpy.read_map(lfi30) healpy.mollview(m30, norm=simple_norm(m30, stretch='log'), cbar=False, title='') img_lfi30, footprint = reproject_from_healpix(fn_lfi30, target_header) badinds = np.where(np.isnan(img_lfi30)) img_lfi30 = fix_nans(img_lfi30) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header), frame_class=EllipticalFrame) ax.imshow(img_lfi30, norm=simple_norm(img_lfi30, min_percent=1, max_percent=99.99, stretch='log')) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) fn_hfi143 = "HFI_SkyMap_143_2048_R2.02_full.fits" if not os.path.exists(fn_hfi143): hfi143 = download_file(f"https://irsa.ipac.caltech.edu/data/Planck/release_2/all-sky-maps/maps/{fn_hfi143}") shutil.move(hfi143, fn_hfi143) hfi143 = fits.open(fn_hfi143) img_hfi143, footprint = reproject_from_healpix(fn_hfi143, target_header) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header), frame_class=EllipticalFrame) ax.imshow(img_hfi143, norm=simple_norm(img_hfi143, min_percent=1, max_percent=99.99, stretch='log')) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) fn_hfi857 = "HFI_SkyMap_857_2048_R2.02_full.fits" if not os.path.exists(fn_hfi857): hfi857 = download_file(f"https://irsa.ipac.caltech.edu/data/Planck/release_2/all-sky-maps/maps/{fn_hfi857}") shutil.move(hfi857, fn_hfi857) hfi857 = fits.open(fn_hfi857) img_hfi857, footprint = reproject_from_healpix(fn_hfi857, target_header) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header), frame_class=EllipticalFrame) ax.imshow(img_hfi857, norm=simple_norm(img_hfi857, min_percent=1, max_percent=99.99, stretch='log')) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) rgb = np.array([simple_norm(img_lfi30, min_percent=0.01, max_percent=99.90, log_a=5e1, stretch='log')(img_lfi30), simple_norm(img_hfi143, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_hfi143), simple_norm(img_hfi857, min_percent=0.01, max_percent=99.90, log_a=5e1, stretch='log')(img_hfi857)]).T.swapaxes(0,1) make_rgb(rgb, 'Planck_RGB_30-143-857') #[Out]# fn_synchrotron = "COM_CompMap_Synchrotron-commander_0256_R2.00.fits" if not os.path.exists(fn_synchrotron): planck_synchrotron = download_file(f"https://irsa.ipac.caltech.edu/data/Planck/release_2/all-sky-maps/maps/component-maps/foregrounds/{fn_synchrotron}") shutil.move(planck_synchrotron, fn_synchrotron) planck_synchrotron = fits.open(fn_synchrotron) img_synchrotron, footprint = reproject_from_healpix(fn_synchrotron, target_header) img_synchrotron = fix_nans(np.where(img_synchrotron<=0, np.nan, img_synchrotron)) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header), frame_class=EllipticalFrame) ax.imshow(img_synchrotron, norm=simple_norm(img_synchrotron, min_percent=1, max_percent=99.99, stretch='log')) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) fn_freefree = "COM_CompMap_freefree-commander_0256_R2.00.fits" if not os.path.exists(fn_freefree): planck_freefree = download_file(f"https://irsa.ipac.caltech.edu/data/Planck/release_2/all-sky-maps/maps/component-maps/foregrounds/{fn_freefree}") shutil.move(planck_freefree, fn_freefree) planck_freefree = fits.open(fn_freefree) img_freefree, footprint = reproject_from_healpix(fn_freefree, target_header) img_freefree = fix_nans(img_freefree) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header), frame_class=EllipticalFrame) ax.imshow(img_freefree, norm=simple_norm(img_freefree, min_percent=1, max_percent=99.99, stretch='log')) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) fn_CO21 = "COM_CompMap_CO21-commander_2048_R2.00.fits" if not os.path.exists(fn_CO21): planck_CO21 = download_file(f"https://irsa.ipac.caltech.edu/data/Planck/release_2/all-sky-maps/maps/component-maps/foregrounds/{fn_CO21}") shutil.move(planck_CO21, fn_CO21) planck_CO21 = fits.open(fn_CO21) img_CO21, footprint = reproject_from_healpix(fn_CO21, target_header) img_CO21 = fix_nans(img_CO21) img_CO21_orion, footprint = reproject_from_healpix(fn_CO21, target_header_orion) img_CO21_orion = fix_nans(img_CO21_orion) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header), frame_class=EllipticalFrame) ax.imshow(img_CO21, norm=simple_norm(img_CO21, min_percent=1, max_percent=99.99, stretch='log')) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) fn_ThermalDust = "COM_CompMap_ThermalDust-commander_2048_R2.00.fits" if not os.path.exists(fn_ThermalDust): planck_ThermalDust = download_file(f"https://irsa.ipac.caltech.edu/data/Planck/release_2/all-sky-maps/maps/component-maps/foregrounds/{fn_ThermalDust}") shutil.move(planck_ThermalDust, fn_ThermalDust) planck_ThermalDust = fits.open(fn_ThermalDust) img_ThermalDust, footprint = reproject_from_healpix(fn_ThermalDust, target_header) img_ThermalDust_orion, footprint = reproject_from_healpix(fn_ThermalDust, target_header_orion) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header), frame_class=EllipticalFrame) ax.imshow(img_ThermalDust, norm=simple_norm(img_ThermalDust, min_percent=1, max_percent=99.99, stretch='log')) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) fn_LAB_HI = 'LAB_fullvel.fits' if not os.path.exists(fn_LAB_HI): LAB_HI = download_file(f"https://lambda.gsfc.nasa.gov/data/foregrounds/HI/{fn_LAB_HI}") shutil.move(LAB_HI, fn_LAB_HI) LAB_HI = fits.open(fn_LAB_HI) fn_HI4PI = 'NHI_HPX.fits' if not os.path.exists(fn_HI4PI): HI4PI = download_file(f"https://lambda.gsfc.nasa.gov/data/foregrounds/HI4PI/{fn_HI4PI}") shutil.move(HI4PI, fn_HI4PI) HI4PI = fits.open(fn_HI4PI) img_HI4PI, footprint = reproject_from_healpix(fn_HI4PI, target_header, field=5) img_HI4PI_orion, footprint = reproject_from_healpix(fn_HI4PI, target_header_orion, field=5) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header), frame_class=EllipticalFrame) ax.imshow(img_HI4PI, norm=simple_norm(img_HI4PI, min_percent=1, max_percent=99.99, stretch='log')) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) img_HI4PI_orion, footprint = reproject_from_healpix(fn_HI4PI, target_header_orion, field=5) plt.figure(figsize=(10,5), dpi=200, frameon=False) ax = plt.subplot(1,1,1, projection=WCS(target_header_orion), frame_class=EllipticalFrame) ax.axis('off') ax.imshow(img_HI4PI_orion, norm=simple_norm(img_HI4PI_orion, min_percent=1, max_percent=99.99, stretch='log')) #ax.coords.grid(color='white') ax.coords['ra'].set_ticklabel(visible=False) ax.coords['dec'].set_ticklabel(visible=False) ax.coords['ra'].set_ticks_visible(False) ax.coords['dec'].set_ticks_visible(False) rgb = np.array([simple_norm(img_HI4PI, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PI), simple_norm(img_CO21, min_percent=0.01, max_percent=99.90, log_a=2e1, stretch='log')(img_CO21), simple_norm(img_ThermalDust, min_percent=0.01, max_percent=99.90, log_a=5e2, stretch='log')(img_ThermalDust)]).T.swapaxes(0,1) make_rgb(rgb, 'HI-CO-Dust') #[Out]# rgb = np.array([simple_norm(img_HI4PI, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PI), simple_norm(img_CO21, min_percent=0.01, max_percent=99.90, log_a=2e1, stretch='log')(img_CO21), simple_norm(img_ThermalDust, min_percent=0.01, max_percent=99.90, log_a=5e2, stretch='log')(img_ThermalDust)]).T.swapaxes(0,1) make_rgb(rgb, 'HI-CO-Dust_m0.25', -0.25, layernames=('HI', 'CO', 'Dust')) #[Out]# # different colors, overwrites previous rgb = np.array([simple_norm(img_HI4PI, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PI), simple_norm(img_CO21, min_percent=0.01, max_percent=99.90, log_a=2e1, stretch='log')(img_CO21), simple_norm(img_ThermalDust, min_percent=0.01, max_percent=99.90, log_a=5e2, stretch='log')(img_ThermalDust)]).T.swapaxes(0,1) make_rgb(rgb, 'HI-CO-Dust_m0.35', -0.35, layernames=('HI', 'CO', 'Dust')) #[Out]# # zoom on lower-left corner rgb = np.array([simple_norm(img_HI4PI, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PI), simple_norm(img_CO21, min_percent=0.01, max_percent=99.90, log_a=2e1, stretch='log')(img_CO21), simple_norm(img_ThermalDust, min_percent=0.01, max_percent=99.90, log_a=5e2, stretch='log')(img_ThermalDust)]).T.swapaxes(0,1) make_rgb(rgb, 'HI-CO-Dust_m0.35_zoom', -0.35, layernames=('HI', 'CO', 'Dust'), axlims=(240, 720, 120, 360), frame_class=None) get_ipython().run_line_magic('run', 'ACES_Animation.py') #[Out]# # CMZoom rgb = np.array([simple_norm(img_HI4PI, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PI), simple_norm(img_CO21, min_percent=0.01, max_percent=99.90, log_a=2e1, stretch='log')(img_CO21), simple_norm(img_ThermalDust, min_percent=0.01, max_percent=99.90, log_a=5e2, stretch='log')(img_ThermalDust)]).T.swapaxes(0,1) make_rgb(rgb, 'HI-CO-Dust_m0.35_zoomCMZ', -0.35, layernames=('HI', 'CO', 'Dust'), axlims=(240*2, 720*2, 120*2, 360*2), frame_class=None) #[Out]# img_rosat5.shape, img_halpha.shape, img_synchrotron.shape #[Out]# ((960, 1920), (960, 1920), (960, 1920)) rgb = np.array([simple_norm(img_rosat5, min_percent=0.01, max_percent=99.99, log_a=5e1, stretch='log')(img_rosat5), simple_norm(img_halpha, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_halpha), simple_norm(img_synchrotron, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_synchrotron)]).T.swapaxes(0,1) make_rgb(rgb, 'XRay_Halpha_Synchro', layernames=('X-Ray', 'Halpha', 'Synchrotron')) #[Out]# rgb = np.array([simple_norm(img_rosat5, min_percent=0.01, max_percent=99.99, log_a=5e1, stretch='log')(img_rosat5), simple_norm(img_synchrotron, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_synchrotron), simple_norm(img_HI4PI, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PI),] ).T.swapaxes(0,1) make_rgb(rgb, 'XRay_Synchro_HI4PI', layernames=('X-Ray', 'Synchrotron', 'HI4PI', ), hsv_rotation=-0.25) #[Out]# rgb = np.array([simple_norm(img_rosat5, min_percent=0.01, max_percent=99.99, log_a=5e1, stretch='log')(img_rosat5), simple_norm(img_halpha, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_halpha), simple_norm(img_ThermalDust, min_percent=0.01, max_percent=99.99, log_a=5e1, stretch='log')(img_ThermalDust)]).T.swapaxes(0,1) make_rgb(rgb, 'XRay_Halpha_ThermalDust', layernames=('X-Ray', 'Halpha', 'ThermalDust')) #[Out]# rgb = np.array([simple_norm(img_rosat5, min_percent=0.01, max_percent=99.99, log_a=5e1, stretch='log')(img_rosat5), simple_norm(img_halpha, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_halpha), simple_norm(img_CO21, min_percent=0.01, max_percent=99.99, log_a=5e1, stretch='log')(img_CO21)]).T.swapaxes(0,1) make_rgb(rgb, 'XRay_Halpha_CO21', layernames=('X-Ray', 'Halpha', 'CO21')) #[Out]# rgb = np.array([simple_norm(img_rosat5, min_percent=0.01, max_percent=99.99, log_a=1e3, stretch='log')(img_rosat5), simple_norm(img_halpha, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_halpha), simple_norm(img_CO21, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21)]).T.swapaxes(0,1) make_rgb(rgb, 'XRay_Halpha_CO21', layernames=('X-Ray', 'Halpha', 'CO21'), hsv_rotation=-0.25) #[Out]# xray_synchro = (simple_norm(img_rosat5, min_percent=0.001, max_percent=99.999, stretch='linear')(img_rosat5) + simple_norm(img_synchrotron, min_percent=0.001, max_percent=99.999, stretch='linear')(img_synchrotron)) rgb = np.array([simple_norm(xray_synchro, min_percent=0.01, max_percent=99.99, log_a=1e2, stretch='log')(xray_synchro), simple_norm(img_halpha, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_halpha), simple_norm(img_CO21, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21)]).T.swapaxes(0,1) make_rgb(rgb, 'XRaySynchro_Halpha_CO21', layernames=('X-RaySynchro', 'Halpha', 'CO21'), hsv_rotation=-0.25) #[Out]# rgb = np.array([simple_norm(img_rosat5, min_percent=0.01, max_percent=99.99, log_a=2e2, stretch='log')(img_rosat5), simple_norm(img_halpha, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_halpha), simple_norm(img_HI4PI, min_percent=0.01, max_percent=99.99, log_a=2e2, stretch='log')(img_HI4PI)]).T.swapaxes(0,1) make_rgb(rgb, 'XRay_Halpha_HI4PI', layernames=('X-Ray', 'Halpha', 'HI4PI')) #[Out]# rgb = np.array([simple_norm(img_rosat5, min_percent=0.01, max_percent=99.99, log_a=1e2, stretch='log')(img_rosat5), simple_norm(img_halpha, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_halpha), simple_norm(img_freefree, min_percent=0.01, max_percent=99.99, log_a=1e3, stretch='log')(img_freefree)]).T.swapaxes(0,1) make_rgb(rgb, 'XRay_Halpha_freefree', layernames=('X-Ray', 'Halpha', 'freefree')) #[Out]# rgb = np.array([simple_norm(img_rosat5, min_percent=0.01, max_percent=99.99, log_a=5e2, stretch='log')(img_rosat5), simple_norm(img_halpha, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_halpha), simple_norm(img_freefree, min_percent=0.01, max_percent=99.99, log_a=1e2, stretch='log')(img_freefree)]).T.swapaxes(0,1) make_rgb(rgb, 'XRay_Halpha_freefree', layernames=('X-Ray', 'Halpha', 'freefree'), hsv_rotation=-0.25) #[Out]# target_header_cmz = fits.Header.fromstring(""" NAXIS = 2 NAXIS1 = 1920 NAXIS2 = 960 CTYPE1 = 'GLON-MOL' CRPIX1 = 960.5 CRVAL1 = 0.0 CDELT1 = -0.05 CUNIT1 = 'deg ' CTYPE2 = 'GLAT-MOL' CRPIX2 = 480.5 CRVAL2 = 0.0 CDELT2 = 0.05 CUNIT2 = 'deg ' COORDSYS= 'icrs ' """, sep='\n') target_header_innergal = fits.Header.fromstring(""" NAXIS = 2 NAXIS1 = 3840 NAXIS2 = 960 CTYPE1 = 'GLON-MOL' CRPIX1 = 1920.5 CRVAL1 = 0.0 CDELT1 = -0.05 CUNIT1 = 'deg ' CTYPE2 = 'GLAT-MOL' CRPIX2 = 480.5 CRVAL2 = 0.0 CDELT2 = 0.05 CUNIT2 = 'deg ' COORDSYS= 'icrs ' """, sep='\n') img_HI4PI_cmz, footprint = reproject_from_healpix(fn_HI4PI, target_header_cmz, field=5) img_ThermalDust_cmz, footprint = reproject_from_healpix(fn_ThermalDust, target_header_cmz) img_CO21_cmz, footprint = reproject_from_healpix(fn_CO21, target_header_cmz) img_HI4PI_innergal, footprint = reproject_from_healpix(fn_HI4PI, target_header_innergal, field=5) img_ThermalDust_innergal, footprint = reproject_from_healpix(fn_ThermalDust, target_header_innergal) img_CO21_innergal, footprint = reproject_from_healpix(fn_CO21, target_header_innergal) get_ipython().run_line_magic('run', 'ACES_Animation.py') rgb = np.array([simple_norm(img_HI4PI_cmz, min_percent=0.01, max_percent=99.99, log_a=1e1, stretch='log')(img_HI4PI_cmz), simple_norm(img_CO21_cmz, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21_cmz), simple_norm(img_ThermalDust_cmz, min_percent=0.01, max_percent=99.99, log_a=5e2, stretch='log')(img_ThermalDust_cmz)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header_cmz),) #frame_class=EllipticalFrame) ax.imshow(rgb_scaled) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) ax.set_xlim(240*2, 720*2) ax.set_ylim(120*2, 360*2) #[Out]# (240.0, 720.0) from astropy import units as u hi_cube_fn = '/orange/adamginsburg/cmz/hi/mcluregriffiths/GC.hi.tb.allgal.fits' #if os.path.exists("GC.hi.tb.allgal.fits.gz"): # hi_cube_fh = fits.open("GC.hi.tb.allgal.fits.gz") #else: # print("Not downloaded.") # #hi_cube_fh = download_file("http://www.atnf.csiro.au/research/HI/sgps/SGPS_1/GC.hi.tb.allgal.fits.gz") hi_cube = SpectralCube.read(hi_cube_fn, use_dask=True) peak_hi = hi_cube.max(axis=0) mean_hi = hi_cube.mean(axis=0) pl.imshow(peak_hi.value) #[Out]# pl.imshow(mean_hi.value, vmin=0) #[Out]# from matplotlib.colors import ListedColormap cmap = pl.cm.gray gray_transparent = cmap(np.arange(cmap.N)) gray_transparent[:,-1] = np.linspace(0, 1, cmap.N) gray_transparent = ListedColormap(gray_transparent) cmap = pl.cm.Oranges orange_transparent = cmap(np.arange(cmap.N)) orange_transparent[:,-1] = np.linspace(0, 1, cmap.N) orange_transparent = ListedColormap(orange_transparent) cmap = pl.cm.Oranges_r orange_transparent_r = cmap(np.arange(cmap.N)) orange_transparent_r[:,-1] = np.linspace(0, 1, cmap.N) orange_transparent_r = ListedColormap(orange_transparent_r) orange_transparent2 = cmap(np.arange(cmap.N)) orange_transparent2[:,-1] = np.linspace(0, 1, cmap.N) orange_transparent2[:,:-1] = [0.95, 0.35, 0.03] orange_transparent2 = ListedColormap(orange_transparent2) if os.path.exists("DHT02_Center_interp.fits"): cube_dame12co = SpectralCube.read('DHT02_Center_interp.fits', use_dask=True) else: pass peak_dameco = cube_dame12co.max(axis=0) mean_dameco = cube_dame12co.mean(axis=0) if os.path.exists('GalacticCenterMosaic_downsampled4x_2kms.fits'): co_cube = SpectralCube.read('GalacticCenterMosaic_downsampled4x_2kms.fits', use_dask=True) peak_co = co_cube.max(axis=0) mean_co = co_cube.mean(axis=0) rgb = np.array([simple_norm(img_HI4PI_cmz, min_percent=0.01, max_percent=99.99, log_a=1e1, stretch='log')(img_HI4PI_cmz), simple_norm(img_CO21_cmz, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21_cmz), simple_norm(img_ThermalDust_cmz, min_percent=0.01, max_percent=99.99, log_a=5e2, stretch='log')(img_ThermalDust_cmz)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header_cmz),) #frame_class=EllipticalFrame) ax.imshow(rgb_scaled) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) ax.set_xlim(240*2, 720*2) ax.set_ylim(120*2, 360*2) ax.imshow(mean_hi.value, transform=ax.get_transform(hi_cube.wcs.celestial), norm=simple_norm(mean_hi.value, min_percent=2, max_percent=99.9, stretch='asinh'), cmap=orange_transparent) #[Out]# rgb = np.array([simple_norm(img_HI4PI_cmz, min_percent=0.01, max_percent=99.99, log_a=1e1, stretch='log')(img_HI4PI_cmz), simple_norm(img_CO21_cmz, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21_cmz), simple_norm(img_ThermalDust_cmz, min_percent=0.01, max_percent=99.99, log_a=5e2, stretch='log')(img_ThermalDust_cmz)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header_cmz),) #frame_class=EllipticalFrame) ax.imshow(rgb_scaled) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) ax.set_xlim(240*2, 720*2) ax.set_ylim(120*2, 360*2) ax.imshow(mean_dameco.value, norm=simple_norm(mean_dameco.value, min_percent=2, max_percent=99.9, stretch='asinh'), transform=ax.get_transform(cube_dame12co.wcs.celestial), cmap=orange_transparent) pl.savefig("HI4PI_CO_Dust_withDameCOOrangeOverlay.png", bbox_inches='tight') rgb = np.array([simple_norm(img_HI4PI_cmz, min_percent=0.01, max_percent=99.99, log_a=1e1, stretch='log')(img_HI4PI_cmz), simple_norm(img_CO21_cmz, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21_cmz), simple_norm(img_ThermalDust_cmz, min_percent=0.01, max_percent=99.99, log_a=5e2, stretch='log')(img_ThermalDust_cmz)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header_cmz),) #frame_class=EllipticalFrame) ax.imshow(rgb_scaled) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) ax.set_xlim(240*2, 720*2) ax.set_ylim(120*2, 360*2) ax.imshow(peak_co.value, norm=simple_norm(peak_co.value, min_percent=0.1, max_percent=99.9, stretch='asinh'), transform=ax.get_transform(co_cube.wcs.celestial), cmap=orange_transparent) pl.savefig("HI4PI_CO_Dust_withJCMTCOOrange.png", bbox_inches='tight') atlasgal = fits.open('/orange/adamginsburg/galactic_plane_surveys/atlasgal/MOSAICS/apex_planck.fits') print(atlasgal[0].data.shape) rgb = np.array([simple_norm(img_HI4PI_cmz, min_percent=0.01, max_percent=99.99, log_a=1e1, stretch='log')(img_HI4PI_cmz), simple_norm(img_CO21_cmz, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21_cmz), simple_norm(img_ThermalDust_cmz, min_percent=0.01, max_percent=99.99, log_a=5e2, stretch='log')(img_ThermalDust_cmz)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header_cmz),) #frame_class=EllipticalFrame) ax.imshow(rgb_scaled) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) ax.set_xlim(240*2, 720*2) ax.set_ylim(120*2, 360*2) ax.imshow(atlasgal[0].data, norm=simple_norm(atlasgal[0].data, min_percent=0.1, max_percent=99.9, stretch='asinh'), transform=ax.get_transform(WCS(atlasgal[0].header)), cmap='gray') pl.savefig("HI4PI_CO_Dust_withATLASGAL.png", bbox_inches='tight') rgb = np.array([simple_norm(img_HI4PI_cmz, min_percent=0.01, max_percent=99.99, log_a=1e1, stretch='log')(img_HI4PI_cmz), simple_norm(img_CO21_cmz, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21_cmz), simple_norm(img_ThermalDust_cmz, min_percent=0.01, max_percent=99.99, log_a=5e2, stretch='log')(img_ThermalDust_cmz)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header_cmz),) #frame_class=EllipticalFrame) ax.imshow(rgb_scaled) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) ax.set_xlim(240*2, 720*2) ax.set_ylim(120*2, 360*2) ax.imshow(atlasgal[0].data, norm=simple_norm(atlasgal[0].data, min_percent=5, max_percent=99.9, stretch='asinh'), transform=ax.get_transform(WCS(atlasgal[0].header)), cmap=gray_transparent) pl.savefig("HI4PI_CO_Dust_withATLASGALtransparent.png", bbox_inches='tight') rgb = np.array([simple_norm(img_HI4PI_cmz, min_percent=0.01, max_percent=99.99, log_a=1e1, stretch='log')(img_HI4PI_cmz), simple_norm(img_CO21_cmz, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21_cmz), simple_norm(img_ThermalDust_cmz, min_percent=0.01, max_percent=99.99, log_a=5e2, stretch='log')(img_ThermalDust_cmz)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header_cmz),) #frame_class=EllipticalFrame) ax.imshow(rgb_scaled) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) ax.set_xlim(240*2, 720*2) ax.set_ylim(120*2, 360*2) ax.imshow(atlasgal[0].data, norm=simple_norm(atlasgal[0].data, min_percent=5, max_percent=99.9, stretch='asinh'), transform=ax.get_transform(WCS(atlasgal[0].header)), cmap=orange_transparent) pl.savefig("HI4PI_CO_Dust_withATLASGALOrangetransparent.png", bbox_inches='tight') from astroquery.vizier import Vizier from astropy.coordinates import SkyCoord almaimf_cat = Vizier.get_catalogs('J/A+A/662/A9/objects')[0] almaimf_crds = SkyCoord(almaimf_cat['RAJ2000'], almaimf_cat['DEJ2000'], frame='fk5', unit=(u.h, u.deg)) rgb = np.array([simple_norm(img_HI4PI_cmz, min_percent=0.01, max_percent=99.99, log_a=1e1, stretch='log')(img_HI4PI_cmz), simple_norm(img_CO21_cmz, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21_cmz), simple_norm(img_ThermalDust_cmz, min_percent=0.01, max_percent=99.99, log_a=5e2, stretch='log')(img_ThermalDust_cmz)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header_cmz),) #frame_class=EllipticalFrame) ax.imshow(rgb_scaled) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) ax.set_xlim(240*2, 720*2) ax.set_ylim(120*2, 360*2) ax.imshow(atlasgal[0].data, norm=simple_norm(atlasgal[0].data, min_percent=5, max_percent=99.9, stretch='asinh'), transform=ax.get_transform(WCS(atlasgal[0].header)), cmap=orange_transparent) ax.scatter(almaimf_crds.galactic.l, almaimf_crds.galactic.b, marker='*', edgecolor='purple', facecolor='none', transform=ax.get_transform('galactic')) pl.savefig("HI4PI_CO_Dust_withATLASGALOrangetransparent_ALMA-IMFtargets.png", bbox_inches='tight') rgb = np.array([simple_norm(img_HI4PI_innergal, min_percent=0.01, max_percent=99.99, log_a=1e1, stretch='log')(img_HI4PI_innergal), simple_norm(img_CO21_innergal, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21_innergal), simple_norm(img_ThermalDust_innergal, min_percent=0.01, max_percent=99.99, log_a=5e2, stretch='log')(img_ThermalDust_innergal)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header_innergal),) #frame_class=EllipticalFrame) ax.imshow(rgb_scaled) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) ax.set_xlim(540, 3840-540) ax.set_ylim(0, 960) ax.imshow(atlasgal[0].data, norm=simple_norm(atlasgal[0].data, min_percent=5, max_percent=99.9, stretch='asinh'), transform=ax.get_transform(WCS(atlasgal[0].header)), cmap=orange_transparent) pl.savefig("HI4PI_CO_Dust_innergal_withATLASGALOrangetransparent.png", bbox_inches='tight') ax.scatter(almaimf_crds.galactic.l, almaimf_crds.galactic.b, marker='*', edgecolor='purple', facecolor='none', transform=ax.get_transform('galactic')) pl.savefig("HI4PI_CO_Dust_innergal_withATLASGALOrangetransparent_ALMA-IMFtargets.png", bbox_inches='tight') import regions from astropy import coordinates rect = regions.RectangleSkyRegion( coordinates.SkyCoord(0*u.deg, 0*u.deg, frame='galactic'), width=np.abs((ax.get_xlim()[1] - ax.get_xlim()[0])*target_header_innergal['CDELT1'])*u.deg, height=(ax.get_ylim()[1] - ax.get_ylim()[0])*target_header_innergal['CDELT2']*u.deg, ) # different colors, overwrites previous rgb = np.array([simple_norm(img_HI4PI, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PI), simple_norm(img_CO21, min_percent=0.01, max_percent=99.90, log_a=2e1, stretch='log')(img_CO21), simple_norm(img_ThermalDust, min_percent=0.01, max_percent=99.90, log_a=5e2, stretch='log')(img_ThermalDust)]).T.swapaxes(0,1) ax = make_rgb(rgb, 'HI-CO-Dust_m0.35', -0.35, layernames=('HI', 'CO', 'Dust'), do_layers=False) prect = rect.to_pixel(ax.wcs) prect.visual['edgecolor'] = 'w' innergalplot = prect.plot(ax=ax) pl.savefig("HI-CO-Dust_m0.35_RGB_innergal_zoombox.png", bbox_inches='tight') innergaltext = ax.text(0, 25, "Galactic Plane", transform=ax.get_transform('world'), color='w', horizontalalignment='center') cmzrect = regions.RectangleSkyRegion( coordinates.SkyCoord(0*u.deg, 0*u.deg, frame='galactic'), width=12*u.deg, height=4.8*u.deg, ) pcmzrect = cmzrect.to_pixel(ax.wcs) pcmzrect.visual['edgecolor'] = 'w' cmzplot = pcmzrect.plot(ax=ax) cmztext = ax.text(0, 5, "Central Molecular Zone", transform=ax.get_transform('world'), color='w', horizontalalignment='center') orionrect = regions.RectangleSkyRegion( coordinates.SkyCoord(209*u.deg, -19*u.deg, frame='galactic'), width=15*u.deg, height=15*u.deg, ) porionrect = orionrect.to_pixel(ax.wcs) porionrect.visual['edgecolor'] = 'w' oriplot = porionrect.plot(ax=ax) oritext = ax.text(209, -7, "Orion", transform=ax.get_transform('world'), color='w', horizontalalignment='center') pl.savefig("HI-CO-Dust_m0.35_RGB_labeled_zones.png", bbox_inches='tight') innergalplot.set_visible(False) cmzplot.set_visible(False) innergaltext.set_visible(False) cmztext.set_visible(False) pl.savefig("HI-CO-Dust_m0.35_RGB_labeled_orion.png", bbox_inches='tight') oritext.set_visible(False) oriplot.set_visible(False) m31 = coordinates.SkyCoord.from_name('M31') andromedarect = regions.RectangleSkyRegion( m31, width=5*u.deg, height=5*u.deg, ) pandromedarect = andromedarect.to_pixel(ax.wcs) pandromedarect.visual['edgecolor'] = 'w' andromedaplot = pandromedarect.plot(ax=ax) andromedatext = ax.text(m31.galactic.l.deg, m31.galactic.b.deg+4, "Andromeda", transform=ax.get_transform('world'), color='w', horizontalalignment='center') pl.savefig("HI-CO-Dust_m0.35_RGB_labeled_andromeda.png", bbox_inches='tight') andromedatext.set_visible(False) andromedaplot.set_visible(False) taurus = coordinates.SkyCoord.from_name('taurus') taurusrect = regions.RectangleSkyRegion( taurus, width=15*u.deg, height=15*u.deg, ) ptaurusrect = taurusrect.to_pixel(ax.wcs) ptaurusrect.visual['edgecolor'] = 'w' taurusplot = ptaurusrect.plot(ax=ax) taurustext = ax.text(taurus.galactic.l.deg-5, taurus.galactic.b.deg+10, "Taurus", transform=ax.get_transform('world'), color='w', horizontalalignment='center') pl.savefig("HI-CO-Dust_m0.35_RGB_labeled_taurus.png", bbox_inches='tight') taurustext.set_visible(False) taurusplot.set_visible(False) # different colors, overwrites previous rgb = np.array([simple_norm(img_HI4PI_orion, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PIorion), simple_norm(img_CO21orion, min_percent=0.01, max_percent=99.90, log_a=2e1, stretch='log')(img_CO21orion), simple_norm(img_ThermalDust_orion, min_percent=0.01, max_percent=99.90, log_a=5e2, stretch='log')(img_ThermalDust_orion)]).T.swapaxes(0,1) ax = make_rgb(rgb, 'HI-CO-Dust_m0.35_orion', -0.35, header=target_header_orion, layernames=('HI', 'CO', 'Dust'), do_layers=False) orionrect = regions.RectangleSkyRegion( coordinates.SkyCoord(209*u.deg, -19*u.deg, frame='galactic'), width=15*u.deg, height=15*u.deg, ) porionrect = orionrect.to_pixel(ax.wcs) porionrect.visual['edgecolor'] = 'w' oriplot = porionrect.plot(ax=ax) oritext = ax.text(83.82, 5, "Orion", transform=ax.get_transform('world'), color='w', horizontalalignment='center') pl.savefig("HI-CO-Dust_m0.35_orion_RGB_labeled_orion.png", bbox_inches='tight') # different colors, overwrites previous rgb = np.array([simple_norm(img_HI4PI_orion, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PI_orion), simple_norm(img_CO21_orion, min_percent=0.01, max_percent=99.90, log_a=2e1, stretch='log')(img_CO21_orion), simple_norm(img_ThermalDust_orion, min_percent=0.01, max_percent=99.90, log_a=5e2, stretch='log')(img_ThermalDust_orion)]).T.swapaxes(0,1) ax = make_rgb(rgb, 'HI-CO-Dust_m0.35_orion', -0.35, header=target_header_orion, layernames=('HI', 'CO', 'Dust'), do_layers=False) orionrect = regions.RectangleSkyRegion( coordinates.SkyCoord(209*u.deg, -19*u.deg, frame='galactic'), width=15*u.deg, height=15*u.deg, ) porionrect = orionrect.to_pixel(ax.wcs) porionrect.visual['edgecolor'] = 'w' oriplot = porionrect.plot(ax=ax) oritext = ax.text(83.82, 5, "Orion", transform=ax.get_transform('world'), color='w', horizontalalignment='center') pl.savefig("HI-CO-Dust_m0.35_orion_RGB_labeled_orion.png", bbox_inches='tight') target_header_oricloud = fits.Header.fromstring(""" NAXIS = 2 NAXIS1 = 1920 NAXIS2 = 1920 CTYPE1 = 'RA---TAN' CRPIX1 = 960.5 CRVAL1 = 83.82 CDELT1 = -0.05 CUNIT1 = 'deg ' CTYPE2 = 'DEC--TAN' CRPIX2 = 960.5 CRVAL2 = -5.39 CDELT2 = 0.05 CUNIT2 = 'deg ' COORDSYS= 'icrs ' """, sep='\n') img_HI4PI_oricloud, footprint = reproject_from_healpix(fn_HI4PI, target_header_oricloud, field=5) img_ThermalDust_oricloud, footprint = reproject_from_healpix(fn_ThermalDust, target_header_oricloud) img_CO21_oricloud, footprint = reproject_from_healpix(fn_CO21, target_header_oricloud) # different colors, overwrites previous rgb = np.array([simple_norm(img_HI4PI, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PI_oricloud), simple_norm(img_CO21, min_percent=0.01, max_percent=99.90, log_a=2e1, stretch='log')(img_CO21_oricloud), simple_norm(img_ThermalDust, min_percent=0.01, max_percent=99.90, log_a=5e2, stretch='log')(img_ThermalDust_oricloud)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) plt.figure(figsize=(10,5), dpi=200, frameon=False) ax = plt.subplot(1,1,1, projection=WCS(target_header_oricloud),) ax.imshow(rgb_scaled) ax.axis('off') oricloudrect = regions.RectangleSkyRegion( coordinates.SkyCoord(209*u.deg, -19*u.deg, frame='galactic'), width=15*u.deg, height=15*u.deg, ) poricloudrect = oricloudrect.to_pixel(ax.wcs) poricloudrect.visual['edgecolor'] = 'w' oriplot = poricloudrect.plot(ax=ax) oritext = ax.text(80, 4, "Orion", transform=ax.get_transform('world'), color='w', horizontalalignment='center') ax.axis([480,1440,480,1440]) pl.savefig("HI-CO-Dust_m0.35_oricloud_RGB_labeled_oricloud.png", bbox_inches='tight', pad_inches=0) target_header_omc = fits.Header.fromstring(""" NAXIS = 2 NAXIS1 = 1920 NAXIS2 = 1920 CTYPE1 = 'RA---TAN' CRPIX1 = 960.5 CRVAL1 = 83.82 CDELT1 = -0.0005 CUNIT1 = 'deg ' CTYPE2 = 'DEC--TAN' CRPIX2 = 960.5 CRVAL2 = -5.39 CDELT2 = 0.0005 CUNIT2 = 'deg ' COORDSYS= 'icrs ' """, sep='\n') 1920*0.0005 #[Out]# 0.96 fn_carmaorion_12co = "mom0_12co_pix_2_Tmb.fits" fn_carmaorion_13co = "mom0_13co_pix_2_Tmb.fits" if not os.path.exists(fn_carmaorion_12co): dl=download_file("https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/6Q26PN/MRUINN") shutil.move(dl, fn_carmaorion_12co) dl=download_file("https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/6Q26PN/AIYBES") shutil.move(dl, fn_carmaorion_13co) fn_orion_mustang = "MUSTANG_OrionClean9.0-24nov08.fits" if not os.path.exists(fn_orion_mustang): dl = download_file("https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/LZELK1/MZNCT1") shutil.move(dl, fn_orion_mustang) fn_orion_bolocam = 'orionfruitawf2.fits' img_carmaori_12co_omc, footprint = reproject.reproject_interp( (fits.getdata(fn_carmaorion_12co).squeeze(), WCS(fits.getheader(fn_carmaorion_12co)).celestial), target_header_omc) img_carmaori_13co_omc, footprint = reproject.reproject_interp( (fits.getdata(fn_carmaorion_13co).squeeze(), WCS(fits.getheader(fn_carmaorion_13co)).celestial), target_header_omc) #img_orion_mustang_omc, footprint = reproject.reproject_interp( # (fits.getdata(fn_orion_mustang).squeeze(), # WCS(fits.getheader(fn_orion_mustang)).celestial), # target_header_omc) img_orion_bolocam_omc, footprint = reproject.reproject_interp( (fits.getdata(fn_orion_bolocam).squeeze(), WCS(fits.getheader(fn_orion_bolocam)).celestial), target_header_omc) pl.imshow(img_carmaori_12co_omc) #[Out]# img_orion_bolocam_omc[img_orion_bolocam_omc<-0.5] = 0 pl.imshow(img_orion_bolocam_omc, norm=simple_norm(img_orion_bolocam_omc, min_percent=1, max_percent=99.95, log_a=5e2, stretch='log')) #[Out]# rgb = np.array([simple_norm(img_carmaori_12co_omc, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_carmaori_12co_omc), simple_norm(img_carmaori_13co_omc, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_carmaori_13co_omc), simple_norm(img_orion_bolocam_omc, min_percent=1, max_percent=99.95, log_a=5e2, stretch='log')(img_orion_bolocam_omc)]).T.swapaxes(0,1) pl.imshow(rgb) #[Out]# # ALMA imaging - not that useful fn_orion4d_cont = "Orion4D_ISF_Continuum_CLEAN_b0.5_30kltaper_DEEP.image.pbcor.fits" if not os.path.exists(fn_orion4d_cont): dl = download_file("https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/I1GQER/DHGPFZ") shutil.move(dl, fn_orion4d_cont) fn_trapezium = 'Trapezium_GEMS_mosaic_redblueorange_normed_small_contrast_bright.png' if not os.path.exists(fn_trapezium): dl = download_file(f'https://keflavich.github.io/talks/assets_orion/{fn_trapezium}') shutil.move(dl, fn_trapezium) # different colors, overwrites previous rgb = np.array([simple_norm(img_HI4PI, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PI_oricloud), simple_norm(img_CO21, min_percent=0.01, max_percent=99.90, log_a=2e1, stretch='log')(img_CO21_oricloud), simple_norm(img_ThermalDust, min_percent=0.01, max_percent=99.90, log_a=5e2, stretch='log')(img_ThermalDust_oricloud)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) plt.figure(figsize=(10,5), dpi=200, frameon=False) ax = plt.subplot(1,2,1, projection=WCS(target_header_oricloud),) ax.imshow(rgb_scaled) ax.axis('off') ax.axis([600,1280,600,1280]) ax.axis([720,1160,720,1160]) ax.axis([800,1080,800,1080]) rgb = np.array([simple_norm(img_carmaori_12co_omc, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_carmaori_12co_omc), simple_norm(img_carmaori_13co_omc, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_carmaori_13co_omc), simple_norm(img_orion_bolocam_omc, min_percent=1, max_percent=99.5, log_a=5e2, stretch='log')(img_orion_bolocam_omc)]).T.swapaxes(0,1) ww = WCS(target_header_omc) axins = inset_axes(ax, width=4, height=4, bbox_to_anchor=[0.0, 0.0, 0.9, 1], bbox_transform=fig.transFigure, loc='center right', axes_class=astropy.visualization.wcsaxes.core.WCSAxes, axes_kwargs=dict(wcs=ww)) axins.imshow(rgb) axins.axis('off') mark_inset_generic(axins, ax, data=rgb[:,:,0], loc2=2, loc1=3, edgecolor='white') pl.savefig("HI-CO-Dust_m0.35_oricloud-rezoom_RGB.png", bbox_inches='tight', pad_inches=0) from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, inset_axes, mark_inset rgb = np.array([simple_norm(img_HI4PI, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PI_oricloud), simple_norm(img_CO21, min_percent=0.01, max_percent=99.90, log_a=2e1, stretch='log')(img_CO21_oricloud), simple_norm(img_ThermalDust, min_percent=0.01, max_percent=99.90, log_a=5e2, stretch='log')(img_ThermalDust_oricloud)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) plt.figure(figsize=(10,5), dpi=200, frameon=False) ax = plt.subplot(1,2,1, projection=WCS(target_header_oricloud),) ax.imshow(rgb_scaled) ax.axis('off') ax.axis([600,1280,600,1280]) ax.axis([720,1160,720,1160]) ax.axis([800,1080,800,1080]) rgb = np.array([simple_norm(img_carmaori_12co_omc, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_carmaori_12co_omc), simple_norm(img_carmaori_13co_omc, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_carmaori_13co_omc), simple_norm(img_orion_bolocam_omc, min_percent=1, max_percent=99.5, log_a=5e2, stretch='log')(img_orion_bolocam_omc)]).T.swapaxes(0,1) ww = WCS(target_header_omc) axins = inset_axes(ax, width=4, height=4, bbox_to_anchor=[0.0, 0.0, 0.9, 1], bbox_transform=fig.transFigure, loc='center right', axes_class=astropy.visualization.wcsaxes.core.WCSAxes, axes_kwargs=dict(wcs=ww)) axins.imshow(rgb) axins.axis('off') mark_inset_generic(axins, ax, data=rgb[:,:,0], loc2=2, loc1=3, edgecolor='white') pl.savefig("HI-CO-Dust_m0.35_oricloud-rezoom_RGB.png", bbox_inches='tight', pad_inches=0) rgb = np.array([simple_norm(img_HI4PI, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PI_oricloud), simple_norm(img_CO21, min_percent=0.01, max_percent=99.90, log_a=2e1, stretch='log')(img_CO21_oricloud), simple_norm(img_ThermalDust, min_percent=0.01, max_percent=99.90, log_a=5e2, stretch='log')(img_ThermalDust_oricloud)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) fig =plt.figure(figsize=(10,5), dpi=200, frameon=False) ax = plt.subplot(1,2,1, projection=WCS(target_header_oricloud),) ax.imshow(rgb_scaled) ax.axis('off') ax.axis([600,1280,600,1280]) ax.axis([720,1160,720,1160]) ax.axis([800,1080,800,1080]) rgb = np.array([simple_norm(img_carmaori_12co_omc, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_carmaori_12co_omc), simple_norm(img_carmaori_13co_omc, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_carmaori_13co_omc), simple_norm(img_orion_bolocam_omc, min_percent=1, max_percent=99.5, log_a=5e2, stretch='log')(img_orion_bolocam_omc)]).T.swapaxes(0,1) ww = WCS(target_header_omc) axins = inset_axes(ax, width=4, height=4, bbox_to_anchor=[0.0, 0.0, 0.9, 1], bbox_transform=fig.transFigure, loc='center right', axes_class=astropy.visualization.wcsaxes.core.WCSAxes, axes_kwargs=dict(wcs=ww)) axins.imshow(rgb) axins.axis('off') mark_inset_generic(axins, ax, data=rgb[:,:,0], loc2=2, loc1=3, edgecolor='white') pl.savefig("HI-CO-Dust_m0.35_oricloud-rezoom_RGB.png", bbox_inches='tight', pad_inches=0) from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, inset_axes, mark_inset from mpl_toolkits.axes_grid1.inset_locator import TransformedBbox, BboxPatch, BboxConnector, Path from astropy import wcs from matplotlib.transforms import Bbox from matplotlib.patches import Polygon, PathPatch, ConnectionPatch def mark_inset_generic(axins, parent_ax, data, loc1=1, loc2=3, loc1in=None, loc2in=None, edgecolor='b', zorder1=100, zorder2=100, polyzorder=1): bl = axins.wcs.pixel_to_world(0, 0) br = axins.wcs.pixel_to_world(data.shape[1], 0) tl = axins.wcs.pixel_to_world(0, data.shape[0]) # x,y not y,x tr = axins.wcs.pixel_to_world(data.shape[1], data.shape[0]) # x,y not y,x fig = parent_ax.get_figure() #frame = parent_ax.wcs.wcs.radesys.lower() #frame = ax.wcs.world_axis_physical_types[0].split(".")[1] frame = wcs.utils.wcs_to_celestial_frame(ax.wcs) blt = bl.transform_to(frame) brt = br.transform_to(frame) tlt = tl.transform_to(frame) trt = tr.transform_to(frame) xys = [parent_ax.wcs.wcs_world2pix([[crd.spherical.lon.deg, crd.spherical.lat.deg]],0)[0] for crd in (trt, tlt, blt, brt, trt)] markinkwargs = dict(fc='none', ec=edgecolor) ppoly = Polygon(xys, fill=False, zorder=polyzorder, **markinkwargs) parent_ax.add_patch(ppoly) corners = parent_ax.transData.inverted().transform(xys) axcorners = [(1,1), (0,1), (0,0), (1,0)] corners = [(crd.spherical.lon.deg, crd.spherical.lat.deg) for crd in (trt, tlt, blt, brt)] if loc1in is None: loc1in = loc1 if loc2in is None: loc2in = loc2 con1 = ConnectionPatch(xyA=axcorners[loc1-1], coordsA='axes fraction', axesA=axins, xyB=corners[loc1in-1], coordsB=parent_ax.get_transform('world'), axesB=parent_ax, linestyle='-', color=edgecolor, zorder=zorder1) con2 = ConnectionPatch(xyA=axcorners[loc2-1], coordsA='axes fraction', axesA=axins, xyB=corners[loc2in-1], coordsB=parent_ax.get_transform('world'), axesB=parent_ax, linestyle='-', color=edgecolor, zorder=zorder2) fig.add_artist(con1) fig.add_artist(con2) return con1, con2 rgb = np.array([simple_norm(img_HI4PI, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PI_oricloud), simple_norm(img_CO21, min_percent=0.01, max_percent=99.90, log_a=2e1, stretch='log')(img_CO21_oricloud), simple_norm(img_ThermalDust, min_percent=0.01, max_percent=99.90, log_a=5e2, stretch='log')(img_ThermalDust_oricloud)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) fig =plt.figure(figsize=(10,5), dpi=200, frameon=False) ax = plt.subplot(1,2,1, projection=WCS(target_header_oricloud),) ax.imshow(rgb_scaled) ax.axis('off') ax.axis([600,1280,600,1280]) ax.axis([720,1160,720,1160]) ax.axis([800,1080,800,1080]) rgb = np.array([simple_norm(img_carmaori_12co_omc, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_carmaori_12co_omc), simple_norm(img_carmaori_13co_omc, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_carmaori_13co_omc), simple_norm(img_orion_bolocam_omc, min_percent=1, max_percent=99.5, log_a=5e2, stretch='log')(img_orion_bolocam_omc)]).T.swapaxes(0,1) ww = WCS(target_header_omc) axins = inset_axes(ax, width=4, height=4, bbox_to_anchor=[0.0, 0.0, 0.9, 1], bbox_transform=fig.transFigure, loc='center right', axes_class=astropy.visualization.wcsaxes.core.WCSAxes, axes_kwargs=dict(wcs=ww)) axins.imshow(rgb) axins.axis('off') mark_inset_generic(axins, ax, data=rgb[:,:,0], loc2=2, loc1=3, edgecolor='white') pl.savefig("HI-CO-Dust_m0.35_oricloud-rezoom_RGB.png", bbox_inches='tight', pad_inches=0) import regions from astropy import coordinates from astropy import units as u, constants import astropy.visualization.wcsaxes from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, inset_axes, mark_inset aces12m = fits.open('/orange/adamginsburg/ACES/mosaics/12m_continuum_mosaic.fits') aces12mwcs = WCS(aces12m[0].header) aces12m[0].data[aces12m[0].data < -0.0005] = 0 fig = plt.figure(figsize=(10,5), frameon=False) ax = plt.subplot(1,1,1, projection=aces12mwcs) ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) ax.axis('off') sgrb2m = regions.CircleSkyRegion( coordinates.SkyCoord(000.6672*u.deg, -00.03640*u.deg, frame='galactic'), radius=25*u.arcsec) sgrb2n = regions.CircleSkyRegion( coordinates.SkyCoord(000.6773*u.deg, -00.029*u.deg, frame='galactic'), radius=25*u.arcsec) sgrb2m.visual['edgecolor'] = 'r' sgrb2n.visual['edgecolor'] = 'r' merge = psgrb2n & psgrb2m merge_mask = merge.to_mask() slcs_merge,_ = merge_mask.get_overlap_slices(aces12m[0].data.shape) # moved down here for debugging... ax.imshow(aces12m[0].data, norm=simple_norm(aces12m[0].data, stretch='log', min_percent=None, max_percent=None, min_cut=-0.0005, max_cut=0.05,), cmap='gray', ) psgrb2m = sgrb2m.to_pixel(ax.wcs) point_sgrb2m = psgrb2m.plot(ax=ax) msgrb2m = psgrb2m.to_mask() slcs_sgrb2m,_ = msgrb2m.get_overlap_slices(aces12m[0].data.shape) bbox = [0.01, 0, 0.2, 1.05] axins1 = zoomed_inset_axes(ax, 30, loc='upper left', bbox_to_anchor=bbox, bbox_transform=fig.transFigure, axes_class=astropy.visualization.wcsaxes.core.WCSAxes, axes_kwargs=dict(wcs=ax.wcs) ) axins1.axis('off') axins1.imshow(aces12m[0].data, norm=simple_norm(aces12m[0].data, stretch='log', min_percent=None, max_percent=None, min_cut=-0.0005, max_cut=0.5,), cmap='gray' ) axins1.axis(msgrb2m.bbox.extent) mark1 = mark_inset(ax, axins1, loc1=1, loc2=3, edgecolor='r') point_sgrb2m_in = sgrb2m.to_pixel(axins1.wcs).plot(ax=axins1) axins1.coords['glat'].set_ticklabel(visible=False) axins1.coords['glon'].set_ticklabel(visible=False) axins1.coords['glat'].set_ticks_visible(False) axins1.coords['glon'].set_ticks_visible(False) pl.savefig("aces_mark_sgrb2m.png", bbox_inches='tight', dpi=200) point_sgrb2m.set_visible(False) point_sgrb2m_in.set_visible(False) #axins1.set_visible(False) for pp in mark1: pp.set_visible(False) psgrb2n = sgrb2n.to_pixel(ax.wcs) psgrb2n.visual['edgecolor'] = 'r' point_sgrb2n = psgrb2n.plot(ax=ax) msgrb2n = psgrb2n.to_mask() slcs_sgrb2n,_ = msgrb2n.get_overlap_slices(aces12m[0].data.shape) #axins2 = zoomed_inset_axes(ax, 30, loc='upper left', bbox_to_anchor=bbox, bbox_transform=fig.transFigure, # axes_class=astropy.visualization.wcsaxes.core.WCSAxes, # axes_kwargs=dict(wcs=ax.wcs)) axins1.imshow(aces12m[0].data, norm=simple_norm(aces12m[0].data, stretch='log', min_percent=None, max_percent=None, min_cut=-0.0005, max_cut=0.5,), cmap='gray' ) axins1.axis(msgrb2n.bbox.extent) mark2 = mark_inset(ax, axins1, loc1=1, loc2=3, edgecolor='r') point_sgrb2n_in = sgrb2n.to_pixel(axins1.wcs).plot(ax=axins1) #axins2.coords['glat'].set_ticklabel(visible=False) #axins2.coords['glon'].set_ticklabel(visible=False) #axins2.coords['glat'].set_ticks_visible(False) #axins2.coords['glon'].set_ticks_visible(False) pl.savefig("aces_mark_sgrb2n.png", bbox_inches='tight', dpi=200) point_sgrb2n.set_visible(False) point_sgrb2n_in.set_visible(False) for pp in mark2: pp.set_visible(False) axins1.set_visible(False) axins3 = zoomed_inset_axes(ax, 15, loc='upper left', bbox_to_anchor=[0.05, 0, 0.2, 1], bbox_transform=fig.transFigure, axes_class=astropy.visualization.wcsaxes.core.WCSAxes, axes_kwargs=dict(wcs=ax.wcs)) axins3.axis('off') axins3.imshow(aces12m[0].data, norm=simple_norm(aces12m[0].data, stretch='log', min_percent=None, max_percent=None, min_cut=-0.0005, max_cut=0.5,), cmap='gray' ) axins3.axis(merge_mask.bbox.extent) mark3 = mark_inset(ax, axins3, loc1=1, loc2=3, edgecolor='r') point_sgrb2n_in = sgrb2n.to_pixel(axins3.wcs).plot(ax=axins3) point_sgrb2m_in = sgrb2m.to_pixel(axins3.wcs).plot(ax=axins3) axins3.coords['glat'].set_ticklabel(visible=False) axins3.coords['glon'].set_ticklabel(visible=False) axins3.coords['glat'].set_ticks_visible(False) axins3.coords['glon'].set_ticks_visible(False) pl.savefig("aces_mark_sgrb2mANDn.png", bbox_inches='tight', dpi=200) # different colors, overwrites previous rgb = np.array([simple_norm(img_HI4PI, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PI_oricloud), simple_norm(img_CO21, min_percent=0.01, max_percent=99.90, log_a=2e1, stretch='log')(img_CO21_oricloud), simple_norm(img_ThermalDust, min_percent=0.01, max_percent=99.90, log_a=5e2, stretch='log')(img_ThermalDust_oricloud)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) fig = plt.figure(figsize=(10,5), dpi=200, frameon=False) ax = plt.subplot(1,1,1, projection=WCS(target_header_oricloud),) ax.imshow(rgb_scaled) ax.axis('off') ax.axis([600,1280,600,1280]) ax.axis([720,1160,720,1160]) rgb = np.array(PIL.Image.open(fn_trapezium))[::-1, :, :] avm = pyavm.AVM.from_image(fn_trapezium) ww = avm.to_wcs() axins = inset_axes(ax, width=4, height=4, bbox_to_anchor=[0.8, 0.0, 0.3, 1], bbox_transform=fig.transFigure, loc='center right', axes_class=astropy.visualization.wcsaxes.core.WCSAxes, axes_kwargs=dict(wcs=ww)) axins.imshow(rgb) axins.coords['ra'].set_ticklabel(visible=False) axins.coords['dec'].set_ticklabel(visible=False) axins.coords['ra'].set_ticks_visible(False) axins.coords['dec'].set_ticks_visible(False) mark_inset_generic(axins, ax, data=rgb[:,:,0], loc2=2, loc1=3, edgecolor='w') pl.savefig("HI-CO-Dust_m0.35_oricloud-rezoom-trapezium_RGB.png", bbox_inches='tight', pad_inches=0) theta1c = SkyCoord.from_name('* tet01 Ori C') theta1c #[Out]# trapezium_header = """ WCSAXES = 2 / Number of coordinate axes CRPIX1 = 1020.4109 / Pixel coordinate of reference point CRPIX2 = 1020.4288 / Pixel coordinate of reference point PC1_1 = -3.1E-05 / Coordinate transformation matrix element PC2_2 = 3.1E-05 / Coordinate transformation matrix element CDELT1 = 1.0 / [deg] Coordinate increment at reference point CDELT2 = 1.0 / [deg] Coordinate increment at reference point CUNIT1 = 'deg' / Units of coordinate increment and value CUNIT2 = 'deg' / Units of coordinate increment and value CTYPE1 = 'RA---TAN' / Right ascension, gnomonic projection CTYPE2 = 'DEC--TAN' / Declination, gnomonic projection CRVAL1 = 83.81251 / [deg] Coordinate value at reference point CRVAL2 = -5.3631591 / [deg] Coordinate value at reference point LONPOLE = 180.0 / [deg] Native longitude of celestial pole LATPOLE = -5.3631591 / [deg] Native latitude of celestial pole MJDREF = 0.0 / [d] MJD of fiducial time RADESYS = 'FK5' / Equatorial coordinate system EQUINOX = 2000.0 / [yr] Equinox of equatorial coordinates """ rgb = np.array([simple_norm(img_carmaori_12co_omc, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_carmaori_12co_omc), simple_norm(img_carmaori_13co_omc, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_carmaori_13co_omc), simple_norm(img_orion_bolocam_omc, min_percent=1, max_percent=99.5, log_a=5e2, stretch='log')(img_orion_bolocam_omc)]).T.swapaxes(0,1) ww = WCS(target_header_omc) plt.figure(figsize=(10,5), dpi=200, frameon=False) ax = plt.subplot(1,2,1, projection=ww,) ax.imshow(rgb) ax.axis('off') ax.axis([480, 1440, 480, 1440]) #ax.scatter(theta1c.ra, theta1c.dec, transform=ax.get_transform('world'), # marker='*', color='y') rgb = np.array(PIL.Image.open(fn_trapezium))[::-1, :, :] avm = pyavm.AVM.from_image(fn_trapezium) ww = avm.to_wcs() ww = WCS(fits.Header.fromstring(trapezium_header, sep='\n')) axins = inset_axes(ax, width=4, height=4, bbox_to_anchor=[0.0, 0.0, 0.9, 1], bbox_transform=fig.transFigure, loc='center right', axes_class=astropy.visualization.wcsaxes.core.WCSAxes, axes_kwargs=dict(wcs=ww)) axins.imshow(rgb) axins.axis('off') #axins.scatter(theta1c.ra, theta1c.dec, transform=axins.get_transform('world'), # marker='*', color='y') mark_inset_generic(axins, ax, data=rgb[:,:,0], loc2=2, loc1=3, edgecolor='white') pl.savefig("HI-CO-Dust_m0.35_oricloud-rezoom-trapezium_RGB.png", bbox_inches='tight', pad_inches=0) fn_trapezium_alma = 'Trapezium_GEMS_mosaic_redblueorange_ALMA_normed_small_contrast_bright.png' rgb = np.array([simple_norm(img_carmaori_12co_omc, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_carmaori_12co_omc), simple_norm(img_carmaori_13co_omc, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_carmaori_13co_omc), simple_norm(img_orion_bolocam_omc, min_percent=1, max_percent=99.5, log_a=5e2, stretch='log')(img_orion_bolocam_omc)]).T.swapaxes(0,1) ww = WCS(target_header_omc) plt.figure(figsize=(10,5), dpi=200, frameon=False) ax = plt.subplot(1,2,1, projection=ww,) ax.imshow(rgb) ax.axis('off') ax.axis([480, 1440, 480, 1440]) #ax.scatter(theta1c.ra, theta1c.dec, transform=ax.get_transform('world'), # marker='*', color='y') rgb = np.array(PIL.Image.open(fn_trapezium_alma))[::-1, :, :] ww = WCS(fits.Header.fromstring(trapezium_header, sep='\n')) axins = inset_axes(ax, width=4, height=4, bbox_to_anchor=[0.0, 0.0, 0.9, 1], bbox_transform=fig.transFigure, loc='center right', axes_class=astropy.visualization.wcsaxes.core.WCSAxes, axes_kwargs=dict(wcs=ww)) axins.imshow(rgb) axins.axis('off') #axins.scatter(theta1c.ra, theta1c.dec, transform=axins.get_transform('world'), # marker='*', color='y') mark_inset_generic(axins, ax, data=rgb[:,:,0], loc2=2, loc1=3, edgecolor='white') pl.savefig("HI-CO-Dust_m0.35_oricloud-rezoom-trapezium-ALMA_RGB.png", bbox_inches='tight', pad_inches=0) rgb = np.array([simple_norm(img_HI4PI_cmz, min_percent=0.01, max_percent=99.99, log_a=1e1, stretch='log')(img_HI4PI_cmz), simple_norm(img_CO21_cmz, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21_cmz), simple_norm(img_ThermalDust_cmz, min_percent=0.01, max_percent=99.99, log_a=5e2, stretch='log')(img_ThermalDust_cmz)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header_cmz),) #frame_class=EllipticalFrame) ax.imshow(rgb_scaled) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) ax.set_xlim(240*2, 720*2) ax.set_ylim(120*2, 360*2) ax.imshow(atlasgal[0].data, norm=simple_norm(atlasgal[0].data, min_percent=5, max_percent=99.9, stretch='asinh'), transform=ax.get_transform(WCS(atlasgal[0].header)), cmap=orange_transparent) cmzrect = regions.RectangleSkyRegion( coordinates.SkyCoord(0*u.deg, 0*u.deg, frame='galactic'), width=12*u.deg, height=4.8*u.deg, ) pcmzrect = cmzrect.to_pixel(ax.wcs) pcmzrect.visual['edgecolor'] = 'w' pcmzrect.plot(ax=ax) pl.savefig("HI4PI_CO_Dust_withATLASGALOrangetransparent_cmzbox.png", bbox_inches='tight') hnco_cube = SpectralCube.read("CMZ_3mm_HNCO.fits", use_dask=True) hnco_cube.shape #[Out]# (327, 165, 765) target_header_cmz_hires = fits.Header.fromstring(""" NAXIS = 2 NAXIS1 = 1200 NAXIS2 = 480 CTYPE1 = 'GLON-CAR' CRPIX1 = 600.5 CRVAL1 = 0.0 CDELT1 = -0.01 CUNIT1 = 'deg ' CTYPE2 = 'GLAT-CAR' CRPIX2 = 240.5 CRVAL2 = 0.0 CDELT2 = 0.01 CUNIT2 = 'deg ' COORDSYS= 'icrs ' """, sep='\n') hi_cube.spectral_slab(100*u.km/u.s, 300*u.km/u.s).mean(axis=0).quicklook() hi_cube.spectral_slab(-300*u.km/u.s, -100*u.km/u.s).mean(axis=0).quicklook() mean_hi_highv = (hi_cube.spectral_slab(-300*u.km/u.s, -100*u.km/u.s).mean(axis=0) + hi_cube.spectral_slab(100*u.km/u.s, 300*u.km/u.s).mean(axis=0)) fn_meerkat = 'MeerKAT_Galactic_Centre_1284MHz-StokesI.fits' if not os.path.exists(fn_meerkat): meerkat = download_file(f"https://archive-gw-1.kat.ac.za/public/repository/10.48479/fyst-hj47/data/{fn_meerkat}") shutil.move(meerkat, fn_meerkat) meerkat = fits.open(fn_meerkat) glimpsei4 = fits.open('/orange/adamginsburg/cmz/glimpse_data/GLM_00000+0000_mosaic_I4.fits') img_hi_cmz_hires, footprint = reproject.reproject_interp(mean_hi_highv.hdu, target_header_cmz_hires) img_agal_cmz_hires, footprint = reproject.reproject_interp(atlasgal, target_header_cmz_hires) img_meerkat_cmz_hires, footprint = reproject.reproject_interp((meerkat[0].data.squeeze(), WCS(meerkat[0].header).celestial), target_header_cmz_hires) img_co_cmz_hires, footprint = reproject.reproject_interp(mean_co.hdu, target_header_cmz_hires) for im in (img_hi_cmz_hires, img_agal_cmz_hires, img_co_cmz_hires): im[im==0] = np.nan pl.imshow(img_meerkat_cmz_hires, norm=simple_norm(img_meerkat_cmz_hires, min_percent=1, max_percent=99.99, log_a=5e-1, stretch='log')) #[Out]# pl.imshow(img_hi_cmz_hires, norm=simple_norm(img_hi_cmz_hires, min_percent=25, max_percent=99.99, log_a=5e-1, stretch='linear')) #[Out]# pl.imshow(img_co_cmz_hires, norm=simple_norm(img_co_cmz_hires, min_percent=5, max_percent=99.99, log_a=1e1, stretch='asinh')) #[Out]# pl.imshow(img_agal_cmz_hires, norm=simple_norm(img_agal_cmz_hires, min_percent=1, max_percent=99.9, log_a=3e1, stretch='log')) #[Out]# rgb = np.array([simple_norm(img_hi_cmz_hires, min_percent=25, max_percent=99.99, log_a=4e1, stretch='linear')(img_hi_cmz_hires), simple_norm(img_co_cmz_hires, min_percent=5, max_percent=99.99, log_a=4e1, stretch='asinh')(img_co_cmz_hires), simple_norm(img_agal_cmz_hires, min_percent=1, max_percent=99.99, log_a=4e1, stretch='log')(img_agal_cmz_hires)]).T.swapaxes(0,1) hsv = rgb_to_hsv(np.nan_to_num(rgb)) hsv[:,:,0] += -0.35#.45 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header_cmz_hires),) #frame_class=EllipticalFrame) ax.imshow(rgb_scaled) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) ax.set_xlim(100, 1100); ax.set_ylim(140, 340); pl.savefig("HI_CO_Dust_BigCMZ.png", bbox_inches='tight') rgb = np.array([simple_norm(img_HI4PI_cmz, min_percent=0.01, max_percent=99.99, log_a=1e1, stretch='log')(img_HI4PI_cmz), simple_norm(img_CO21_cmz, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21_cmz), simple_norm(img_ThermalDust_cmz, min_percent=0.01, max_percent=99.99, log_a=5e2, stretch='log')(img_ThermalDust_cmz)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) fig = plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(2,1,2, projection=WCS(target_header_cmz),) #frame_class=EllipticalFrame) ax.imshow(rgb_scaled) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) #ax.set_xlim(240*2, 720*2) ax.set_xlim(0, 1920) ax.set_ylim(120*2, 360*2) ax.imshow(atlasgal[0].data, norm=simple_norm(atlasgal[0].data, min_percent=5, max_percent=99.9, stretch='asinh'), transform=ax.get_transform(WCS(atlasgal[0].header)), cmap=orange_transparent) #cmzrect = regions.RectangleSkyRegion( # coordinates.SkyCoord(0*u.deg, 0*u.deg, frame='galactic'), # width=12*u.deg, # height=4.8*u.deg, #) #pcmzrect = cmzrect.to_pixel(ax.wcs) #pcmzrect.visual['edgecolor'] = 'w' #pcmzrect.plot(ax=ax) rgb = np.array([simple_norm(img_hi_cmz_hires, min_percent=25, max_percent=99.99, log_a=4e1, stretch='linear')(np.nan_to_num(img_hi_cmz_hires)), simple_norm(img_co_cmz_hires, min_percent=5, max_percent=99.99, log_a=4e1, stretch='asinh')(np.nan_to_num(img_co_cmz_hires)), simple_norm(img_agal_cmz_hires, min_percent=1, max_percent=99.99, log_a=4e1, stretch='log')(np.nan_to_num(img_agal_cmz_hires))]).T.swapaxes(0,1) hsv = rgb_to_hsv(np.nan_to_num(rgb)) hsv[:,:,0] += -0.35#.45 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) axins = inset_axes(ax, width=4, height=4, bbox_to_anchor=[0, 0.25, 1, 0.5], bbox_transform=fig.transFigure, loc='lower center', axes_class=astropy.visualization.wcsaxes.core.WCSAxes, axes_kwargs=dict(wcs=WCS(target_header_cmz_hires))) axins.imshow(rgb_scaled) axins.coords['glat'].set_ticklabel(visible=False) axins.coords['glon'].set_ticklabel(visible=False) axins.coords['glat'].set_ticks_visible(False) axins.coords['glon'].set_ticks_visible(False) mark_inset_generic(axins, ax, data=img_hi_cmz_hires, loc1=4, edgecolor='w') pl.savefig("HI4PI_CO_Dust_withATLASGALOrangetransparent_cmzbox_insetzoom.png", bbox_inches='tight') import PIL import pyavm np.array(PIL.Image.open('gc_fullres_6.jpg')).shape #[Out]# (3295, 6473, 3) rgb = np.array([simple_norm(img_hi_cmz_hires, min_percent=25, max_percent=99.99, log_a=4e1, stretch='linear')(np.nan_to_num(img_hi_cmz_hires)), simple_norm(img_co_cmz_hires, min_percent=5, max_percent=99.99, log_a=4e1, stretch='asinh')(np.nan_to_num(img_co_cmz_hires)), simple_norm(img_agal_cmz_hires, min_percent=1, max_percent=99.99, log_a=4e1, stretch='log')(np.nan_to_num(img_agal_cmz_hires))]).T.swapaxes(0,1) hsv = rgb_to_hsv(np.nan_to_num(rgb)) hsv[:,:,0] += -0.35 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) fig = plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(2,1,1, projection=WCS(target_header_cmz_hires),) ax.imshow(rgb_scaled) ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) rgb = np.array(PIL.Image.open('gc_fullres_6.jpg')) ww = WCS(fits.Header.fromtextfile('gc_fullres_6.wcs')) axins = inset_axes(ax, width=4, height=4, bbox_to_anchor=[0.01, 0.2, 1, 0.5], bbox_transform=fig.transFigure, loc='upper center', axes_class=astropy.visualization.wcsaxes.core.WCSAxes, axes_kwargs=dict(wcs=ww)) axins.imshow(rgb[::-1,:,:]) axins.coords['glat'].set_ticklabel(visible=False) axins.coords['glon'].set_ticklabel(visible=False) axins.coords['glat'].set_ticks_visible(False) axins.coords['glon'].set_ticks_visible(False) mark_inset_generic(axins, ax, data=rgb[:,:,0], loc2=2, edgecolor='w') pl.savefig("HI-CO-Dust_zoominto3color.png", bbox_inches='tight') rgb = np.array(PIL.Image.open('gc_fullres_6.jpg'))[::-1,:,:] ww = WCS(fits.Header.fromtextfile('gc_fullres_6.wcs')) fig = plt.figure(figsize=(10,5), dpi=200, frameon=False) ax = plt.subplot(1,1,1, projection=ww) ax.imshow(rgb) ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) aces7m = fits.open('/orange/adamginsburg/ACES/mosaics/7m_continuum_mosaic.fits') aces7mwcs = WCS(aces7m[0].header) aces7mrepr,_ = reproject.reproject_interp(aces7m, ww, shape_out=rgb.shape[:2]) aces7mrepr[aces7mrepr < 0] = np.nan ax.imshow(aces7mrepr, norm=simple_norm(aces7mrepr, stretch='log', max_percent=99.96, min_percent=1), cmap=orange_transparent) pl.savefig("gc_fullres_6_small_withACES7m.png", bbox_inches='tight', pad_inches=0) rgb = np.array(PIL.Image.open('gc_fullres_6.jpg'))[::-1,:,:] ww = WCS(fits.Header.fromtextfile('gc_fullres_6.wcs')) fig = plt.figure(figsize=(10,5), dpi=200, frameon=False) ax = plt.subplot(2,1,2, projection=ww) ax.imshow(rgb) ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) aces7m = fits.open('/orange/adamginsburg/ACES/mosaics/7m_continuum_mosaic.fits') aces7mwcs = WCS(aces7m[0].header) aces7mrepr,_ = reproject.reproject_interp(aces7m, ww, shape_out=rgb.shape[:2]) aces7mrepr[aces7mrepr < 0] = np.nan ax.imshow(aces7mrepr, norm=simple_norm(aces7mrepr, stretch='log', max_percent=99.96, min_percent=1), cmap=orange_transparent) aces12m = fits.open('/orange/adamginsburg/ACES/mosaics/12m_continuum_mosaic.fits') aces12mwcs = WCS(aces12m[0].header) aces12m[0].data[aces12m[0].data < -0.0005] = 0 axins = inset_axes(ax, width=4, height=4, bbox_to_anchor=[0.01, 0.5, 1, 0.5], bbox_transform=fig.transFigure, loc='upper center', axes_class=astropy.visualization.wcsaxes.core.WCSAxes, axes_kwargs=dict(wcs=aces12mwcs)) axins.imshow(aces12m[0].data, norm=simple_norm(aces12m[0].data, stretch='log', min_percent=None, max_percent=None, min_cut=-0.0005, max_cut=0.05,), cmap='gray', ) axins.coords['glat'].set_ticklabel(visible=False) axins.coords['glon'].set_ticklabel(visible=False) axins.coords['glat'].set_ticks_visible(False) axins.coords['glon'].set_ticks_visible(False) mark_inset_generic(axins, ax, data=aces12m[0].data, loc1=4, edgecolor='y') pl.savefig("gc_fullres_6_small_withACES7m_12minset.png", bbox_inches='tight') import regions from astropy import coordinates from astropy import units as u, constants import astropy.visualization.wcsaxes from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, inset_axes, mark_inset aces12m = fits.open('/orange/adamginsburg/ACES/mosaics/12m_continuum_mosaic.fits') aces12mwcs = WCS(aces12m[0].header) aces12m[0].data[aces12m[0].data < -0.0005] = 0 fig = plt.figure(figsize=(10,5), frameon=False) ax = plt.subplot(1,1,1, projection=aces12mwcs) ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) ax.axis('off') sgrb2m = regions.CircleSkyRegion( coordinates.SkyCoord(000.6672*u.deg, -00.03640*u.deg, frame='galactic'), radius=25*u.arcsec) sgrb2n = regions.CircleSkyRegion( coordinates.SkyCoord(000.6773*u.deg, -00.029*u.deg, frame='galactic'), radius=25*u.arcsec) sgrb2m.visual['edgecolor'] = 'r' sgrb2n.visual['edgecolor'] = 'r' merge = psgrb2n & psgrb2m merge_mask = merge.to_mask() slcs_merge,_ = merge_mask.get_overlap_slices(aces12m[0].data.shape) # moved down here for debugging... ax.imshow(aces12m[0].data, norm=simple_norm(aces12m[0].data, stretch='log', min_percent=None, max_percent=None, min_cut=-0.0005, max_cut=0.05,), cmap='gray', ) psgrb2m = sgrb2m.to_pixel(ax.wcs) point_sgrb2m = psgrb2m.plot(ax=ax) msgrb2m = psgrb2m.to_mask() slcs_sgrb2m,_ = msgrb2m.get_overlap_slices(aces12m[0].data.shape) bbox = [0.01, 0, 0.2, 1.05] axins1 = zoomed_inset_axes(ax, 30, loc='upper left', bbox_to_anchor=bbox, bbox_transform=fig.transFigure, axes_class=astropy.visualization.wcsaxes.core.WCSAxes, axes_kwargs=dict(wcs=ax.wcs) ) axins1.axis('off') axins1.imshow(aces12m[0].data, norm=simple_norm(aces12m[0].data, stretch='log', min_percent=None, max_percent=None, min_cut=-0.0005, max_cut=0.5,), cmap='gray' ) axins1.axis(msgrb2m.bbox.extent) mark1 = mark_inset(ax, axins1, loc1=1, loc2=3, edgecolor='r') point_sgrb2m_in = sgrb2m.to_pixel(axins1.wcs).plot(ax=axins1) axins1.coords['glat'].set_ticklabel(visible=False) axins1.coords['glon'].set_ticklabel(visible=False) axins1.coords['glat'].set_ticks_visible(False) axins1.coords['glon'].set_ticks_visible(False) pl.savefig("aces_mark_sgrb2m.png", bbox_inches='tight', dpi=200) point_sgrb2m.set_visible(False) point_sgrb2m_in.set_visible(False) #axins1.set_visible(False) for pp in mark1: pp.set_visible(False) psgrb2n = sgrb2n.to_pixel(ax.wcs) psgrb2n.visual['edgecolor'] = 'r' point_sgrb2n = psgrb2n.plot(ax=ax) msgrb2n = psgrb2n.to_mask() slcs_sgrb2n,_ = msgrb2n.get_overlap_slices(aces12m[0].data.shape) #axins2 = zoomed_inset_axes(ax, 30, loc='upper left', bbox_to_anchor=bbox, bbox_transform=fig.transFigure, # axes_class=astropy.visualization.wcsaxes.core.WCSAxes, # axes_kwargs=dict(wcs=ax.wcs)) axins1.imshow(aces12m[0].data, norm=simple_norm(aces12m[0].data, stretch='log', min_percent=None, max_percent=None, min_cut=-0.0005, max_cut=0.5,), cmap='gray' ) axins1.axis(msgrb2n.bbox.extent) mark2 = mark_inset(ax, axins1, loc1=1, loc2=3, edgecolor='r') point_sgrb2n_in = sgrb2n.to_pixel(axins1.wcs).plot(ax=axins1) #axins2.coords['glat'].set_ticklabel(visible=False) #axins2.coords['glon'].set_ticklabel(visible=False) #axins2.coords['glat'].set_ticks_visible(False) #axins2.coords['glon'].set_ticks_visible(False) pl.savefig("aces_mark_sgrb2n.png", bbox_inches='tight', dpi=200) point_sgrb2n.set_visible(False) point_sgrb2n_in.set_visible(False) for pp in mark2: pp.set_visible(False) axins1.set_visible(False) axins3 = zoomed_inset_axes(ax, 15, loc='upper left', bbox_to_anchor=[0.05, 0, 0.2, 1], bbox_transform=fig.transFigure, axes_class=astropy.visualization.wcsaxes.core.WCSAxes, axes_kwargs=dict(wcs=ax.wcs)) axins3.axis('off') axins3.imshow(aces12m[0].data, norm=simple_norm(aces12m[0].data, stretch='log', min_percent=None, max_percent=None, min_cut=-0.0005, max_cut=0.5,), cmap='gray' ) axins3.axis(merge_mask.bbox.extent) mark3 = mark_inset(ax, axins3, loc1=1, loc2=3, edgecolor='r') point_sgrb2n_in = sgrb2n.to_pixel(axins3.wcs).plot(ax=axins3) point_sgrb2m_in = sgrb2m.to_pixel(axins3.wcs).plot(ax=axins3) axins3.coords['glat'].set_ticklabel(visible=False) axins3.coords['glon'].set_ticklabel(visible=False) axins3.coords['glat'].set_ticks_visible(False) axins3.coords['glon'].set_ticks_visible(False) pl.savefig("aces_mark_sgrb2mANDn.png", bbox_inches='tight', dpi=200) import regions from astropy import coordinates from astropy import units as u, constants import astropy.visualization.wcsaxes from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, inset_axes, mark_inset aces12m = fits.open('/orange/adamginsburg/ACES/mosaics/12m_continuum_mosaic.fits') aces12mwcs = WCS(aces12m[0].header) aces12m[0].data[aces12m[0].data < -0.0005] = 0 fig = plt.figure(figsize=(10,5), frameon=False) ax = plt.subplot(1,1,1, projection=aces12mwcs) ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) ax.axis('off') sgrb2m = regions.CircleSkyRegion( coordinates.SkyCoord(000.6672*u.deg, -00.03640*u.deg, frame='galactic'), radius=25*u.arcsec) sgrb2n = regions.CircleSkyRegion( coordinates.SkyCoord(000.6773*u.deg, -00.029*u.deg, frame='galactic'), radius=25*u.arcsec) sgrb2m.visual['edgecolor'] = 'r' sgrb2n.visual['edgecolor'] = 'r' psgrb2m = sgrb2m.to_pixel(ax.wcs) point_sgrb2m = psgrb2m.plot(ax=ax) msgrb2m = psgrb2m.to_mask() slcs_sgrb2m,_ = msgrb2m.get_overlap_slices(aces12m[0].data.shape) psgrb2n = sgrb2n.to_pixel(ax.wcs) psgrb2n.visual['edgecolor'] = 'r' point_sgrb2n = psgrb2n.plot(ax=ax) msgrb2n = psgrb2n.to_mask() slcs_sgrb2n,_ = msgrb2n.get_overlap_slices(aces12m[0].data.shape) merge = psgrb2n & psgrb2m merge_mask = merge.to_mask() slcs_merge,_ = merge_mask.get_overlap_slices(aces12m[0].data.shape) # moved down here for debugging... ax.imshow(aces12m[0].data, norm=simple_norm(aces12m[0].data, stretch='log', min_percent=None, max_percent=None, min_cut=-0.0005, max_cut=0.05,), cmap='gray', ) bbox = [0.01, 0, 0.2, 1.05] axins1 = zoomed_inset_axes(ax, 30, loc='upper left', bbox_to_anchor=bbox, bbox_transform=fig.transFigure, axes_class=astropy.visualization.wcsaxes.core.WCSAxes, axes_kwargs=dict(wcs=ax.wcs) ) axins1.axis('off') axins1.imshow(aces12m[0].data, norm=simple_norm(aces12m[0].data, stretch='log', min_percent=None, max_percent=None, min_cut=-0.0005, max_cut=0.5,), cmap='gray' ) axins1.axis(msgrb2m.bbox.extent) mark1 = mark_inset(ax, axins1, loc1=1, loc2=3, edgecolor='r') point_sgrb2m_in = sgrb2m.to_pixel(axins1.wcs).plot(ax=axins1) axins1.coords['glat'].set_ticklabel(visible=False) axins1.coords['glon'].set_ticklabel(visible=False) axins1.coords['glat'].set_ticks_visible(False) axins1.coords['glon'].set_ticks_visible(False) pl.savefig("aces_mark_sgrb2m.png", bbox_inches='tight', dpi=200) point_sgrb2m.set_visible(False) point_sgrb2m_in.set_visible(False) #axins1.set_visible(False) for pp in mark1: pp.set_visible(False) #axins2 = zoomed_inset_axes(ax, 30, loc='upper left', bbox_to_anchor=bbox, bbox_transform=fig.transFigure, # axes_class=astropy.visualization.wcsaxes.core.WCSAxes, # axes_kwargs=dict(wcs=ax.wcs)) axins1.imshow(aces12m[0].data, norm=simple_norm(aces12m[0].data, stretch='log', min_percent=None, max_percent=None, min_cut=-0.0005, max_cut=0.5,), cmap='gray' ) axins1.axis(msgrb2n.bbox.extent) mark2 = mark_inset(ax, axins1, loc1=1, loc2=3, edgecolor='r') point_sgrb2n_in = sgrb2n.to_pixel(axins1.wcs).plot(ax=axins1) #axins2.coords['glat'].set_ticklabel(visible=False) #axins2.coords['glon'].set_ticklabel(visible=False) #axins2.coords['glat'].set_ticks_visible(False) #axins2.coords['glon'].set_ticks_visible(False) pl.savefig("aces_mark_sgrb2n.png", bbox_inches='tight', dpi=200) point_sgrb2n.set_visible(False) point_sgrb2n_in.set_visible(False) for pp in mark2: pp.set_visible(False) axins1.set_visible(False) axins3 = zoomed_inset_axes(ax, 15, loc='upper left', bbox_to_anchor=[0.05, 0, 0.2, 1], bbox_transform=fig.transFigure, axes_class=astropy.visualization.wcsaxes.core.WCSAxes, axes_kwargs=dict(wcs=ax.wcs)) axins3.axis('off') axins3.imshow(aces12m[0].data, norm=simple_norm(aces12m[0].data, stretch='log', min_percent=None, max_percent=None, min_cut=-0.0005, max_cut=0.5,), cmap='gray' ) axins3.axis(merge_mask.bbox.extent) mark3 = mark_inset(ax, axins3, loc1=1, loc2=3, edgecolor='r') point_sgrb2n_in = sgrb2n.to_pixel(axins3.wcs).plot(ax=axins3) point_sgrb2m_in = sgrb2m.to_pixel(axins3.wcs).plot(ax=axins3) axins3.coords['glat'].set_ticklabel(visible=False) axins3.coords['glon'].set_ticklabel(visible=False) axins3.coords['glat'].set_ticks_visible(False) axins3.coords['glon'].set_ticks_visible(False) pl.savefig("aces_mark_sgrb2mANDn.png", bbox_inches='tight', dpi=200) # FAILED attempt to show with 12m rgb = np.array(PIL.Image.open('gc_fullres_6.jpg'))[::-1,:,:] ww = WCS(fits.Header.fromtextfile('gc_fullres_6.wcs')) fig = plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=ww) ax.imshow(rgb) ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) aces12m = fits.open('/orange/adamginsburg/ACES/mosaics/12m_continuum_mosaic.fits') aces12m[0].data[aces12m[0].data < -0.0005] = 0 aces12mwcs = WCS(aces12m[0].header) aces12mrepr,_ = reproject.reproject_interp(aces12m, ww, shape_out=rgb.shape[:2]) aces12mrepr[aces12mrepr < -5e-4] = 0 aces12mrepr[aces12mrepr < 1e-4] = np.nan ax.imshow(aces12mrepr, norm=simple_norm(aces12mrepr, stretch='log', max_percent=99.5, min_cut=1e-4), cmap=orange_transparent) pl.savefig("gc_fullres_6_small_withACES12m.png", bbox_inches='tight') rgb = np.array([simple_norm(img_HI4PI_cmz, min_percent=0.01, max_percent=99.99, log_a=1e1, stretch='log')(img_HI4PI_cmz), simple_norm(img_CO21_cmz, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21_cmz), simple_norm(img_ThermalDust_cmz, min_percent=0.01, max_percent=99.99, log_a=5e2, stretch='log')(img_ThermalDust_cmz)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) plt.figure(figsize=(10,5), dpi=200) ax = plt.subplot(1,1,1, projection=WCS(target_header_cmz),) #frame_class=EllipticalFrame) ax.imshow(rgb_scaled) #ax.coords.grid(color='white') ax.coords['glat'].set_ticklabel(visible=False) ax.coords['glon'].set_ticklabel(visible=False) ax.coords['glat'].set_ticks_visible(False) ax.coords['glon'].set_ticks_visible(False) ax.set_xlim(240*2, 720*2) ax.set_ylim(120*2, 360*2) ax.imshow(meerkat[0].data.squeeze(), norm=simple_norm(meerkat[0].data, min_percent=0.1, max_percent=99.9, stretch='log'), transform=ax.get_transform(WCS(meerkat[0].header).celestial), interpolation='none', cmap='gray') #[Out]# aces12m = fits.open('/orange/adamginsburg/ACES/mosaics/12m_continuum_mosaic.fits') aces12mwcs = WCS(aces12m[0].header) if os.path.exists('meerkat_projectedto_aces.fits'): img_meerkat_ACES = fits.getdata('meerkat_projectedto_aces.fits') else: img_meerkat_ACES,_ = reproject.reproject_interp((meerkat[0].data.squeeze(), WCS(meerkat[0].header).celestial), aces12mwcs.celestial, shape_out=aces12m[0].data.shape) fits.PrimaryHDU(data=img_meerkat_ACES, header=aces12m[0].header).writeto('meerkat_projectedto_aces.fits') if os.path.exists('glimpsei4_projectedto_aces.fits'): img_glimpsei4_ACES = fits.getdata('glimpsei4_projectedto_aces.fits') else: img_glimpsei4_ACES,_ = reproject.reproject_interp(glimpsei4, aces12mwcs.celestial, shape_out=aces12m[0].data.shape) fits.PrimaryHDU(data=img_glimpsei4_ACES, header=aces12m[0].header).writeto('glimpsei4_projectedto_aces.fits') rgb = np.array([simple_norm(img_meerkat_ACES, min_percent=1, max_percent=99.95, log_a=4e1, stretch='log')(np.nan_to_num(img_meerkat_ACES)), simple_norm(aces12m[0].data.squeeze(), min_percent=1, max_percent=99.95, log_a=5e1, stretch='log')(np.nan_to_num(aces12m[0].data.squeeze())), simple_norm(img_glimpsei4_ACES, min_percent=1, max_percent=99.9, log_a=4e1, stretch='log')(np.nan_to_num(img_glimpsei4_ACES))]).T.swapaxes(0,1) pl.imshow(rgb[4000:6000,12000:14000,:], origin='lower') #[Out]# pl.imshow(rgb, origin='lower') #[Out]# img_meerkat_ACES.shape, img_glimpsei4_ACES.shape, aces12m[0].data.shape #[Out]# ((10200, 27800), (10200, 27800), (10200, 27800)) target_header_innercmz = fits.Header.fromstring(""" NAXIS = 2 NAXIS1 = 3000 NAXIS2 = 1000 CTYPE1 = 'GLON-CAR' CRPIX1 = 1500.5 CRVAL1 = 0.0 CDELT1 = -0.001 CUNIT1 = 'deg ' CTYPE2 = 'GLAT-CAR' CRPIX2 = 500.5 CRVAL2 = 0.0 CDELT2 = 0.001 CUNIT2 = 'deg ' COORDSYS= 'icrs ' """, sep='\n') # attempt progressive zoom import numpy as np import matplotlib.animation as animation rgb = np.array([simple_norm(img_HI4PI_cmz, min_percent=0.01, max_percent=99.99, log_a=1e1, stretch='log')(img_HI4PI_cmz), simple_norm(img_CO21_cmz, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21_cmz), simple_norm(img_ThermalDust_cmz, min_percent=0.01, max_percent=99.99, log_a=5e2, stretch='log')(img_ThermalDust_cmz)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) fig = pl.figure(figsize=(10,5), dpi=200, frameon=False) ax = pl.subplot(1,1,1, projection=WCS(target_header_cmz),) #frame_class=EllipticalFrame) ax.imshow(rgb_scaled) ax.axis('off') #ax.set_ylim(120*2, 360*2) axlims = ax.axis() cx = np.mean(ax.get_xlim()) cy = np.mean(ax.get_ylim()) dy0 = np.abs(np.diff(ax.get_ylim())) dx0 = np.abs(np.diff(ax.get_xlim())) nframes = 90 maxzoom = 5 def animate(n, nframes=nframes): # lol at math. # n = 1 -> denom = 1 # n = 63 -> denom = 5 # y = m x + b # 5 = m 63 + b # 1 = m 1 + b # 5 = m 63 + 1 - m # m = 4/62 # b = 58/62 denom = (maxzoom-1)/(nframes-1) * n + (nframes - maxzoom)/(nframes - 1) dy = dy0 * 1/denom dx = dx0 * 1/denom ax.axis((cx-dx, cx+dx, cy-dy, cy+dy)) print('.', end='') return fig, anim = animation.FuncAnimation(fig, animate, frames=nframes, repeat_delay=5000, interval=50) anim.save('zoom_anim_HI-CO-Dust_m0.35.gif') rgb = np.array([simple_norm(img_HI4PI_innergal, min_percent=0.01, max_percent=99.99, log_a=1e1, stretch='log')(img_HI4PI_innergal), simple_norm(img_CO21_innergal, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21_innergal), simple_norm(img_ThermalDust_innergal, min_percent=0.01, max_percent=99.99, log_a=5e2, stretch='log')(img_ThermalDust_innergal)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) rgb_scaled[rgb_scaled > 1] = 1 rgb_scaled[rgb_scaled < 0] = 0 rgb_full = np.array([simple_norm(img_HI4PI, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PI), simple_norm(img_CO21, min_percent=0.01, max_percent=99.90, log_a=2e1, stretch='log')(img_CO21), simple_norm(img_ThermalDust, min_percent=0.01, max_percent=99.90, log_a=5e2, stretch='log')(img_ThermalDust)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb_full) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_full_scaled = hsv_to_rgb(hsv) rgb_full_scaled[rgb_full_scaled > 1] = 1 rgb_full_scaled[rgb_full_scaled < 0] = 0 fig = pl.figure(figsize=(10,5), dpi=200, frameon=False) ax = pl.subplot(1,1,1, projection=WCS(target_header), frame_class=EllipticalFrame) ax.imshow(rgb_full_scaled) ax.axis('off') nframes = 80 maxzoom = 5 zoomfac = np.hstack([np.linspace(1, 5, nframes//2), np.linspace(1, 5, nframes//2)]) zoomfac = np.linspace(1, maxzoom, nframes) dy0 = np.abs(np.diff(ax.get_ylim()))/2 dx0 = np.abs(np.diff(ax.get_xlim()))/2 dxdy = [dx0, dy0] cx = np.mean(ax.get_xlim()) cy = np.mean(ax.get_ylim()) def animate(n, nframes=nframes): dx0, dy0 = dxdy #if n < nframes // 2: # ax = fig.gca() # cx = np.mean(ax.get_xlim()) # cy = np.mean(ax.get_ylim()) #else: # fig.clf() # ax = pl.subplot(1,1,1, projection=WCS(target_header_innergal),) # ax.imshow(rgb_scaled) # ax.axis('off') # cx = rgb_scaled.shape[1]/2 # cy = rgb_scaled.shape[0]/2 # #if n == nframes // 2: # dy0 = rgb_scaled.shape[0] # dx0 = rgb_scaled.shape[1] # dxdy[:] = (dx0, dy0) # print(f"\ndx0,dy0: {dxdy}") if n == nframes // 2: ax.imshow(atlasgal[0].data, norm=simple_norm(atlasgal[0].data, min_percent=5, max_percent=99.9, stretch='asinh'), transform=ax.get_transform(WCS(atlasgal[0].header)), cmap=orange_transparent) dy = dy0 / zoomfac[n] dx = dx0 / zoomfac[n] ax.axis((cx-dx, cx+dx, cy-dy, cy+dy)) print(f'{n}:{zoomfac[n]:0.1f}|', end='') return fig, anim = animation.FuncAnimation(fig, animate, frames=nframes, repeat_delay=5000, interval=50) anim.save('zoom_anim_molltorect_HI-CO-Dust_m0.35.gif') rgb_full_scaled.shape #[Out]# (960, 1920, 3) 360/1920 #[Out]# 0.1875 1920/2 #[Out]# 960.0 0.5 / 0.1875 #[Out]# 2.6666666666666665 np.linspace(1, 1/24., 120), 1/16. #[Out]# (array([1. , 0.99194678, 0.98389356, 0.97584034, 0.96778711, #[Out]# 0.95973389, 0.95168067, 0.94362745, 0.93557423, 0.92752101, #[Out]# 0.91946779, 0.91141457, 0.90336134, 0.89530812, 0.8872549 , #[Out]# 0.87920168, 0.87114846, 0.86309524, 0.85504202, 0.8469888 , #[Out]# 0.83893557, 0.83088235, 0.82282913, 0.81477591, 0.80672269, #[Out]# 0.79866947, 0.79061625, 0.78256303, 0.7745098 , 0.76645658, #[Out]# 0.75840336, 0.75035014, 0.74229692, 0.7342437 , 0.72619048, #[Out]# 0.71813725, 0.71008403, 0.70203081, 0.69397759, 0.68592437, #[Out]# 0.67787115, 0.66981793, 0.66176471, 0.65371148, 0.64565826, #[Out]# 0.63760504, 0.62955182, 0.6214986 , 0.61344538, 0.60539216, #[Out]# 0.59733894, 0.58928571, 0.58123249, 0.57317927, 0.56512605, #[Out]# 0.55707283, 0.54901961, 0.54096639, 0.53291317, 0.52485994, #[Out]# 0.51680672, 0.5087535 , 0.50070028, 0.49264706, 0.48459384, #[Out]# 0.47654062, 0.46848739, 0.46043417, 0.45238095, 0.44432773, #[Out]# 0.43627451, 0.42822129, 0.42016807, 0.41211485, 0.40406162, #[Out]# 0.3960084 , 0.38795518, 0.37990196, 0.37184874, 0.36379552, #[Out]# 0.3557423 , 0.34768908, 0.33963585, 0.33158263, 0.32352941, #[Out]# 0.31547619, 0.30742297, 0.29936975, 0.29131653, 0.28326331, #[Out]# 0.27521008, 0.26715686, 0.25910364, 0.25105042, 0.2429972 , #[Out]# 0.23494398, 0.22689076, 0.21883754, 0.21078431, 0.20273109, #[Out]# 0.19467787, 0.18662465, 0.17857143, 0.17051821, 0.16246499, #[Out]# 0.15441176, 0.14635854, 0.13830532, 0.1302521 , 0.12219888, #[Out]# 0.11414566, 0.10609244, 0.09803922, 0.08998599, 0.08193277, #[Out]# 0.07387955, 0.06582633, 0.05777311, 0.04971989, 0.04166667]), #[Out]# 0.0625) rgb = np.array([simple_norm(img_HI4PI_innergal, min_percent=0.01, max_percent=99.99, log_a=1e1, stretch='log')(img_HI4PI_innergal), simple_norm(img_CO21_innergal, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21_innergal), simple_norm(img_ThermalDust_innergal, min_percent=0.01, max_percent=99.99, log_a=5e2, stretch='log')(img_ThermalDust_innergal)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) rgb_scaled[rgb_scaled > 1] = 1 rgb_scaled[rgb_scaled < 0] = 0 rgb_full = np.array([simple_norm(img_HI4PI, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PI), simple_norm(img_CO21, min_percent=0.01, max_percent=99.90, log_a=2e1, stretch='log')(img_CO21), simple_norm(img_ThermalDust, min_percent=0.01, max_percent=99.90, log_a=5e2, stretch='log')(img_ThermalDust)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb_full) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_full_scaled = hsv_to_rgb(hsv) rgb_full_scaled[rgb_full_scaled > 1] = 1 rgb_full_scaled[rgb_full_scaled < 0] = 0 fig = pl.figure(figsize=(10,5), dpi=200, frameon=False) ax = pl.subplot(1,1,1, projection=WCS(target_header), frame_class=EllipticalFrame) ax.imshow(rgb_full_scaled) ax.axis('off') nframes = 180 zoomfac = np.geomspace(1, 1/180., nframes) dy0 = np.abs(np.diff(ax.get_ylim()))/2 dx0 = np.abs(np.diff(ax.get_xlim()))/2 dxdy = [dx0, dy0] cx = np.mean(ax.get_xlim()) - 2.6 # zoom in between Sgr A and Sgr B2 cy = np.mean(ax.get_ylim()) rgbcmz = np.array(PIL.Image.open('gc_fullres_6.jpg'))[::-1,:,:] wwcmz = WCS(fits.Header.fromtextfile('gc_fullres_6.wcs')) def animate(n, nframes=nframes): dx0, dy0 = dxdy if n == 40: ax.imshow(atlasgal[0].data, norm=simple_norm(atlasgal[0].data, min_percent=5, max_percent=99.9, stretch='asinh'), transform=ax.get_transform(WCS(atlasgal[0].header)), cmap=orange_transparent) if n == 80: ax.imshow(rgbcmz, transform=ax.get_transform(wwcmz)) dy = dy0 * zoomfac[n] dx = dx0 * zoomfac[n] ax.axis((cx-dx, cx+dx, cy-dy, cy+dy)) print(f'{n}:{zoomfac[n]:0.1f}|', end='') return fig, anim = animation.FuncAnimation(fig, animate, frames=nframes, repeat_delay=5000, interval=50) anim.save('zoom_anim_cmz_linear_HI-CO-Dust_m0.35.gif') rgb = np.array([simple_norm(img_HI4PI_innergal, min_percent=0.01, max_percent=99.99, log_a=1e1, stretch='log')(img_HI4PI_innergal), simple_norm(img_CO21_innergal, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21_innergal), simple_norm(img_ThermalDust_innergal, min_percent=0.01, max_percent=99.99, log_a=5e2, stretch='log')(img_ThermalDust_innergal)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) rgb_scaled[rgb_scaled > 1] = 1 rgb_scaled[rgb_scaled < 0] = 0 rgb_full = np.array([simple_norm(img_HI4PI, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PI), simple_norm(img_CO21, min_percent=0.01, max_percent=99.90, log_a=2e1, stretch='log')(img_CO21), simple_norm(img_ThermalDust, min_percent=0.01, max_percent=99.90, log_a=5e2, stretch='log')(img_ThermalDust)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb_full) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_full_scaled = hsv_to_rgb(hsv) rgb_full_scaled[rgb_full_scaled > 1] = 1 rgb_full_scaled[rgb_full_scaled < 0] = 0 fig = pl.figure(figsize=(10,5), dpi=200, frameon=False) ax = pl.subplot(1,1,1, projection=WCS(target_header), frame_class=EllipticalFrame) ax.imshow(rgb_full_scaled) ax.axis('off') nframes = 80 maxzoom = 5 zoomfac = np.linspace(1, 1/maxzoom, nframes) zoomfac = np.geomspace(1, 1/maxzoom, nframes) dy0 = np.abs(np.diff(ax.get_ylim()))/2 dx0 = np.abs(np.diff(ax.get_xlim()))/2 dxdy = [dx0, dy0] cx = np.mean(ax.get_xlim()) cy = np.mean(ax.get_ylim()) def animate(n, nframes=nframes): dx0, dy0 = dxdy #if n < nframes // 2: # ax = fig.gca() # cx = np.mean(ax.get_xlim()) # cy = np.mean(ax.get_ylim()) #else: # fig.clf() # ax = pl.subplot(1,1,1, projection=WCS(target_header_innergal),) # ax.imshow(rgb_scaled) # ax.axis('off') # cx = rgb_scaled.shape[1]/2 # cy = rgb_scaled.shape[0]/2 # #if n == nframes // 2: # dy0 = rgb_scaled.shape[0] # dx0 = rgb_scaled.shape[1] # dxdy[:] = (dx0, dy0) # print(f"\ndx0,dy0: {dxdy}") if n == nframes // 2: ax.imshow(atlasgal[0].data, norm=simple_norm(atlasgal[0].data, min_percent=5, max_percent=99.9, stretch='asinh'), transform=ax.get_transform(WCS(atlasgal[0].header)), cmap=orange_transparent) dy = dy0 * zoomfac[n] dx = dx0 * zoomfac[n] ax.axis((cx-dx, cx+dx, cy-dy, cy+dy)) print(f'{n}:{zoomfac[n]:0.1f}|', end='') return fig, anim = animation.FuncAnimation(fig, animate, frames=nframes, repeat_delay=5000, interval=50) anim.save('zoom_anim_molltorect_linear_HI-CO-Dust_m0.35.gif') # convert zoom_anim_molltorect_HI-CO-Dust_m0.35.gif \( +clone -set delay 500 \) +swap +delete zoom_anim_molltorect_HI-CO-Dust_m0.35_with_pause.gif sgrc = SkyCoord.from_name('Sgr C') sgrc #[Out]# ax.wcs.world_to_pixel(sgrb2), ax.wcs.world_to_pixel(sgra), ax.wcs.world_to_pixel(sgrc) rgb = np.array([simple_norm(img_HI4PI_innergal, min_percent=0.01, max_percent=99.99, log_a=1e1, stretch='log')(img_HI4PI_innergal), simple_norm(img_CO21_innergal, min_percent=0.01, max_percent=99.99, log_a=3e1, stretch='log')(img_CO21_innergal), simple_norm(img_ThermalDust_innergal, min_percent=0.01, max_percent=99.99, log_a=5e2, stretch='log')(img_ThermalDust_innergal)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_scaled = hsv_to_rgb(hsv) rgb_scaled[rgb_scaled > 1] = 1 rgb_scaled[rgb_scaled < 0] = 0 rgb_full = np.array([simple_norm(img_HI4PI, min_percent=0.01, max_percent=99.99, log_a=2e1, stretch='log')(img_HI4PI), simple_norm(img_CO21, min_percent=0.01, max_percent=99.90, log_a=2e1, stretch='log')(img_CO21), simple_norm(img_ThermalDust, min_percent=0.01, max_percent=99.90, log_a=5e2, stretch='log')(img_ThermalDust)]).T.swapaxes(0,1) hsv = rgb_to_hsv(rgb_full) hsv[:,:,0] += -0.35 # 0.25 = 90/360 hsv[:,:,0] = hsv[:,:,0] % 1 rgb_full_scaled = hsv_to_rgb(hsv) rgb_full_scaled[rgb_full_scaled > 1] = 1 rgb_full_scaled[rgb_full_scaled < 0] = 0 fig = pl.figure(figsize=(10,5), dpi=200, frameon=False) ax = pl.subplot(1,1,1, projection=WCS(target_header), frame_class=EllipticalFrame) ax.imshow(rgb_full_scaled) ax.axis('off') nframes = 180 zoomfac = np.geomspace(1, 1/180., nframes) dy0 = np.abs(np.diff(ax.get_ylim()))/2 dx0 = np.abs(np.diff(ax.get_xlim()))/2 dxdy = [dx0, dy0] cx = np.mean(ax.get_xlim()) - 2.6 # zoom in between Sgr A and Sgr B2 cy = np.mean(ax.get_ylim()) rgbcmz = np.array(PIL.Image.open('gc_fullres_6.jpg'))[::-1,:,:] wwcmz = WCS(fits.Header.fromtextfile('gc_fullres_6.wcs')) def animate(n, nframes=nframes): dx0, dy0 = dxdy if n == 40: ax.imshow(atlasgal[0].data, norm=simple_norm(atlasgal[0].data, min_percent=5, max_percent=99.9, stretch='asinh'), transform=ax.get_transform(WCS(atlasgal[0].header)), cmap=orange_transparent) if n == 120: ax.imshow(rgbcmz, transform=ax.get_transform(wwcmz)) dy = dy0 * zoomfac[n] dx = dx0 * zoomfac[n] ax.axis((cx-dx, cx+dx, cy-dy, cy+dy)) print(f'{n}:{zoomfac[n]:0.1f}|', end='') return fig, anim = animation.FuncAnimation(fig, animate, frames=nframes, repeat_delay=5000, interval=50) anim.save('zoom_anim_cmz_linear_HI-CO-Dust_m0.35.gif')