######################################################## # Started Logging At: 2022-02-03 07:22:08 ######################################################## ######################################################## # # Started Logging At: 2022-02-03 07:22:08 ######################################################## from astropy.utils.data import download_file hi_data = download_file("http://www.atnf.csiro.au/research/HI/sgps/SGPS_1/GC.hi.tb.allgal.fits.gz") get_ipython().run_line_magic('pinfo', 'download_file') if os.path.exists("GC.hi.tb.allgal.fits.gz"): hi_data = SpectralCube.read("GC.hi.tb.allgal.fits.gz") else: hi_data = download_file("http://www.atnf.csiro.au/research/HI/sgps/SGPS_1/GC.hi.tb.allgal.fits.gz") from astropy.utils.data import download_file from spectral_cube import SpectralCube if os.path.exists("GC.hi.tb.allgal.fits.gz"): hi_data = SpectralCube.read("GC.hi.tb.allgal.fits.gz") else: hi_data = download_file("http://www.atnf.csiro.au/research/HI/sgps/SGPS_1/GC.hi.tb.allgal.fits.gz") import os from astropy.utils.data import download_file from astropy.io import fits from spectral_cube import SpectralCube if os.path.exists("GC.hi.tb.allgal.fits.gz"): fh = fits.open("GC.hi.tb.allgal.fits.gz") hi_cube = SpectralCube.read(fh) else: hi_cube = download_file("http://www.atnf.csiro.au/research/HI/sgps/SGPS_1/GC.hi.tb.allgal.fits.gz") hi_cube #[Out]# SpectralCube with shape=(800, 1050, 1050) and unit=K: #[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 # CO hi_cube.beam #[Out]# Beam: BMAJ=145.00000476828 arcsec BMIN=145.00000476828 arcsec BPA=0.0 deg lv = cube.mean(axis=1) lv = hi_cube.mean(axis=1) 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_cube.mean(axis=1) pl.imshow(lv) get_ipython().run_line_magic('matplotlib', 'inline') import pylab as pl pl.rcParams['figure.facecolor']='w' pl.imshow(lv) #[Out]# pl.imshow(lv.value) #[Out]# from astropy.visualization import simple_norm pl.imshow(lv.value, norm=simple_norm(lv.value, stretch='asinh')) #[Out]# pl.imshow(lv.value, norm=simple_norm(lv.value, stretch='asinh', max_percent=99.95)) #[Out]# pl.imshow(lv.value, norm=simple_norm(lv.value, stretch='asinh', max_percent=99.9)) #[Out]# pl.imshow(lv.value, norm=simple_norm(lv.value, stretch='asinh', max_percent=99.9, min_percent=1)) #[Out]# pl.imshow(lv.value, norm=simple_norm(lv.value, stretch='asinh', max_percent=99.9, min_percent=0.1)) #[Out]# pl.imshow(lv.value, norm=simple_norm(lv.value, stretch='asinh', max_percent=99., min_percent=0.1)) #[Out]# lv = hi_cube.max(axis=1) pl.imshow(lv.value, norm=simple_norm(lv.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_hi_mx = hi_cube.max(axis=1) pl.imshow(lv_hi_mx.value, norm=simple_norm(lv_hi_mx.value, stretch='asinh', max_percent=99., min_percent=0.1)) #[Out]# 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]# co_rep = lv_co_mx.reproject(lv_hi_mx.header) co_rep import reproject import reproject.interpolation.core import reproject import reproject.interpolation.core from astropy.wcs import WCS reproject.interpolation.core._validate_wcs(lv_co_mx.wcs, lv_hi_mx.wcs) reproject.interpolation.core._validate_wcs(lv_co_mx.wcs, lv_hi_mx.wcs, [100,100]) hi_cube.with_spectral_unit(u.km/u.s, velocity_convention='radio') from astropy import units as u hi_cube.with_spectral_unit(u.km/u.s, velocity_convention='radio') #[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: VRAD-W2F unit_s: km / s range: -309.463 km / s: 349.127 km / 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]# reproject.interpolation.core._validate_wcs(lv_co_mx.wcs, lv_hi_mx.wcs, [100,100]) hi_cube = hi_cube.with_spectral_unit(u.km/u.s, velocity_convention='radio') 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]# reproject.interpolation.core._validate_wcs(lv_co_mx.wcs, lv_hi_mx.wcs, [100,100]) reproject.interpolation.core.efficient_pixel_to_pixel_with_roundtrip(lv_co_mx.wcs, lv_hi_mx.wcs, [5,5]) reproject.interpolation.core.efficient_pixel_to_pixel_with_roundtrip(lv_co_mx.wcs, lv_hi_mx.wcs, np.array([[5,5],[7,6]])) reproject.interpolation.core.efficient_pixel_to_pixel_with_roundtrip(lv_co_mx.wcs, lv_hi_mx.wcs, np.array([[5,5],[7,6], [8,4])) reproject.interpolation.core.efficient_pixel_to_pixel_with_roundtrip(lv_co_mx.wcs, lv_hi_mx.wcs, np.array([[5,5],[7,6], [8,4]])) reproject.interpolation.core.efficient_pixel_to_pixel_with_roundtrip(lv_co_mx.wcs, lv_hi_mx.wcs, *pixel_out[::-1]) shape_out = [20,21] pixel_out = np.meshgrid(*[np.arange(size, dtype=float) for size in shape_out], indexing='ij', sparse=False, copy=False) reproject.interpolation.core.efficient_pixel_to_pixel_with_roundtrip(lv_co_mx.wcs, lv_hi_mx.wcs, *pixel_out[::-1]) #[Out]# [array([[-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326], #[Out]# [-37119.71513923, -37118.62942493, -37117.54371063, #[Out]# -37116.45799633, -37115.37228204, -37114.28656774, #[Out]# -37113.20085344, -37112.11513914, -37111.02942484, #[Out]# -37109.94371054, -37108.85799625, -37107.77228195, #[Out]# -37106.68656765, -37105.60085335, -37104.51513905, #[Out]# -37103.42942475, -37102.34371046, -37101.25799616, #[Out]# -37100.17228186, -37099.08656756, -37098.00085326]]), #[Out]# array([[132.5576015 , 132.5576015 , 132.5576015 , 132.5576015 , #[Out]# 132.5576015 , 132.5576015 , 132.5576015 , 132.5576015 , #[Out]# 132.5576015 , 132.5576015 , 132.5576015 , 132.5576015 , #[Out]# 132.5576015 , 132.5576015 , 132.5576015 , 132.5576015 , #[Out]# 132.5576015 , 132.5576015 , 132.5576015 , 132.5576015 , #[Out]# 132.5576015 ], #[Out]# [133.76902369, 133.76902369, 133.76902369, 133.76902369, #[Out]# 133.76902369, 133.76902369, 133.76902369, 133.76902369, #[Out]# 133.76902369, 133.76902369, 133.76902369, 133.76902369, #[Out]# 133.76902369, 133.76902369, 133.76902369, 133.76902369, #[Out]# 133.76902369, 133.76902369, 133.76902369, 133.76902369, #[Out]# 133.76902369], #[Out]# [134.98045395, 134.98045395, 134.98045395, 134.98045395, #[Out]# 134.98045395, 134.98045395, 134.98045395, 134.98045395, #[Out]# 134.98045395, 134.98045395, 134.98045395, 134.98045395, #[Out]# 134.98045395, 134.98045395, 134.98045395, 134.98045395, #[Out]# 134.98045395, 134.98045395, 134.98045395, 134.98045395, #[Out]# 134.98045395], #[Out]# [136.19189229, 136.19189229, 136.19189229, 136.19189229, #[Out]# 136.19189229, 136.19189229, 136.19189229, 136.19189229, #[Out]# 136.19189229, 136.19189229, 136.19189229, 136.19189229, #[Out]# 136.19189229, 136.19189229, 136.19189229, 136.19189229, #[Out]# 136.19189229, 136.19189229, 136.19189229, 136.19189229, #[Out]# 136.19189229], #[Out]# [137.40333871, 137.40333871, 137.40333871, 137.40333871, #[Out]# 137.40333871, 137.40333871, 137.40333871, 137.40333871, #[Out]# 137.40333871, 137.40333871, 137.40333871, 137.40333871, #[Out]# 137.40333871, 137.40333871, 137.40333871, 137.40333871, #[Out]# 137.40333871, 137.40333871, 137.40333871, 137.40333871, #[Out]# 137.40333871], #[Out]# [138.6147932 , 138.6147932 , 138.6147932 , 138.6147932 , #[Out]# 138.6147932 , 138.6147932 , 138.6147932 , 138.6147932 , #[Out]# 138.6147932 , 138.6147932 , 138.6147932 , 138.6147932 , #[Out]# 138.6147932 , 138.6147932 , 138.6147932 , 138.6147932 , #[Out]# 138.6147932 , 138.6147932 , 138.6147932 , 138.6147932 , #[Out]# 138.6147932 ], #[Out]# [139.82625578, 139.82625578, 139.82625578, 139.82625578, #[Out]# 139.82625578, 139.82625578, 139.82625578, 139.82625578, #[Out]# 139.82625578, 139.82625578, 139.82625578, 139.82625578, #[Out]# 139.82625578, 139.82625578, 139.82625578, 139.82625578, #[Out]# 139.82625578, 139.82625578, 139.82625578, 139.82625578, #[Out]# 139.82625578], #[Out]# [141.03772642, 141.03772642, 141.03772642, 141.03772642, #[Out]# 141.03772642, 141.03772642, 141.03772642, 141.03772642, #[Out]# 141.03772642, 141.03772642, 141.03772642, 141.03772642, #[Out]# 141.03772642, 141.03772642, 141.03772642, 141.03772642, #[Out]# 141.03772642, 141.03772642, 141.03772642, 141.03772642, #[Out]# 141.03772642], #[Out]# [142.24920515, 142.24920515, 142.24920515, 142.24920515, #[Out]# 142.24920515, 142.24920515, 142.24920515, 142.24920515, #[Out]# 142.24920515, 142.24920515, 142.24920515, 142.24920515, #[Out]# 142.24920515, 142.24920515, 142.24920515, 142.24920515, #[Out]# 142.24920515, 142.24920515, 142.24920515, 142.24920515, #[Out]# 142.24920515], #[Out]# [143.46069195, 143.46069195, 143.46069195, 143.46069195, #[Out]# 143.46069195, 143.46069195, 143.46069195, 143.46069195, #[Out]# 143.46069195, 143.46069195, 143.46069195, 143.46069195, #[Out]# 143.46069195, 143.46069195, 143.46069195, 143.46069195, #[Out]# 143.46069195, 143.46069195, 143.46069195, 143.46069195, #[Out]# 143.46069195], #[Out]# [144.67218683, 144.67218683, 144.67218683, 144.67218683, #[Out]# 144.67218683, 144.67218683, 144.67218683, 144.67218683, #[Out]# 144.67218683, 144.67218683, 144.67218683, 144.67218683, #[Out]# 144.67218683, 144.67218683, 144.67218683, 144.67218683, #[Out]# 144.67218683, 144.67218683, 144.67218683, 144.67218683, #[Out]# 144.67218683], #[Out]# [145.88368979, 145.88368979, 145.88368979, 145.88368979, #[Out]# 145.88368979, 145.88368979, 145.88368979, 145.88368979, #[Out]# 145.88368979, 145.88368979, 145.88368979, 145.88368979, #[Out]# 145.88368979, 145.88368979, 145.88368979, 145.88368979, #[Out]# 145.88368979, 145.88368979, 145.88368979, 145.88368979, #[Out]# 145.88368979], #[Out]# [147.09520082, 147.09520082, 147.09520082, 147.09520082, #[Out]# 147.09520082, 147.09520082, 147.09520082, 147.09520082, #[Out]# 147.09520082, 147.09520082, 147.09520082, 147.09520082, #[Out]# 147.09520082, 147.09520082, 147.09520082, 147.09520082, #[Out]# 147.09520082, 147.09520082, 147.09520082, 147.09520082, #[Out]# 147.09520082], #[Out]# [148.30671993, 148.30671993, 148.30671993, 148.30671993, #[Out]# 148.30671993, 148.30671993, 148.30671993, 148.30671993, #[Out]# 148.30671993, 148.30671993, 148.30671993, 148.30671993, #[Out]# 148.30671993, 148.30671993, 148.30671993, 148.30671993, #[Out]# 148.30671993, 148.30671993, 148.30671993, 148.30671993, #[Out]# 148.30671993], #[Out]# [149.51824712, 149.51824712, 149.51824712, 149.51824712, #[Out]# 149.51824712, 149.51824712, 149.51824712, 149.51824712, #[Out]# 149.51824712, 149.51824712, 149.51824712, 149.51824712, #[Out]# 149.51824712, 149.51824712, 149.51824712, 149.51824712, #[Out]# 149.51824712, 149.51824712, 149.51824712, 149.51824712, #[Out]# 149.51824712], #[Out]# [150.72978238, 150.72978238, 150.72978238, 150.72978238, #[Out]# 150.72978238, 150.72978238, 150.72978238, 150.72978238, #[Out]# 150.72978238, 150.72978238, 150.72978238, 150.72978238, #[Out]# 150.72978238, 150.72978238, 150.72978238, 150.72978238, #[Out]# 150.72978238, 150.72978238, 150.72978238, 150.72978238, #[Out]# 150.72978238], #[Out]# [151.94132573, 151.94132573, 151.94132573, 151.94132573, #[Out]# 151.94132573, 151.94132573, 151.94132573, 151.94132573, #[Out]# 151.94132573, 151.94132573, 151.94132573, 151.94132573, #[Out]# 151.94132573, 151.94132573, 151.94132573, 151.94132573, #[Out]# 151.94132573, 151.94132573, 151.94132573, 151.94132573, #[Out]# 151.94132573], #[Out]# [153.15287715, 153.15287715, 153.15287715, 153.15287715, #[Out]# 153.15287715, 153.15287715, 153.15287715, 153.15287715, #[Out]# 153.15287715, 153.15287715, 153.15287715, 153.15287715, #[Out]# 153.15287715, 153.15287715, 153.15287715, 153.15287715, #[Out]# 153.15287715, 153.15287715, 153.15287715, 153.15287715, #[Out]# 153.15287715], #[Out]# [154.36443664, 154.36443664, 154.36443664, 154.36443664, #[Out]# 154.36443664, 154.36443664, 154.36443664, 154.36443664, #[Out]# 154.36443664, 154.36443664, 154.36443664, 154.36443664, #[Out]# 154.36443664, 154.36443664, 154.36443664, 154.36443664, #[Out]# 154.36443664, 154.36443664, 154.36443664, 154.36443664, #[Out]# 154.36443664], #[Out]# [155.57600422, 155.57600422, 155.57600422, 155.57600422, #[Out]# 155.57600422, 155.57600422, 155.57600422, 155.57600422, #[Out]# 155.57600422, 155.57600422, 155.57600422, 155.57600422, #[Out]# 155.57600422, 155.57600422, 155.57600422, 155.57600422, #[Out]# 155.57600422, 155.57600422, 155.57600422, 155.57600422, #[Out]# 155.57600422]])] cube = SpectralCube.read('/orange/adamginsburg/sgrb2/2013.1.00269.S/HC3N_TP_7m_12m_feather.fits', use_dask=True) cube = SpectralCube.read('/orange/adamginsburg/sgrb2/2013.1.00269.S/HC3N_TP_7m_12m_feather_ds.fits', use_dask=True) m0 = cube.moment0() m0.quicklook() pl.figure() ax = pl.subplot(projection=m0.wcs) pl.figure() ax = pl.subplot(projection=m0.wcs) ax.imshow(m0.value, norm=simple_norm(m0.value, stretch='asinh', max_percent=99.995), cmap='gray') #[Out]# pl.figure() ax = pl.subplot(projection=m0.wcs) ax.imshow(m0.value, norm=simple_norm(m0.value, stretch='asinh', max_percent=99.995), cmap='gray_r') #[Out]# pl.figure(figsize=(20,20)) ax = pl.subplot(projection=m0.wcs) ax.imshow(m0.value, norm=simple_norm(m0.value, stretch='asinh', max_percent=99.995), cmap='gray_r') #[Out]# pl.figure(figsize=(20,20)) ax = pl.subplot(projection=m0.wcs) ax.imshow(m0.value, norm=simple_norm(m0.value, stretch='asinh', max_percent=99.995, min_percent=0.01), cmap='gray_r') #[Out]# cube = SpectralCube.read('/orange/adamginsburg/sgrb2/2013.1.00269.S/HC3N_TP_7m_12m_feather_ds.fits', use_dask=True) print(cube) cube = SpectralCube.read("/orange/adamginsburg/sgrb2/2013.1.00269.S/merge/SgrB2_b3_7M_12M.HC3N.r2.fits", use_dask=True) print(cube) cube = SpectralCube.read('/orange/adamginsburg/sgrb2/2013.1.00269.S/HC3N_TP_7m_12m_feather_ds.fits', use_dask=True) print(cube) cube = SpectralCube.read("/orange/adamginsburg/sgrb2/2013.1.00269.S/tp_merge/HC3N_TP_7m_12m_feather.fits", use_dask=True) print(cube) cube = SpectralCube.read('/orange/adamginsburg/sgrb2/2013.1.00269.S/HC3N_TP_7m_12m_feather_ds.fits', use_dask=True) print(cube) cube = SpectralCube.read("/orange/adamginsburg/sgrb2/2013.1.00269.S/tp_merge/HC3N_TP_7m_12m_feather_r05.fits", use_dask=True) print(cube) cube = SpectralCube.read("/orange/adamginsburg/sgrb2/2013.1.00269.S/tp_merge/HC3N_TP_7m_12m_feather.fits", use_dask=True) print(cube) m0 = cube.moment0() pl.figure(figsize=(20,20)) ax = pl.subplot(projection=m0.wcs) ax.imshow(m0.value, norm=simple_norm(m0.value, stretch='asinh', max_percent=99.995, min_percent=0.01), cmap='gray_r') #[Out]# pl.figure(figsize=(20,20)) ax = pl.subplot(projection=m0.wcs) ax.imshow(m0.value, norm=simple_norm(m0.value, stretch='asinh', max_percent=99.5, min_percent=0.01), cmap='gray_r') #[Out]# pl.figure(figsize=(20,20)) ax = pl.subplot(projection=m0.wcs) ax.imshow(m0.value, norm=simple_norm(m0.value, stretch='asinh', max_percent=99.5, min_percent=0.0001), cmap='gray_r') #[Out]# pl.figure(figsize=(20,20)) ax = pl.subplot(projection=m0.wcs) ax.imshow(m0.value, norm=simple_norm(m0.value, stretch='log', max_percent=99.5, min_percent=0.0001), cmap='gray_r') #[Out]# pl.figure(figsize=(20,20)) ax = pl.subplot(projection=m0.wcs) ax.imshow(m0.value, norm=simple_norm(m0.value, stretch='log', max_percent=99., min_percent=0.0001), cmap='gray_r') #[Out]# pl.figure(figsize=(20,20)) ax = pl.subplot(projection=m0.wcs) ax.imshow(m0.value, norm=simple_norm(m0.value, stretch='log', max_percent=90., min_percent=0.0001), cmap='gray_r') #[Out]# pl.figure(figsize=(20,20)) ax = pl.subplot(projection=m0.wcs) ax.imshow(m0.value, norm=simple_norm(m0.value, stretch='log', max_percent=90., min_percent=5), cmap='gray_r') #[Out]# pl.figure(figsize=(20,20)) ax = pl.subplot(projection=m0.wcs) ax.imshow(m0.value, norm=simple_norm(m0.value, stretch='log', max_percent=99.9, min_percent=5), cmap='gray_r') #[Out]# get_ipython().run_line_magic('ls', '-lh /orange/adamginsburg/sgrb2/2013.1.00269.S/tp_merge/H') get_ipython().run_line_magic('ls', '-lh /orange/adamginsburg/sgrb2/2013.1.00269.S/tp_merge/') get_ipython().run_line_magic('ls', '-lh /orange/adamginsburg/sgrb2/2013.1.00269.S/tp_merge/HC3N_TP_7m_12m_feather') m0 = cube.spectral_slab(-8*u.km/u.s, 138*u.km/u.s).moment0() pl.figure(figsize=(20,20)) ax = pl.subplot(projection=m0.wcs) ax.imshow(m0.value, norm=simple_norm(m0.value, stretch='log', max_percent=99.95, min_percent=1), cmap='gray_r') #[Out]# pl.figure(figsize=(20,20)) ax = pl.subplot(projection=m0.wcs) ax.imshow(m0.value, norm=simple_norm(m0.value, stretch='log', max_percent=99.95, min_percent=2), cmap='gray_r') #[Out]# pl.figure(figsize=(20,20)) ax = pl.subplot(projection=m0.wcs) ax.imshow(m0.value, norm=simple_norm(m0.value, stretch='log', max_percent=99.95, min_percent=5), cmap='gray_r') #[Out]# pl.figure(figsize=(20,20)) ax = pl.subplot(projection=m0.wcs) ax.imshow(m0.value, norm=simple_norm(m0.value, stretch='log', max_percent=99.995, min_percent=5), cmap='gray_r') #[Out]# pl.figure(figsize=(20,20)) ax = pl.subplot(projection=m0.wcs) ax.imshow(m0.value, norm=simple_norm(m0.value, stretch='asinh', max_percent=99.995, min_percent=5), cmap='gray_r') #[Out]# import reproject import reproject.interpolation.core from astropy.wcs import WCS from scipy.ndimage import map_coordinates array = lv_co_mx.value wcs_in = lv_co_mx.wcs wcs_out = lv_hi_mx.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) 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 array = lv_co_mx.value wcs_in = lv_co_mx.wcs wcs_out = lv_hi_mx.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) 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 import reproject import reproject.interpolation.core from astropy.wcs import WCS from scipy.ndimage import map_coordinates from import reproject.interpolation.core import efficient_pixel_to_pixel_with_roundtrip 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 array = lv_co_mx.value wcs_in = lv_co_mx.wcs wcs_out = lv_hi_mx.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) 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 ax.imshow(co_rep) # hack to get CO onto HI grid array = lv_co_mx.value wcs_in = lv_co_mx.wcs wcs_out = lv_hi_mx.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) 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 # hack to get CO onto HI grid array = lv_co_mx.value wcs_in = lv_co_mx.wcs wcs_out = lv_hi_mx.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) 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 ax.imshow(co_rep) #[Out]# pl.imshow(co_rep) #[Out]# co_rep #[Out]# array([[nan, nan, nan, ..., nan, nan, nan], #[Out]# [nan, nan, nan, ..., nan, nan, nan], #[Out]# [nan, nan, nan, ..., nan, nan, nan], #[Out]# ..., #[Out]# [nan, nan, nan, ..., nan, nan, nan], #[Out]# [nan, nan, nan, ..., nan, nan, nan], #[Out]# [nan, nan, nan, ..., nan, nan, nan]]) pixel_out #[Out]# [array([ 0., 0., 0., ..., 799., 799., 799.]), #[Out]# array([0.000e+00, 1.000e+00, 2.000e+00, ..., 1.047e+03, 1.048e+03, #[Out]# 1.049e+03])] pixel_in #[Out]# array([[ -109.46341694, -109.46341694, -109.46341694, ..., #[Out]# 549.12706715, 549.12706715, 549.12706715], #[Out]# [34189.21091585, 34190.13196847, 34191.05302109, ..., #[Out]# 35153.55300992, 35154.47406254, 35155.39511517]]) pixel_in.shape #[Out]# (2, 840000) pl.plot(pixel_in[0,:], pixel_in[1,:], ',') #[Out]# [] wcs_out.wcs_pix2world(0,0,0) #[Out]# [array(5.10416655), array(-309463.41693863)] wcs_in.wcs_pix2world(0,0,0) #[Out]# [array(365.99027771), array(-200000.)] wcs_out.wcs_pix2world(500,0,0) #[Out]# [array(0.24305555), array(-309463.41693863)] wcs_out.wcs_pix2world(500,200,0) #[Out]# [array(0.24305555), array(-144338.05955191)] wcs_out.wcs_pix2world(500,500,0) #[Out]# [array(0.24305555), array(103009.54346129)] wcs_out.wcs_pix2world(500,400,0) #[Out]# [array(0.24305555), array(20605.68360877)] wcs_in.wcs_world2pix(-0.24, 20605, 0) #[Out]# [array(34695.50038418), array(220.605)] wcs_in.wcs_world2pix(360-0.24, 20605, 0) #[Out]# [array(590.23684236), array(220.605)] wcs_out.wcs_pix2world(600,400,0) #[Out]# [array(-0.72916665), array(20605.68360877)] wcs_in #[Out]# WCS Keywords #[Out]# #[Out]# Number of WCS axes: 2 #[Out]# CTYPE : 'GLON' 'VRAD' #[Out]# CRVAL : 360.0 -200000.0 #[Out]# CRPIX : 568.5 1.0 #[Out]# PC1_1 PC1_2 : 1.0 0.0 #[Out]# PC2_1 PC2_2 : 0.0 1.0 #[Out]# CDELT : -0.010555555436732 1000.0 #[Out]# NAXIS : 0 0 wcs_in.wcs.crval[0] #[Out]# 360.0 wcs_in.wcs.crval[0] = 0 wcs_in.wcs_pix2world(0,0,0) #[Out]# [array(5.99027771), array(-200000.)] # 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 = lv_hi_mx.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) 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.plot(pixel_in[0,:], pixel_in[1,:], ',') #[Out]# [] pl.imshow(co_rep) #[Out]# # hack to get CO 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 = lv_hi_mx.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) 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_hnco_mx.value wcs_in = lv_hnco_mx.wcs #wcs_in.wcs.crval[0] = 0 # need to wrap at 0 wcs_out = lv_hi_mx.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) 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]# pl.plot(pixel_in[0,:], pixel_in[1,:], ',') #[Out]# [] pixel_in #[Out]# array([[ nan, nan, nan, ..., #[Out]# nan, nan, nan], #[Out]# [-1149.2498712 , -1146.33320478, -1143.41653835, ..., #[Out]# 1904.49987193, 1907.41653835, 1910.33320478]]) wcs_in #[Out]# WCS Keywords #[Out]# #[Out]# Number of WCS axes: 2 #[Out]# CTYPE : 'GLON' 'VOPT' #[Out]# CRVAL : 0.0 46.6396151787 #[Out]# CRPIX : 383.0 164.5 #[Out]# PC1_1 PC1_2 : 1.0 0.0 #[Out]# PC2_1 PC2_2 : 0.0 1.0 #[Out]# CDELT : -0.00333333353753 1837.8933618 #[Out]# NAXIS : 0 0 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: VRAD-W2F unit_s: km / s range: -309.463 km / s: 349.127 km / 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]# 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 = lv_hi_mx.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) 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]# hnco_rep = lv_hnco_mx.reproject(lv_hi_mx.header) hnco_rep # 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 = lv_hi_mx.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) 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.plot(pixel_in[0,:], pixel_in[1,:], ',') #[Out]# [] wcs_in.wcs_pix2world(50,50,0) #[Out]# [array(1.10666673), array(-208554.25694912)] wcs_in.wcs_world2pix(0.24, 50000) wcs_in.wcs_world2pix(0.24, 50000, 0) #[Out]# [array(310.00000441), array(190.67968377)] shape_out #[Out]# (800, 1050) pixel_out #[Out]# [array([ 0., 0., 0., ..., 799., 799., 799.]), #[Out]# array([0.000e+00, 1.000e+00, 2.000e+00, ..., 1.047e+03, 1.048e+03, #[Out]# 1.049e+03])] efficient_pixel_to_pixel_with_roundtrip(wcs_out, wcs_in, *pixel_out[::-1])[::-1] #[Out]# [array([nan, nan, nan, ..., nan, nan, nan]), #[Out]# array([-1149.2498712 , -1146.33320478, -1143.41653835, ..., #[Out]# 1904.49987193, 1907.41653835, 1910.33320478])] wcs_out.pixel_to_world(pixel_out[0,:], pixel_out[1,:]) wcs_out.pixel_to_world(pixel_out[0], pixel_out[1]) #[Out]# [, #[Out]# ] pinx, piny = wcs_out.pixel_to_world(pixel_out[0], pixel_out[1]) pinx, piny #[Out]# (, #[Out]# ) wcs_in.world_to_pixel(pinx, piny) #[Out]# (array([-1149.2498712 , -1149.2498712 , -1149.2498712 , ..., #[Out]# 1181.16659944, 1181.16659944, 1181.16659944]), #[Out]# array([nan, nan, nan, ..., nan, nan, nan])) 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) # 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 = lv_hi_mx.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) 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 = lv_hi_mx.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) 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 # 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')[:,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) 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 wcs_in.wcs_pix2world(50,50,0) #[Out]# [array(1.10666673), array(-208554.25694912)] pinx, piny = wcs_out.pixel_to_world(pixel_out[0], pixel_out[1]) pinx, piny #[Out]# (, #[Out]# ) wcs_in.world_to_pixel(pinx, piny) #[Out]# (array([-1149.2498712 , -1149.2498712 , -1149.2498712 , ..., #[Out]# 1181.16659944, 1181.16659944, 1181.16659944]), #[Out]# array([nan, nan, nan, ..., nan, nan, nan])) wcs_in.wcs_world2pix(0.24, 50000, 0) #[Out]# [array(310.00000441), array(190.67968377)] efficient_pixel_to_pixel_with_roundtrip(wcs_out, wcs_in, *pixel_out[::-1])[::-1] #[Out]# [array([nan, nan, nan, ..., nan, nan, nan]), #[Out]# array([-1149.2498712 , -1146.33320478, -1143.41653835, ..., #[Out]# 1904.49987193, 1907.41653835, 1910.33320478])] wcs_in #[Out]# WCS Keywords #[Out]# #[Out]# Number of WCS axes: 2 #[Out]# CTYPE : 'GLON' 'VOPT' #[Out]# CRVAL : 0.0 46.6396151787 #[Out]# CRPIX : 383.0 164.5 #[Out]# PC1_1 PC1_2 : 1.0 0.0 #[Out]# PC2_1 PC2_2 : 0.0 1.0 #[Out]# CDELT : -0.00333333353753 1837.8933618 #[Out]# NAXIS : 0 0 pl.plot(pixel_in[0,:], pixel_in[1,:], ',') #[Out]# [] hi_cube.wcs.wcs.restfrq #[Out]# 1420405758.58 # 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.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) 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 # 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) 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 wcs_in.world_to_pixel(pinx, piny) #[Out]# (array([-1149.2498712 , -1149.2498712 , -1149.2498712 , ..., #[Out]# 1181.16659944, 1181.16659944, 1181.16659944]), #[Out]# array([nan, nan, nan, ..., nan, nan, nan])) wcs_in.wcs_pix2world(50,50,0) #[Out]# [array(1.10666673), array(-208554.25694912)] pinx, piny = wcs_out.pixel_to_world(pixel_out[0], pixel_out[1]) pinx, piny #[Out]# (, #[Out]# ) wcs_in.world_to_pixel(pinx, piny) #[Out]# (array([-1149.2498712 , -1149.2498712 , -1149.2498712 , ..., #[Out]# 1181.16659944, 1181.16659944, 1181.16659944]), #[Out]# array([nan, nan, nan, ..., nan, nan, nan])) pinx, piny = wcs_out.pixel_to_world(pixel_out[1], pixel_out[0]) pinx, piny #[Out]# (, #[Out]# ) wcs_in.world_to_pixel(pinx, piny) #[Out]# (array([-1149.2498712 , -1146.33320478, -1143.41653835, ..., #[Out]# 1904.49987193, 1907.41653835, 1910.33320478]), #[Out]# array([nan, nan, nan, ..., nan, nan, nan])) # 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) xc,yc = wcs_out.wcs_pix2world(pixel_out[1], pixel_out[0], 0) pixel_in = wcs_in.wcs_world2pix(xc, yc, 0) 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.plot(pixel_in[0,:], pixel_in[1,:], ',') pixel_in #[Out]# [array([-1149.2498712 , -1146.33320478, -1143.41653835, ..., #[Out]# 1904.49987193, 1907.41653835, 1910.33320478]), #[Out]# array([ -4.73116403, -4.73116403, -4.73116403, ..., 354.02612841, #[Out]# 354.02612841, 354.02612841])] pl.plot(pixel_in[0], pixel_in[1], ',') #[Out]# [] pl.imshow(hnco_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) xc,yc = wcs_out.wcs_pix2world(pixel_out[0], pixel_out[1], 0) pixel_in = wcs_in.wcs_world2pix(xc, yc, 0) 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 pixel_in #[Out]# [array([-1149.2498712 , -1149.2498712 , -1149.2498712 , ..., #[Out]# 1181.16659944, 1181.16659944, 1181.16659944]), #[Out]# array([ -4.73116403, -4.28215615, -3.83314828, ..., 465.38008151, #[Out]# 465.82908938, 466.27809726])] pl.plot(pixel_in[0], pixel_in[1], ',') #[Out]# [] pl.imshow(hnco_rep) #[Out]# hnco_cube #[Out]# DaskSpectralCube with shape=(327, 165, 765) and unit=K and chunk size (327, 165, 621): #[Out]# n_x: 765 type_x: GLON-SIN unit_x: deg range: 1.815115 deg: 359.268218 deg #[Out]# n_y: 165 type_y: GLAT-SIN unit_y: deg range: -0.314991 deg: 0.231678 deg #[Out]# n_s: 327 type_s: VOPT unit_s: m / s range: -300448.925 m / s: 298704.311 m / s pixel_out #[Out]# [array([ 0., 0., 0., ..., 799., 799., 799.]), #[Out]# array([0.000e+00, 1.000e+00, 2.000e+00, ..., 1.047e+03, 1.048e+03, #[Out]# 1.049e+03])] shape_out #[Out]# (800, 1050) lv_hi_mx.shape #[Out]# (800, 1050) # 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) 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]# 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 hnco_cube #[Out]# DaskSpectralCube with shape=(327, 165, 765) and unit=K and chunk size (327, 165, 621): #[Out]# n_x: 765 type_x: GLON-SIN unit_x: deg range: 1.815115 deg: 359.268218 deg #[Out]# n_y: 165 type_y: GLAT-SIN unit_y: deg range: -0.314991 deg: 0.231678 deg #[Out]# n_s: 327 type_s: VOPT unit_s: m / s range: -300448.925 m / s: 298704.311 m / s rgb = [lv_hi_mx.value, co_rep, hnco_rep] rgb = np.array([lv_hi_mx.value, co_rep, hnco_rep]).swapaxes(0,2).swapaxes(0,1) pl.imshow(rgb) #[Out]# norm_hi = simple_norm(lv_hi_mx) norm_co = simple_norm(co_rep) norm_hnco = simple_norm(hnco_rep) rgb = np.array([norm_hi(lv_hi_mx.value), norm_co(co_rep), norm_nhco(hnco_rep)]).swapaxes(0,2).swapaxes(0,1) norm_hi = simple_norm(lv_hi_mx) norm_co = simple_norm(co_rep) norm_hnco = simple_norm(hnco_rep) rgb = np.array([norm_hi(lv_hi_mx.value), norm_co(co_rep), norm_hnco(hnco_rep)]).swapaxes(0,2).swapaxes(0,1) pl.imshow(rgb) #[Out]# pl.figure(figsize=(16,10)) ax = pl.subplot(projection=lv_hi_mx.wcs) ax.imshow(rgb) #[Out]# norm_hi = simple_norm(lv_hi_mx, stretch='asinh', max_percent=99., min_percent=0.1) norm_co = simple_norm(co_rep) norm_hnco = simple_norm(hnco_rep) rgb = np.array([norm_hi(lv_hi_mx.value), norm_co(co_rep), norm_hnco(hnco_rep)]).swapaxes(0,2).swapaxes(0,1) pl.figure(figsize=(16,10)) ax = pl.subplot(projection=lv_hi_mx.wcs) ax.imshow(rgb) #[Out]# pl.figure(figsize=(16,10)) slc = [slice(100,700), slice(300,700)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(rgb[slc[0],slc[1],:]) #[Out]# norm_hi = simple_norm(lv_hi_mx, stretch='asinh', max_percent=99., min_percent=0.1) 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) rgb = np.array([norm_hi(lv_hi_mx.value), norm_co(co_rep), norm_hnco(hnco_rep)]).swapaxes(0,2).swapaxes(0,1) pl.figure(figsize=(16,10)) slc = [slice(100,700), slice(300,700)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(rgb[slc[0],slc[1],:]) #[Out]# norm_hi = simple_norm(lv_hi_mx, stretch='asinh', max_percent=99., min_percent=0.1) 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) rgb = np.array([norm_hi(lv_hi_mx.value), norm_hnco(hnco_rep), norm_co(co_rep),]).swapaxes(0,2).swapaxes(0,1) pl.figure(figsize=(16,10)) slc = [slice(100,700), slice(300,700)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(rgb[slc[0],slc[1],:]) #[Out]# norm_hi = simple_norm(lv_hi_mx, stretch='asinh', max_percent=99., min_percent=0.1) norm_co = simple_norm(co_rep, stretch='asinh', max_percent=99.5, min_percent=0.1) norm_hnco = simple_norm(hnco_rep, stretch='asinh', max_percent=99.995, min_percent=0.1) rgb = np.array([norm_hi(lv_hi_mx.value), norm_hnco(hnco_rep), norm_co(co_rep),]).swapaxes(0,2).swapaxes(0,1) norm_hi = simple_norm(lv_hi_mx, stretch='asinh', max_percent=99., min_percent=0.1) norm_co = simple_norm(co_rep, stretch='asinh', max_percent=99.5, min_percent=0.1) norm_hnco = simple_norm(hnco_rep, stretch='asinh', max_percent=99.995, min_percent=0.1) rgb = np.array([norm_hi(lv_hi_mx.value), norm_hnco(hnco_rep), norm_co(co_rep),]).swapaxes(0,2).swapaxes(0,1) pl.figure(figsize=(16,10)) slc = [slice(100,700), slice(300,700)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(rgb[slc[0],slc[1],:]) #[Out]# norm_hi = simple_norm(lv_hi_mx, stretch='asinh', max_percent=99., min_percent=1) norm_co = simple_norm(co_rep, stretch='asinh', max_percent=99.5, min_percent=0.1) norm_hnco = simple_norm(hnco_rep, stretch='asinh', max_percent=99.995, min_percent=0.1) rgb = np.array([norm_hi(lv_hi_mx.value), norm_hnco(hnco_rep), norm_co(co_rep),]).swapaxes(0,2).swapaxes(0,1) norm_hi = simple_norm(lv_hi_mx, stretch='asinh', max_percent=99., min_percent=1) norm_co = simple_norm(co_rep, stretch='asinh', max_percent=99.5, min_percent=0.1) norm_hnco = simple_norm(hnco_rep, stretch='asinh', max_percent=99.995, min_percent=0.1) rgb = np.array([norm_hi(lv_hi_mx.value), norm_hnco(hnco_rep), norm_co(co_rep),]).swapaxes(0,2).swapaxes(0,1) pl.figure(figsize=(16,10)) slc = [slice(100,700), slice(300,700)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(rgb[slc[0],slc[1],:]) #[Out]# 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]# def make_hsv_fig(img): hsvim = rgb_to_hsv(img) fig = pl.figure(figsize=(16,16), facecolor='k') for ii,angle in enumerate(np.linspace(0,120,9)): pl.subplot(3,3,ii+1) hsvim[:,:,0] = (hsvim[:,:,0]+45/360.) hsvim[:,:,0][hsvim[:,:,0] > 1] -= 1 pl.imshow(hsv_to_rgb(hsvim), origin='lower') pl.xticks([]) pl.yticks([]) pl.tight_layout() return fig make_hsv_fig(rgb) from matplotlib.colors import hsv_to_rgb, rgb_to_hsv def make_hsv_fig(img): hsvim = rgb_to_hsv(img) fig = pl.figure(figsize=(16,16), facecolor='k') for ii,angle in enumerate(np.linspace(0,120,9)): pl.subplot(3,3,ii+1) hsvim[:,:,0] = (hsvim[:,:,0]+45/360.) hsvim[:,:,0][hsvim[:,:,0] > 1] -= 1 pl.imshow(hsv_to_rgb(hsvim), origin='lower') pl.xticks([]) pl.yticks([]) pl.tight_layout() return fig make_hsv_fig(rgb) #[Out]#
norm_hi = simple_norm(lv_hi_mx, stretch='asinh', max_percent=99., min_percent=1) norm_co = simple_norm(co_rep, stretch='asinh', max_percent=99.5, min_percent=0.1) norm_hnco = simple_norm(hnco_rep, stretch='asinh', max_percent=99.995, min_percent=0.1) rgb = np.array([norm_hi(lv_hi_mx.value), norm_co(co_rep), norm_hnco(hnco_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 def make_hsv_fig(img): hsvim = rgb_to_hsv(img) fig = pl.figure(figsize=(16,16), facecolor='k') for ii,angle in enumerate(np.linspace(0,120,9)): pl.subplot(3,3,ii+1) hsvim[:,:,0] = (hsvim[:,:,0]+45/360.) hsvim[:,:,0][hsvim[:,:,0] > 1] -= 1 pl.imshow(hsv_to_rgb(hsvim), origin='lower') pl.xticks([]) pl.yticks([]) pl.tight_layout() return fig make_hsv_fig(rgb) #[Out]#
norm_hi = simple_norm(lv_hi_mx, stretch='log', max_percent=99., min_percent=1) norm_co = simple_norm(co_rep, stretch='asinh', max_percent=99.5, min_percent=0.1) norm_hnco = simple_norm(hnco_rep, stretch='asinh', max_percent=99.995, min_percent=0.1) rgb = np.array([norm_hi(lv_hi_mx.value), norm_co(co_rep), norm_hnco(hnco_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]# norm_hi = simple_norm(lv_hi_mx, stretch='log', max_percent=99., min_percent=5) norm_co = simple_norm(co_rep, stretch='asinh', max_percent=99.5, min_percent=0.1) norm_hnco = simple_norm(hnco_rep, stretch='asinh', max_percent=99.995, min_percent=0.1) rgb = np.array([norm_hi(lv_hi_mx.value), norm_co(co_rep), norm_hnco(hnco_rep), ]).swapaxes(0,2).swapaxes(0,1) norm_hi = simple_norm(lv_hi_mx, stretch='log', max_percent=99., min_percent=5) norm_co = simple_norm(co_rep, stretch='asinh', max_percent=99.5, min_percent=0.1) norm_hnco = simple_norm(hnco_rep, stretch='asinh', max_percent=99.995, min_percent=0.1) rgb = np.array([norm_hi(lv_hi_mx.value), norm_co(co_rep), norm_hnco(hnco_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]# norm_hi = simple_norm(lv_hi_mx, stretch='log', max_percent=99., min_percent=10) norm_co = simple_norm(co_rep, stretch='asinh', max_percent=99.5, min_percent=0.1) norm_hnco = simple_norm(hnco_rep, stretch='asinh', max_percent=99.995, min_percent=0.1) rgb = np.array([norm_hi(lv_hi_mx.value), norm_co(co_rep), norm_hnco(hnco_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]# norm_hi = simple_norm(lv_hi_mx, stretch='log', max_percent=99., min_percent=15) norm_co = simple_norm(co_rep, stretch='asinh', max_percent=99.5, min_percent=0.1) norm_hnco = simple_norm(hnco_rep, stretch='asinh', max_percent=99.995, min_percent=0.1) rgb = np.array([norm_hi(lv_hi_mx.value), norm_co(co_rep), norm_hnco(hnco_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]# 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.5, min_percent=0.1) norm_hnco = simple_norm(hnco_rep, stretch='asinh', max_percent=99.995, min_percent=0.1) rgb = np.array([norm_hi(lv_hi_mx.value), norm_co(co_rep), norm_hnco(hnco_rep), ]).swapaxes(0,2).swapaxes(0,1) 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.5, min_percent=0.1) norm_hnco = simple_norm(hnco_rep, stretch='asinh', max_percent=99.995, min_percent=0.1) rgb = np.array([norm_hi(lv_hi_mx.value), norm_co(co_rep), norm_hnco(hnco_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]# 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.5, min_percent=0.1) norm_hnco = simple_norm(hnco_rep, stretch='asinh', max_percent=99.995, min_percent=0.1) rgb = np.array([norm_hi(lv_hi_mx.value)*0.9, norm_co(co_rep), norm_hnco(hnco_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]# 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.5, min_percent=0.1) norm_hnco = simple_norm(hnco_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), norm_hnco(hnco_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]# 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.5, min_percent=0.1) norm_hnco = simple_norm(hnco_rep, stretch='asinh', max_percent=99.995, min_percent=0.1) rgb = np.array([norm_hi(lv_hi_mx.value)*0.7, norm_hnco(hnco_rep), norm_co(co_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 def make_hsv_fig(img): hsvim = rgb_to_hsv(img) fig = pl.figure(figsize=(16,16), facecolor='k') for ii,angle in enumerate(np.linspace(0,120,9)): pl.subplot(3,3,ii+1) hsvim[:,:,0] = (hsvim[:,:,0]+45/360.) hsvim[:,:,0][hsvim[:,:,0] > 1] -= 1 pl.imshow(hsv_to_rgb(hsvim), origin='lower') pl.xticks([]) pl.yticks([]) pl.tight_layout() return fig make_hsv_fig(rgb) #[Out]#
from matplotlib.colors import hsv_to_rgb, rgb_to_hsv def make_hsv_fig(img): hsvim = rgb_to_hsv(img) fig = pl.figure(figsize=(10,10), facecolor='k') for ii,angle in enumerate(np.linspace(0,120,9)): pl.subplot(3,3,ii+1) hsvim[:,:,0] = (hsvim[:,:,0]+45/360.) hsvim[:,:,0][hsvim[:,:,0] > 1] -= 1 pl.imshow(hsv_to_rgb(hsvim), origin='lower') pl.xticks([]) pl.yticks([]) pl.tight_layout() return fig make_hsv_fig(rgb) #[Out]#
######################################################## # Started Logging At: 2022-02-03 08:34:18 ######################################################## ######################################################## # # Started Logging At: 2022-02-03 08:34:19 ######################################################## 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 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]# 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 # hack to get CO onto HI grid array = lv_dameco_mx.value wcs_in = lv_dameco_mx.wcs wcs_in.wcs.crval[0] = 0 # need to wrap at 0 wcs_out = lv_hi_max.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) 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 # hack to get CO onto HI grid array = lv_dameco_mx.value wcs_in = lv_dameco_mx.wcs 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) 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]# # hack to get CO onto HI grid array = lv_dameco_mx.value wcs_in = lv_dameco_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) 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]# # hack to get CO onto HI grid array = lv_dameco_mx.value wcs_in = lv_dameco_mx.wcs 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]# pixel_in #[Out]# [array([-237491.24137549, -237491.24137549, -237491.24137549, ..., #[Out]# 269042.63798444, 269042.63798444, 269042.63798444]), #[Out]# array([-40.8333324 , -40.75555462, -40.67777685, ..., 40.59999907, #[Out]# 40.67777685, 40.75555462])] lv_dameco_mx.wcs #[Out]# WCS Keywords #[Out]# #[Out]# Number of WCS axes: 2 #[Out]# CTYPE : 'GLON' 'VOPT' #[Out]# CRVAL : 0.0 -319.2394 #[Out]# CRPIX : 1.0 1.0 #[Out]# PC1_1 PC1_2 : 1.0 0.0 #[Out]# PC2_1 PC2_2 : 0.0 1.0 #[Out]# CDELT : -0.125 1.300364 #[Out]# NAXIS : 0 0 # hack to get CO onto HI grid array = lv_dameco_mx.value wcs_in = lv_dameco_mx.wcs #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]# # hack to get CO onto HI grid array = lv_dameco_mx.value wcs_in = lv_dameco_mx.wcs #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 lv_dameco_mx.wcs #[Out]# WCS Keywords #[Out]# #[Out]# Number of WCS axes: 2 #[Out]# CTYPE : 'GLON' 'VOPT' #[Out]# CRVAL : 0.0 -319.2394 #[Out]# CRPIX : 1.0 1.0 #[Out]# PC1_1 PC1_2 : 1.0 0.0 #[Out]# PC2_1 PC2_2 : 0.0 1.0 #[Out]# CDELT : -0.125 1.300364 #[Out]# NAXIS : 0 0 pl.imshow(dameco_rep) #[Out]# lv_dameco_mx.wcs #[Out]# WCS Keywords #[Out]# #[Out]# Number of WCS axes: 2 #[Out]# CTYPE : 'GLON' 'VOPT' #[Out]# CRVAL : 0.0 -319.2394 #[Out]# CRPIX : 1.0 1.0 #[Out]# PC1_1 PC1_2 : 1.0 0.0 #[Out]# PC2_1 PC2_2 : 0.0 1.0 #[Out]# CDELT : -0.125 1.300364 #[Out]# NAXIS : 0 0 # hack to get CO onto HI grid array = lv_dameco_mx.value wcs_in = lv_dameco_mx.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 lv_dameco_mx.wcs #[Out]# WCS Keywords #[Out]# #[Out]# Number of WCS axes: 2 #[Out]# CTYPE : 'GLON' 'VOPT' #[Out]# CRVAL : 0.0 -319239.39999999997 #[Out]# CRPIX : 1.0 1.0 #[Out]# PC1_1 PC1_2 : 1.0 0.0 #[Out]# PC2_1 PC2_2 : 0.0 1.0 #[Out]# CDELT : -0.125 1300.364 #[Out]# NAXIS : 0 0 pl.imshow(dameco_rep) #[Out]# # hack to get CO onto HI grid array = lv_dameco_mx.value wcs_in = lv_dameco_mx.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 # hack to get CO onto HI grid array = lv_dameco_mx.value wcs_in = lv_dameco_mx.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 lv_dameco_mx.wcs #[Out]# WCS Keywords #[Out]# #[Out]# Number of WCS axes: 2 #[Out]# CTYPE : 'GLON' 'VOPT' #[Out]# CRVAL : 0.0 -319239399999.99994 #[Out]# CRPIX : 1.0 1.0 #[Out]# PC1_1 PC1_2 : 1.0 0.0 #[Out]# PC2_1 PC2_2 : 0.0 1.0 #[Out]# CDELT : -0.125 1300364000.0 #[Out]# NAXIS : 0 0 pl.imshow(dameco_rep) #[Out]# # hack to get CO onto HI grid array = lv_dameco_mx.value wcs_in = lv_dameco_mx.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 lv_dameco_mx.wcs #[Out]# WCS Keywords #[Out]# #[Out]# Number of WCS axes: 2 #[Out]# CTYPE : 'GLON' 'VOPT' #[Out]# CRVAL : 0.0 -319239399999999.94 #[Out]# CRPIX : 1.0 1.0 #[Out]# PC1_1 PC1_2 : 1.0 0.0 #[Out]# PC2_1 PC2_2 : 0.0 1.0 #[Out]# CDELT : -0.125 1300364000000.0 #[Out]# NAXIS : 0 0 # 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 lv_dameco_mx.wcs #[Out]# WCS Keywords #[Out]# #[Out]# Number of WCS axes: 2 #[Out]# CTYPE : 'GLON' 'VOPT' #[Out]# CRVAL : 0.0 -319239399999999.94 #[Out]# CRPIX : 1.0 1.0 #[Out]# PC1_1 PC1_2 : 1.0 0.0 #[Out]# PC2_1 PC2_2 : 0.0 1.0 #[Out]# CDELT : -0.125 1300364000000.0 #[Out]# NAXIS : 0 0 wcs_in.wcs #[Out]# flag: 137 #[Out]# naxis: 2 #[Out]# crpix: 0x558aa59a6cd0 #[Out]# 1.0000 1.0000 #[Out]# pc: 0x558b1c2731e0 #[Out]# pc[0][]: 1.0000 0.0000 #[Out]# pc[1][]: 0.0000 1.0000 #[Out]# cdelt: 0x558ab2c68610 #[Out]# -0.12500 1300.4 #[Out]# crval: 0x558ab80aed10 #[Out]# 13.000 -3.1924e+05 #[Out]# cunit: 0x558b1c289e40 #[Out]# "deg" #[Out]# "m s-1" #[Out]# ctype: 0x558ab0e2e350 #[Out]# "GLON" #[Out]# "VOPT" #[Out]# lonpole: 0.000000 #[Out]# latpole: 90.000000 #[Out]# restfrq: 0.000000 #[Out]# restwav: 0.000000 #[Out]# npv: 0 #[Out]# npvmax: 64 #[Out]# pv: 0x558ab5e2d7c0 #[Out]# nps: 0 #[Out]# npsmax: 8 #[Out]# ps: 0x558aa9221380 #[Out]# cd: 0x558b36c014e0 #[Out]# cd[0][]: 0.0000 0.0000 #[Out]# cd[1][]: 0.0000 0.0000 #[Out]# crota: 0x558abdb0f590 #[Out]# 0.0000 0.0000 #[Out]# altlin: 1 #[Out]# velref: 0 #[Out]# alt: ' ' #[Out]# colnum: 0 #[Out]# colax: 0x558aac630ac0 #[Out]# 0 0 #[Out]# cname: 0x558aa5844410 #[Out]# UNDEFINED #[Out]# UNDEFINED #[Out]# crder: 0x558aba133190 #[Out]# UNDEFINED UNDEFINED #[Out]# csyer: 0x558aac5a5cb0 #[Out]# UNDEFINED UNDEFINED #[Out]# czphs: 0x558ab8126630 #[Out]# UNDEFINED UNDEFINED #[Out]# cperi: 0x558aa7f9b590 #[Out]# UNDEFINED UNDEFINED #[Out]# wcsname: UNDEFINED #[Out]# timesys: UNDEFINED #[Out]# trefpos: UNDEFINED #[Out]# trefdir: UNDEFINED #[Out]# plephem: UNDEFINED #[Out]# timeunit: UNDEFINED #[Out]# dateref: "1858-11-17" #[Out]# mjdref: 0.000000000 0.000000000 #[Out]# timeoffs: UNDEFINED #[Out]# dateobs: UNDEFINED #[Out]# datebeg: UNDEFINED #[Out]# dateavg: UNDEFINED #[Out]# dateend: UNDEFINED #[Out]# mjdobs: UNDEFINED #[Out]# mjdbeg: UNDEFINED #[Out]# mjdavg: UNDEFINED #[Out]# mjdend: UNDEFINED #[Out]# jepoch: UNDEFINED #[Out]# bepoch: UNDEFINED #[Out]# tstart: UNDEFINED #[Out]# tstop: UNDEFINED #[Out]# xposure: UNDEFINED #[Out]# telapse: UNDEFINED #[Out]# timsyer: UNDEFINED #[Out]# timrder: UNDEFINED #[Out]# timedel: UNDEFINED #[Out]# timepixr: UNDEFINED #[Out]# obsgeo: UNDEFINED UNDEFINED UNDEFINED #[Out]# UNDEFINED UNDEFINED UNDEFINED #[Out]# obsorbit: UNDEFINED #[Out]# radesys: UNDEFINED #[Out]# equinox: UNDEFINED #[Out]# specsys: "LSRK" #[Out]# ssysobs: UNDEFINED #[Out]# velosys: UNDEFINED #[Out]# zsource: UNDEFINED #[Out]# ssyssrc: UNDEFINED #[Out]# velangl: UNDEFINED #[Out]# aux: 0x0 #[Out]# ntab: 0 #[Out]# tab: 0x0 #[Out]# nwtb: 0 #[Out]# wtb: 0x0 #[Out]# types: 0x558aa56e65b0 #[Out]# 2000 3000 #[Out]# lngtyp: "GLON" #[Out]# lattyp: " " #[Out]# lng: 0 #[Out]# lat: -1 #[Out]# spec: 1 #[Out]# cubeface: -1 #[Out]# err: 0x0 #[Out]# lin: (see below) #[Out]# cel: (see below) #[Out]# spc: (see below) #[Out]# m_flag: 137 #[Out]# m_naxis: 2 #[Out]# m_crpix: 0x558aa59a6cd0 (= crpix) #[Out]# m_pc: 0x558b1c2731e0 (= pc) #[Out]# m_cdelt: 0x558ab2c68610 (= cdelt) #[Out]# m_crval: 0x558ab80aed10 (= crval) #[Out]# m_cunit: 0x558b1c289e40 (= cunit) #[Out]# m_ctype: 0x558ab0e2e350 (= ctype) #[Out]# m_pv: 0x558ab5e2d7c0 (= pv) #[Out]# m_ps: 0x558aa9221380 (= ps) #[Out]# m_cd: 0x558b36c014e0 (= cd) #[Out]# m_crota: 0x558abdb0f590 (= crota) #[Out]# #[Out]# m_colax: 0x558aac630ac0 (= colax) #[Out]# m_cname: 0x558aa5844410 (= cname) #[Out]# m_crder: 0x558aba133190 (= crder) #[Out]# m_csyer: 0x558aac5a5cb0 (= csyer) #[Out]# m_czphs: 0x558ab8126630 (= czphs) #[Out]# m_cperi: 0x558aa7f9b590 (= cperi) #[Out]# m_aux: 0x0 (= aux) #[Out]# m_tab: 0x0 (= tab) #[Out]# m_wtb: 0x0 (= wtb) #[Out]# #[Out]# lin.* #[Out]# flag: 137 #[Out]# naxis: 2 #[Out]# crpix: 0x558aa59a6cd0 #[Out]# 1.0000 1.0000 #[Out]# pc: 0x558b1c2731e0 #[Out]# pc[0][]: 1.0000 0.0000 #[Out]# pc[1][]: 0.0000 1.0000 #[Out]# cdelt: 0x558ab2c68610 #[Out]# -0.12500 1300.4 #[Out]# dispre: 0x0 #[Out]# disseq: 0x0 #[Out]# piximg: (nil) #[Out]# imgpix: (nil) #[Out]# i_naxis: 0 #[Out]# unity: 1 #[Out]# affine: 1 #[Out]# simple: 1 #[Out]# err: 0x0 #[Out]# tmpcrd: 0x558aa57037e0 #[Out]# m_flag: 0 #[Out]# m_naxis: 0 #[Out]# m_crpix: 0x0 #[Out]# m_pc: 0x0 #[Out]# m_cdelt: 0x0 #[Out]# m_dispre: 0x0 #[Out]# m_disseq: 0x0 #[Out]# #[Out]# cel.* #[Out]# flag: 0 #[Out]# offset: 0 #[Out]# phi0: 0.000000 #[Out]# theta0: 0.000000 #[Out]# ref: 0.0000 0.0000 9.8765e+107 90.000 #[Out]# prj: (see below) #[Out]# euler: 0.0000 0.0000 0.0000 0.0000 0.0000 #[Out]# latpreq: -1 (UNDEFINED) #[Out]# isolat: 0 #[Out]# err: 0x0 #[Out]# #[Out]# prj.* #[Out]# flag: 0 #[Out]# code: " " #[Out]# r0: 0.000000 #[Out]# pv: (not used) #[Out]# phi0: UNDEFINED #[Out]# theta0: UNDEFINED #[Out]# bounds: 7 #[Out]# #[Out]# name: "undefined" #[Out]# category: 0 (undefined) #[Out]# pvrange: 0 #[Out]# simplezen: 0 #[Out]# equiareal: 0 #[Out]# conformal: 0 #[Out]# global: 0 #[Out]# divergent: 0 #[Out]# x0: 0.000000 #[Out]# y0: 0.000000 #[Out]# err: 0x0 #[Out]# w[]: 0.0000 0.0000 0.0000 0.0000 0.0000 #[Out]# 0.0000 0.0000 0.0000 0.0000 0.0000 #[Out]# m: 0 #[Out]# n: 0 #[Out]# prjx2s: 0x0 #[Out]# prjs2x: 0x0 #[Out]# #[Out]# spc.* #[Out]# flag: 0 #[Out]# type: " " #[Out]# code: " " #[Out]# crval: UNDEFINED #[Out]# restfrq: 0.000000 #[Out]# restwav: 0.000000 #[Out]# pv: (not used) #[Out]# w: 0.0000 0.0000 0.0000 (remainder unused) #[Out]# isGrism: 0 #[Out]# err: 0x0 #[Out]# spxX2P: 0x0 #[Out]# spxP2S: 0x0 #[Out]# spxS2P: 0x0 #[Out]# spxP2X: 0x0 wcs_in #[Out]# WCS Keywords #[Out]# #[Out]# Number of WCS axes: 2 #[Out]# CTYPE : 'GLON' 'VOPT' #[Out]# CRVAL : 13.0 -319239.39999999997 #[Out]# CRPIX : 1.0 1.0 #[Out]# PC1_1 PC1_2 : 1.0 0.0 #[Out]# PC2_1 PC2_2 : 0.0 1.0 #[Out]# CDELT : -0.125 1300.364 #[Out]# NAXIS : 0 0 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.5, 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), norm_hnco(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]# 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.5, 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), 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]# 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), norm_dameco(dameco_rep), ]).swapaxes(0,2).swapaxes(0,1) 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), 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]# 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 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), 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]# pl.figure(figsize=(16,10)) slc = [slice(50,750), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(rgb[slc[0],slc[1],:]) #[Out]# pl.figure(figsize=(16,10)) slc = [slice(100,750), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(rgb[slc[0],slc[1],:]) #[Out]# pl.figure(figsize=(16,10)) slc = [slice(50,750), slice(100,900)] ax = pl.subplot(projection=lv_hi_mx.wcs[slc]) ax.imshow(rgb[slc[0],slc[1],:]) #[Out]# 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]# 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.9, norm_co(co_rep), 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]# 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), 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]# 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]# 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]# 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]# from matplotlib.colors import hsv_to_rgb, rgb_to_hsv hsvim = rgb_to_hsv(rgb) hsvim[:,:,0] = (hsvim[:,:,0]+15/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]# from matplotlib.colors import hsv_to_rgb, rgb_to_hsv hsvim = rgb_to_hsv(rgb) hsvim[:,:,0] = (hsvim[:,:,0]+-15/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]# from matplotlib.colors import hsv_to_rgb, rgb_to_hsv hsvim = rgb_to_hsv(rgb) hsvim[:,:,0] = (hsvim[:,:,0]-15/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]# 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.) pl.imshow(hsv_to_rgb(hsv_hi)) #[Out]# pl.imshow(rotate_color(norm_co(co_rep)), -30) def rotate_color(img, angle): 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_co(co_rep)), -30) pl.imshow(rotate_color(norm_co(co_rep), -30)) #[Out]# pl.imshow(rotate_color(norm_co(co_rep), 30)) #[Out]# pl.imshow(rotate_color(norm_co(co_rep), -60)) #[Out]# pl.imshow(rotate_color(norm_co(co_rep), 90)) #[Out]# pl.imshow(rotate_color(norm_co(co_rep), 150)) #[Out]# pl.imshow(rotate_color(norm_co(co_rep), 180)) #[Out]# pl.imshow(rotate_color(norm_co(co_rep), 270)) #[Out]# pl.imshow(rotate_color(norm_dameco(dameco_rep), 200)) #[Out]# pl.imshow(rotate_color(norm_hi(lv_hi_mx.value)*0.7, 45)) #[Out]# pl.imshow(rotate_color(norm_dameco(dameco_rep), 200) + rotate_color(norm_co(co_rep), 270) + rotate_color(norm_hi(lv_hi_mx.value)*0.7, 45)) #[Out]# pl.imshow(rotate_color(norm_dameco(dameco_rep), 200)*0.8 + rotate_color(norm_co(co_rep), 270) + rotate_color(norm_hi(lv_hi_mx.value)*0.7, 45)) #[Out]# new_rgb = rotate_color(norm_dameco(dameco_rep), 200)*0.8 + rotate_color(norm_co(co_rep), 270) + rotate_color(norm_hi(lv_hi_mx.value)*0.7, 45) pl.imshow(new_rgb[slc[0], slc[1], :]) #[Out]# new_rgb = rotate_color(norm_dameco(dameco_rep), 200)*0.8 + rotate_color(norm_co(co_rep), 270) + rotate_color(norm_hi(lv_hi_mx.value)*0.7, 45) pl.imshow(new_rgb[slc[0], slc[1], :], origin='lower') #[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.imshow(new_rgb[slc[0], slc[1], :], origin='lower') #[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=(16,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') #[Out]# 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]# 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]# 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=(16,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') #[Out]#