######################################################## # Started Logging At: 2022-04-22 11:50:16 ######################################################## ######################################################## # # Started Logging At: 2022-04-22 11:50:17 ######################################################## ######################################################## # Started Logging At: 2022-04-22 13:44:40 ######################################################## ######################################################## # # Started Logging At: 2022-04-22 13:44:40 ######################################################## get_ipython().run_line_magic('matplotlib', 'inline') import pylab as pl pl.rcParams['figure.facecolor']='w' import os from astropy.utils.data import download_file from astropy.io import fits from spectral_cube import SpectralCube from astropy import units as u if os.path.exists("GC.hi.tb.allgal.fits.gz"): fh = fits.open("GC.hi.tb.allgal.fits.gz") hi_cube = SpectralCube.read(fh, use_dask=True) else: hi_cube = download_file("http://www.atnf.csiro.au/research/HI/sgps/SGPS_1/GC.hi.tb.allgal.fits.gz") hi_cube #[Out]# DaskSpectralCube with shape=(800, 1050, 1050) and unit=K and chunk size (200, 210, 210): #[Out]# n_x: 1050 type_x: GLON-CAR unit_x: deg range: 5.104167 deg: 354.905556 deg #[Out]# n_y: 1050 type_y: GLAT-CAR unit_y: deg range: -5.104167 deg: 5.094444 deg #[Out]# n_s: 800 type_s: VOPT unit_s: m / s range: -309144.300 m / s: 349534.122 m / s lv_hi_mx = hi_cube.max(axis=1) from astropy.visualization import simple_norm pl.imshow(lv_hi_mx.value, norm=simple_norm(lv_hi_mx.value, stretch='asinh', max_percent=99., min_percent=0.1)) #[Out]# if os.path.exists('GalacticCenterMosaic_downsampled4x_2kms.fits'): co_cube = SpectralCube.read('GalacticCenterMosaic_downsampled4x_2kms.fits', use_dask=True) else: #https://sedigism.mpifr-bonn.mpg.de/SEDIGISM_DATACUBES_DR1c/G000_13CO21_Tmb_DR1.fits raise NotImplementedError lv_co_mx = co_cube.max(axis=1) pl.imshow(lv_co_mx.value, norm=simple_norm(lv_co_mx.value, stretch='asinh', max_percent=99.5, min_percent=0.1)) #[Out]# if os.path.exists("CMZ_3mm_HNCO.fits"): hnco_cube = SpectralCube.read("CMZ_3mm_HNCO.fits", use_dask=True) else: raise NotImplementedError lv_hnco_mx = hnco_cube.max(axis=1) pl.imshow(lv_hnco_mx.value, norm=simple_norm(lv_hnco_mx.value, stretch='asinh', max_percent=99.5, min_percent=0.1)) #[Out]# if os.path.exists("DHT02_Center_interp.fits"): cube_12co = SpectralCube.read('DHT02_Center_interp.fits', use_dask=True) else: pass # https://dataverse.harvard.edu/api/access/datafile/:persistentId?persistentId=doi:10.7910/DVN/7RMW6Q/CA0UPU cube_12co #[Out]# DaskSpectralCube with shape=(492, 41, 201) and unit=K and chunk size (492, 41, 201): #[Out]# n_x: 201 type_x: GLON-CAR unit_x: deg range: 13.000000 deg: 348.000000 deg #[Out]# n_y: 41 type_y: GLAT-CAR unit_y: deg range: -2.500000 deg: 2.500000 deg #[Out]# n_s: 492 type_s: VOPT unit_s: m / s range: -319.239 m / s: 319.239 m / s lv_dameco_mx = cube_12co.max(axis=1) pl.imshow(lv_dameco_mx.value, norm=simple_norm(lv_dameco_mx.value, stretch='asinh', max_percent=99.5, min_percent=0.1)) #[Out]# import reproject import reproject.interpolation.core from astropy.wcs import WCS from scipy.ndimage import map_coordinates from reproject.interpolation.core import efficient_pixel_to_pixel_with_roundtrip # hack to get CO onto HI grid array = lv_co_mx.value wcs_in = lv_co_mx.wcs wcs_in.wcs.crval[0] = 0 # need to wrap at 0 wcs_out = hi_cube.with_spectral_unit(u.km/u.s, velocity_convention='radio')[:,0,:].wcs # need to match velo convention shape_out = lv_hi_mx.shape order = 1 array_out = None # Make sure image is floating point array = np.asarray(array, dtype=float) # shape_out must be exact a tuple type shape_out = tuple(shape_out) #_validate_array_out(array_out, array, shape_out) pixel_out = np.meshgrid(*[np.arange(size, dtype=float) for size in shape_out], indexing='ij', sparse=False, copy=False) pixel_out = [p.ravel() for p in pixel_out] # For each pixel in the ouput array, get the pixel value in the input WCS pixel_in = efficient_pixel_to_pixel_with_roundtrip(wcs_out, wcs_in, *pixel_out[::-1])[::-1] pixel_in = np.array(pixel_in) if array_out is not None: array_out.shape = (array_out.size,) else: array_out = np.empty(shape_out).ravel() # Interpolate array on to the pixels coordinates in pixel_in map_coordinates(array, pixel_in, order=order, cval=np.nan, mode='constant', output=array_out,).reshape(shape_out) array_out.shape = shape_out co_rep = array_out get_ipython().run_line_magic('matplotlib', 'inline') import pylab as pl pl.rcParams['figure.facecolor']='w' import numpy as np # hack to get CO onto HI grid array = lv_co_mx.value wcs_in = lv_co_mx.wcs wcs_in.wcs.crval[0] = 0 # need to wrap at 0 wcs_out = hi_cube.with_spectral_unit(u.km/u.s, velocity_convention='radio')[:,0,:].wcs # need to match velo convention shape_out = lv_hi_mx.shape order = 1 array_out = None # Make sure image is floating point array = np.asarray(array, dtype=float) # shape_out must be exact a tuple type shape_out = tuple(shape_out) #_validate_array_out(array_out, array, shape_out) pixel_out = np.meshgrid(*[np.arange(size, dtype=float) for size in shape_out], indexing='ij', sparse=False, copy=False) pixel_out = [p.ravel() for p in pixel_out] # For each pixel in the ouput array, get the pixel value in the input WCS pixel_in = efficient_pixel_to_pixel_with_roundtrip(wcs_out, wcs_in, *pixel_out[::-1])[::-1] pixel_in = np.array(pixel_in) if array_out is not None: array_out.shape = (array_out.size,) else: array_out = np.empty(shape_out).ravel() # Interpolate array on to the pixels coordinates in pixel_in map_coordinates(array, pixel_in, order=order, cval=np.nan, mode='constant', output=array_out,).reshape(shape_out) array_out.shape = shape_out co_rep = array_out pl.imshow(co_rep) #[Out]# # hack to get hnCO onto HI grid array = lv_hnco_mx.value wcs_in = lv_hnco_mx.wcs #wcs_in.wcs.crval[0] = 0 # need to wrap at 0 wcs_out = hi_cube.with_spectral_unit(u.km/u.s, velocity_convention='optical', rest_value=hi_cube.wcs.wcs.restfrq*u.Hz)[:,0,:].wcs shape_out = lv_hi_mx.shape order = 1 array_out = None # Make sure image is floating point array = np.asarray(array, dtype=float) # shape_out must be exact a tuple type shape_out = tuple(shape_out) #_validate_array_out(array_out, array, shape_out) pixel_out = np.meshgrid(*[np.arange(size, dtype=float) for size in shape_out], indexing='ij', sparse=False, copy=False) pixel_out = [p.ravel() for p in pixel_out] # For each pixel in the ouput array, get the pixel value in the input WCS #pixel_in = efficient_pixel_to_pixel_with_roundtrip(wcs_out, wcs_in, *pixel_out[::-1])[::-1] #pixel_in = np.array(pixel_in) # for reasons unclear, this was trickier xc,yc = wcs_out.wcs_pix2world(pixel_out[1], pixel_out[0], 0) xp,yp = wcs_in.wcs_world2pix(xc, yc, 0) pixel_in = [yp,xp] if array_out is not None: array_out.shape = (array_out.size,) else: array_out = np.empty(shape_out).ravel() # Interpolate array on to the pixels coordinates in pixel_in map_coordinates(array, pixel_in, order=order, cval=np.nan, mode='constant', output=array_out,).reshape(shape_out) array_out.shape = shape_out hnco_rep = array_out pl.imshow(hnco_rep) #[Out]# # hack to get CO onto HI grid array = lv_dameco_mx.value wcs_in = cube_12co[:,0,:].wcs wcs_in.wcs.cdelt[1] *= 1000 wcs_in.wcs.crval[1] *= 1000 #wcs_in.wcs.crval[0] = 0 # need to wrap at 0 wcs_out = lv_hi_mx.wcs #hi_cube.with_spectral_unit(u.km/u.s, velocity_convention='radio')[:,0,:].wcs # need to match velo convention shape_out = lv_hi_mx.shape order = 1 array_out = None # Make sure image is floating point array = np.asarray(array, dtype=float) # shape_out must be exact a tuple type shape_out = tuple(shape_out) #_validate_array_out(array_out, array, shape_out) pixel_out = np.meshgrid(*[np.arange(size, dtype=float) for size in shape_out], indexing='ij', sparse=False, copy=False) pixel_out = [p.ravel() for p in pixel_out] # For each pixel in the ouput array, get the pixel value in the input WCS #pixel_in = efficient_pixel_to_pixel_with_roundtrip(wcs_out, wcs_in, *pixel_out[::-1])[::-1] #pixel_in = np.array(pixel_in) xc,yc = wcs_out.wcs_pix2world(pixel_out[1], pixel_out[0], 0) xp,yp = wcs_in.wcs_world2pix(xc, yc, 0) pixel_in = [yp,xp] if array_out is not None: array_out.shape = (array_out.size,) else: array_out = np.empty(shape_out).ravel() # Interpolate array on to the pixels coordinates in pixel_in map_coordinates(array, pixel_in, order=order, cval=np.nan, mode='constant', output=array_out,).reshape(shape_out) array_out.shape = shape_out dameco_rep = array_out pl.imshow(dameco_rep) #[Out]# norm_hi = simple_norm(lv_hi_mx, stretch='log', max_percent=99.999, min_percent=15) norm_co = simple_norm(co_rep, stretch='asinh', max_percent=99.95, min_percent=0.1) norm_hnco = simple_norm(hnco_rep, stretch='asinh', max_percent=99.995, min_percent=0.1) norm_dameco = simple_norm(dameco_rep, stretch='asinh', max_percent=99.995, min_percent=0.1) rgb = np.array([norm_hi(lv_hi_mx.value)*0.7, norm_co(co_rep)*0.8, norm_dameco(dameco_rep), ]).swapaxes(0,2).swapaxes(0,1) pl.figure(figsize=(16,10)) slc = [slice(100,700), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(rgb[slc[0],slc[1],:]) #[Out]# from matplotlib.colors import hsv_to_rgb, rgb_to_hsv hsvim = rgb_to_hsv(rgb) hsvim[:,:,0] = (hsvim[:,:,0]-45/360.) hsvim[:,:,0][hsvim[:,:,0] > 1] -= 1 pl.figure(figsize=(16,10)) slc = [slice(100,700), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(hsv_to_rgb(hsvim)[slc[0],slc[1],:], origin='lower') #[Out]# red_hi = np.array([norm_hi(lv_hi_mx.value)*0.7, np.zeros_like(lv_hi_mx.value), np.zeros_like(lv_hi_mx).value]).swapaxes(0,2).swapaxes(0,1) hsv_hi = rgb_to_hsv(red_hi) hsv_hi[:,:,0] = (hsv_hi[:,:,0]+45/360.) orange_hi = hsv_to_rgb(hsv_hi) pl.imshow(orange_hi) #[Out]# def rotate_color(img, angle): img = np.nan_to_num(img) red_ = np.array([img, np.zeros_like(img), np.zeros_like(img)]).swapaxes(0,2).swapaxes(0,1) hsv_ = rgb_to_hsv(red_) hsv_[:,:,0] = (hsv_[:,:,0]+angle/360.) return hsv_to_rgb(hsv_) pl.imshow(rotate_color(norm_hi(lv_hi_mx.value)*0.7, 45)) #[Out]# pl.imshow(rotate_color(norm_co(co_rep), 270)) #[Out]# pl.imshow(rotate_color(norm_dameco(dameco_rep), 200)) #[Out]# new_rgb = rotate_color(norm_dameco(dameco_rep), 200)*0.7 + rotate_color(norm_co(co_rep), 270)*0.9 + rotate_color(norm_hi(lv_hi_mx.value)*0.6, 45) pl.figure(figsize=(20,10)) slc = [slice(100,700), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(new_rgb[slc[0], slc[1], :], origin='lower') ax.axes['glon'] = 'Galactic Longitude' ax.axes['velo'] = 'Velocity' ax.coords #[Out]# new_rgb = rotate_color(norm_dameco(dameco_rep), 200)*0.7 + rotate_color(norm_co(co_rep), 270)*0.9 + rotate_color(norm_hi(lv_hi_mx.value)*0.6, 45) pl.figure(figsize=(20,10)) slc = [slice(100,700), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(new_rgb[slc[0], slc[1], :], origin='lower') ax.coords['glon'] = 'Galactic Longitude' ax.coords['vopt'] = 'Velocity' new_rgb = rotate_color(norm_dameco(dameco_rep), 200)*0.7 + rotate_color(norm_co(co_rep), 270)*0.9 + rotate_color(norm_hi(lv_hi_mx.value)*0.6, 45) pl.figure(figsize=(20,10)) slc = [slice(100,700), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(new_rgb[slc[0], slc[1], :], origin='lower') ax.coords['glon'].set_label('Galactic Longitude') ax.coords['vopt'].set_label('Velocity') new_rgb = rotate_color(norm_dameco(dameco_rep), 200)*0.7 + rotate_color(norm_co(co_rep), 270)*0.9 + rotate_color(norm_hi(lv_hi_mx.value)*0.6, 45) pl.figure(figsize=(20,10)) slc = [slice(100,700), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(new_rgb[slc[0], slc[1], :], origin='lower') ax.coords['glon'].set_xlabel('Galactic Longitude') ax.coords['vopt'].set_ylabel('Velocity') vax = ax.coords['vopt'] new_rgb = rotate_color(norm_dameco(dameco_rep), 200)*0.7 + rotate_color(norm_co(co_rep), 270)*0.9 + rotate_color(norm_hi(lv_hi_mx.value)*0.6, 45) pl.figure(figsize=(20,10)) slc = [slice(100,700), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(new_rgb[slc[0], slc[1], :], origin='lower') ax.coords['glon'].set_axislabel('Galactic Longitude') ax.coords['vopt'].set_axislabel('Velocity') new_rgb = rotate_color(norm_dameco(dameco_rep), 200)*0.7 + rotate_color(norm_co(co_rep), 270)*0.9 + rotate_color(norm_hi(lv_hi_mx.value)*0.6, 45) pl.figure(figsize=(20,10)) slc = [slice(100,700), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(new_rgb[slc[0], slc[1], :], origin='lower') ax.coords['glon'].set_axislabel('Galactic Longitude') vax = ax.coords['vopt'] vax.set_axislabel('Velocity') vax.set_format_unit(u.km/u.s) new_rgb = rotate_color(norm_dameco(dameco_rep), 200)*0.7 + rotate_color(norm_co(co_rep), 270)*0.9 + rotate_color(norm_hi(lv_hi_mx.value)*0.6, 45) pl.figure(figsize=(20,10)) slc = [slice(100,700), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(new_rgb[slc[0], slc[1], :], origin='lower') ax.coords['glon'].set_axislabel('Galactic Longitude', fontsize=16) vax = ax.coords['vopt'] vax.set_axislabel('Velocity', fontsize=16) vax.set_format_unit(u.km/u.s) new_rgb = rotate_color(norm_dameco(dameco_rep), 200)*0.7 + rotate_color(norm_co(co_rep), 270)*0.9 + rotate_color(norm_hi(lv_hi_mx.value)*0.6, 45) pl.figure(figsize=(20,10)) slc = [slice(100,700), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(new_rgb[slc[0], slc[1], :], origin='lower') glax = ax.coords['glon'] glax.set_axislabel('Galactic Longitude', fontsize=16) glax.set_ticklabel(fontsize=16) vax = ax.coords['vopt'] vax.set_axislabel('Velocity', fontsize=16) vax.set_format_unit(u.km/u.s) vax.set_ticklabel(fontsize=16) new_rgb = rotate_color(norm_dameco(dameco_rep), 200)*0.7 + rotate_color(norm_co(co_rep), 270)*0.9 + rotate_color(norm_hi(lv_hi_mx.value)*0.6, 45) pl.figure(figsize=(20,10)) slc = [slice(100,700), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(new_rgb[slc[0], slc[1], :], origin='lower') glax = ax.coords['glon'] glax.set_axislabel('Galactic Longitude', fontsize=16) glax.set_ticklabel(fontsize=18) vax = ax.coords['vopt'] vax.set_axislabel('Velocity [km/s]', fontsize=18) vax.set_format_unit(u.km/u.s) vax.set_ticklabel(fontsize=16) new_rgb = rotate_color(norm_dameco(dameco_rep), 200)*0.7 + rotate_color(norm_co(co_rep), 270)*0.9 + rotate_color(norm_hi(lv_hi_mx.value)*0.6, 45) pl.figure(figsize=(20,10)) slc = [slice(0,800), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(new_rgb[slc[0], slc[1], :], origin='lower') glax = ax.coords['glon'] glax.set_axislabel('Galactic Longitude', fontsize=16) glax.set_ticklabel(fontsize=18) vax = ax.coords['vopt'] vax.set_axislabel('Velocity [km/s]', fontsize=18) vax.set_format_unit(u.km/u.s) vax.set_ticklabel(fontsize=16) new_rgb = rotate_color(norm_dameco(dameco_rep), 200)*0.7 + rotate_color(norm_co(co_rep), 270)*0.9 + rotate_color(norm_hi(lv_hi_mx.value)*0.6, 45) pl.figure(figsize=(20,10)) slc = [slice(50,750), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(new_rgb[slc[0], slc[1], :], origin='lower') glax = ax.coords['glon'] glax.set_axislabel('Galactic Longitude', fontsize=16) glax.set_ticklabel(fontsize=18) vax = ax.coords['vopt'] vax.set_axislabel('Velocity [km/s]', fontsize=18) vax.set_format_unit(u.km/u.s) vax.set_ticklabel(fontsize=16) new_rgb = rotate_color(norm_dameco(dameco_rep), 200)*0.7 + rotate_color(norm_co(co_rep), 270)*0.9 + rotate_color(norm_hi(lv_hi_mx.value)*0.6, 45) pl.figure(figsize=(26,10)) slc = [slice(50,750), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(new_rgb[slc[0], slc[1], :], origin='lower') glax = ax.coords['glon'] glax.set_axislabel('Galactic Longitude', fontsize=16) glax.set_ticklabel(fontsize=18) vax = ax.coords['vopt'] vax.set_axislabel('Velocity [km/s]', fontsize=18) vax.set_format_unit(u.km/u.s) vax.set_ticklabel(fontsize=16) new_rgb = rotate_color(norm_dameco(dameco_rep), 200)*0.7 + rotate_color(norm_co(co_rep), 270)*0.9 + rotate_color(norm_hi(lv_hi_mx.value)*0.6, 45) pl.figure(figsize=(20,10)) slc = [slice(50,750), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(new_rgb[slc[0], slc[1], :], origin='lower') glax = ax.coords['glon'] glax.set_axislabel('Galactic Longitude', fontsize=16) glax.set_ticklabel(fontsize=18) vax = ax.coords['vopt'] vax.set_axislabel('Velocity [km/s]', fontsize=18) vax.set_format_unit(u.km/u.s) vax.set_ticklabel(fontsize=16) new_rgb = rotate_color(norm_dameco(dameco_rep), 200)*0.7 + rotate_color(norm_co(co_rep), 270)*0.9 + rotate_color(norm_hi(lv_hi_mx.value)*0.6, 45) pl.figure(figsize=(20,10)) slc = [slice(50,750), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(new_rgb[slc[0], slc[1], :], origin='lower') glax = ax.coords['glon'] glax.set_axislabel('Galactic Longitude', fontsize=16) glax.set_ticklabel(fontsize=18) vax = ax.coords['vopt'] vax.set_axislabel('Velocity [km/s]', fontsize=18) vax.set_format_unit(u.km/u.s) vax.set_ticklabel(fontsize=16) pl.savefig("GC_PV_HI_12CO_CO32.png", bbox_inches='tight')