diff --git a/tutorials/spectral-cube-reprojection/SpectralCubeReprojectExample.ipynb b/tutorials/spectral-cube-reprojection/SpectralCubeReprojectExample.ipynb new file mode 100644 index 00000000..8def15ee --- /dev/null +++ b/tutorials/spectral-cube-reprojection/SpectralCubeReprojectExample.ipynb @@ -0,0 +1,646 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b6372d7e", + "metadata": {}, + "source": [ + "# Cube Reprojection Tutorial\n", + "\n", + "## Authors\n", + "Adam Ginsburg, Eric Koch\n", + "\n", + "## Learning Goals\n", + "* reproject a cube spectrally\n", + "* smooth it spectrally\n", + "* reproject it spatially\n", + "\n", + "## Keywords\n", + "spectral cube, radio astronomy, astroquery, units, dask\n", + "\n", + "## Summary\n", + "Spectroscopic cube observations taken at different wavelength can trace the motion of gas or stars using spectral lines, but often lines at different wavelengths give different information.\n", + "For example, one might observe a galaxy in the 21cm line of HI and the 115 GHz line of CO, or a protoplanetary disk in a line of N2H+ and a line of CO, or a galactic disk in the H-alpha and H-beta lines (in absorption or emission).\n", + "In order to compare these data sets pixel-by-pixel, they must be placed onto a common grid with common resolution.\n", + "\n", + "This tutorial shows how to take two spectral cubes observed toward the same part of the sky, but different frequencies, and put them onto the same grid using [spectral-cube](spectral-cube.readthedocs.io).\n", + "\n", + "It uses [astroquery](https://astroquery.readthedocs.io/) to obtain line frequencies from [splatalogue](https://astroquery.readthedocs.io/en/latest/splatalogue/splatalogue.html); this example uses radio-wavelength data for which Splatalogue's molecular line lists are appropriate.\n", + "Finally, it shows how to do the reprojection using [dask](https://dask.org) to enable parallelization." + ] + }, + { + "cell_type": "markdown", + "id": "70381c49", + "metadata": {}, + "source": [ + "## Index \n", + "\n", + " * [Step 1: Download](#Step-1:-Download-the-data)\n", + " * [Step 2: Open files, collect metadata](#Step-2:-Load-the-cubes)\n", + " * [Step 3: Convert to velocity](#Step-3:-Convert-cubes-from-frequency-to-velocity)\n", + " * [Step 4: Spectral Interpolation](#Step-4.-Spectral-Interpolation)\n", + " * [Step 5: Spatial Smoothing](#Step-5.-Spatial-Smoothing)\n", + " * [Step 6: Reprojection](#Step-6.-Reprojection)\n", + " \n", + " \n", + "In this example, we do spectral smoothing and interpolation (step 4) before spatial smoothing and interpolation (step 5), but if you have a varying-resolution cube (with a different beam size for each channel), you have to do spatial smoothing first. For more information see the [spectral-cube documentation](spectral-cube.readthedocs.io)." + ] + }, + { + "cell_type": "markdown", + "id": "48448b34", + "metadata": {}, + "source": [ + "## Step 1: Download the data\n", + "\n", + "(you might not have to do this step, since you may already have data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "32902371", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from astropy.utils.data import download_file" + ] + }, + { + "cell_type": "markdown", + "id": "4dec821a", + "metadata": {}, + "source": [ + "We download two example spectral cubes of a point in the Galactic center from a permalink on the ALMA archives.\n", + "These are moderately large files, with sizes 18 MB and 337 MB.\n", + "\n", + "If you have trouble with these downloads, try changing to a different ALMA server (e.g., almascience.eso.org->almascience.nrao.edu) or increase the timeout. See the [download_file](https://docs.astropy.org/en/stable/api/astropy.utils.data.download_file.html) documentation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9719f71a", + "metadata": {}, + "outputs": [], + "source": [ + "filename_1 = download_file(\"https://almascience.eso.org/dataPortal/member.uid___A001_X1465_X3a33.BrickMaser_sci.spw71.cube.I.manual.image.pbcor.fits\",\n", + " cache=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2a1d4dc7", + "metadata": {}, + "outputs": [], + "source": [ + "filename_2 = download_file(\"https://almascience.eso.org/dataPortal/member.uid___A001_X87d_X141.a_sma1_sci.spw27.cube.I.pbcor.fits\",\n", + " cache=True)" + ] + }, + { + "cell_type": "markdown", + "id": "d801b283", + "metadata": {}, + "source": [ + "## Step 2: Load the cubes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe4ead93", + "metadata": {}, + "outputs": [], + "source": [ + "from spectral_cube import SpectralCube" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "833abd53", + "metadata": {}, + "outputs": [], + "source": [ + "cube1 = SpectralCube.read(filename_1)\n", + "cube1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "04b8ddcd", + "metadata": {}, + "outputs": [], + "source": [ + "cube2 = SpectralCube.read(filename_2)\n", + "cube2" + ] + }, + { + "cell_type": "markdown", + "id": "19595c31", + "metadata": {}, + "source": [ + "The cubes are at different frequencies - 139 and 217 GHz.\n", + "\n", + "The first cube covers the H2CS 4(1,3)-3(1,2) line at 139.483699\tGHz.\n", + "\n", + "The second covers SiO v=5-4 at 217.104984 GHz" + ] + }, + { + "cell_type": "markdown", + "id": "ccc5334f", + "metadata": {}, + "source": [ + "We use the `find_lines` tool to query [splatalogue](https://splatalogue.online/) with [astroquery](https://astroquery.readthedocs.io/en/latest/splatalogue/splatalogue.html) over the spectral range covered by the cube. It returns a table of matching lines. Note that some line names will be repeated because Splatalogue includes several different databases and most chemical species are present in all of these." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fdf8959f", + "metadata": {}, + "outputs": [], + "source": [ + "cube1.find_lines(chemical_name=' H2CS ').show_in_notebook()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ba4fbb6a", + "metadata": {}, + "outputs": [], + "source": [ + "cube2.find_lines(chemical_name='SiO').show_in_notebook()" + ] + }, + { + "cell_type": "markdown", + "id": "8c7234cb", + "metadata": {}, + "source": [ + "## Step 3: Convert cubes from frequency to velocity\n", + "\n", + "To compare the kinematic structure of the target, we need to convert from the observed frequency (which must be in a common reference frame; in this case, it already is) to the doppler velocity." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9b59615f", + "metadata": {}, + "outputs": [], + "source": [ + "from astropy import units as u" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "86aa822b", + "metadata": {}, + "outputs": [], + "source": [ + "cube1vel = cube1.with_spectral_unit(u.km/u.s, velocity_convention='radio', rest_value=139.483699*u.GHz)\n", + "cube1vel" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a1544293", + "metadata": {}, + "outputs": [], + "source": [ + "cube2vel = cube2.with_spectral_unit(u.km/u.s, velocity_convention='radio', rest_value=217.104984*u.GHz)\n", + "cube2vel" + ] + }, + { + "cell_type": "markdown", + "id": "11e7873f", + "metadata": {}, + "source": [ + "From the shape of the cube, we can see the H2CS cube is narrower in velocity, so we'll use that as the target spectral reprojection. However, the SiO cube is the smaller footprint on the sky." + ] + }, + { + "cell_type": "markdown", + "id": "cc6f783c", + "metadata": {}, + "source": [ + "### Create spatial maps of the peak intensity to quickly explore the cubes:\n", + " \n", + "One way to quickly explore the structure in the data cubes is to produce a peak intensity map, or the maximum along the spectral axis (`axis=0`)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bb078eb0", + "metadata": {}, + "outputs": [], + "source": [ + "mx = cube1.max(axis=0)\n", + "mx.quicklook()" + ] + }, + { + "cell_type": "markdown", + "id": "b8e96d66", + "metadata": {}, + "source": [ + "We can do the same thing all on one line (for the other cube this time):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8b035a79", + "metadata": {}, + "outputs": [], + "source": [ + "cube2.max(axis=0).quicklook()" + ] + }, + { + "cell_type": "markdown", + "id": "98ecf111", + "metadata": {}, + "source": [ + "# Step 4. Spectral Interpolation\n", + "\n", + "We can choose to do either the spatial or spectral step first. \n", + "In this case, we choose the spectral step first because the H$_2$CS cube is narrower in velocity (`cube1vel`) and this will reduce the number of channels we need to spatially interpolate over in the next step.\n", + "\n", + "We need to match resolution to the cube with the largest channel width:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e19188f0", + "metadata": {}, + "outputs": [], + "source": [ + "velocity_res_1 = np.diff(cube1vel.spectral_axis)[0]\n", + "velocity_res_2 = np.diff(cube2vel.spectral_axis)[0]\n", + "velocity_res_1, velocity_res_2" + ] + }, + { + "cell_type": "markdown", + "id": "007fac49", + "metadata": {}, + "source": [ + "Next, we will reduce `cube2vel` to have the same spectral range as `cube1vel`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cee190fc", + "metadata": {}, + "outputs": [], + "source": [ + "cube2vel_cutout = cube2vel.spectral_slab(cube1vel.spectral_axis.min(),\n", + " cube1vel.spectral_axis.max())\n", + "cube1vel, cube2vel_cutout" + ] + }, + { + "cell_type": "markdown", + "id": "385b7042", + "metadata": {}, + "source": [ + "Note that it is important for the to-be-interpolated cube, in this case `cube2`, to have pixels bounding `cube1`'s spectral axis, but in this case it does not. If the pixel range doesn't overlap perfectly, it may blank out one of the edge pixels. So, to fix this, we add a little buffer:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "653d7199", + "metadata": {}, + "outputs": [], + "source": [ + "cube2vel_cutout = cube2vel.spectral_slab(cube1vel.spectral_axis.min() - velocity_res_2,\n", + " cube1vel.spectral_axis.max())\n", + "cube1vel, cube2vel_cutout" + ] + }, + { + "cell_type": "markdown", + "id": "c5a0e750", + "metadata": {}, + "source": [ + "Our H2CS cube (`cube1vel`) has broader channels. We need to first smooth `cube2vel` to the broader channel width before doing the spatial reprojection.\n", + "\n", + "To do this, we will spectrally smooth with a Gaussian with width set such that smoothing `cube2vel` will result in the same width as `cube1vel`. We do this by finding the difference in widths when deconvolving the `cube1vel` channel width from `cube2vel`. For further information see the [documentation on smoothing](https://spectral-cube.readthedocs.io/en/latest/smoothing.html#spectral-smoothing).\n", + "\n", + "Note that if we did not do this smoothing step, we would under-sample the `cube2vel` data in the next downsampling step, reducing our signal-to-noise ratio.\n", + "\n", + "We have adopted a width equal to the channel width; the [line spread function](https://help.almascience.org/kb/articles/what-spectral-resolution-will-i-get-for-a-given-channel-spacing) is actually a Hanning-smoothed tophat. We are making a coarse approximation here." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c86d7716", + "metadata": {}, + "outputs": [], + "source": [ + "fwhm_gaussian = (velocity_res_1**2 - velocity_res_2**2)**0.5\n", + "fwhm_gaussian" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "35acafb7", + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.convolution import Gaussian1DKernel\n", + "fwhm_to_sigma = np.sqrt(8*np.log(2))\n", + "# we want the kernel in pixel units, so we force to km/s and take the value\n", + "spectral_smoothing_kernel = Gaussian1DKernel(stddev=fwhm_gaussian.to(u.km/u.s).value / fwhm_to_sigma)" + ] + }, + { + "cell_type": "markdown", + "id": "ffdd178a", + "metadata": {}, + "source": [ + "We then smooth with the kernel. Note that this is doing 420x420 = 176400 smoothing operations on a length-221 spectrum: it will take a little time" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "86ada2d7", + "metadata": {}, + "outputs": [], + "source": [ + "cube2vel_smooth = cube2vel_cutout.spectral_smooth(spectral_smoothing_kernel)" + ] + }, + { + "cell_type": "markdown", + "id": "00f6d313", + "metadata": {}, + "source": [ + "Now that we've done spectral smoothing, we can resample the spectral axis of `cube2vel_smooth` to match `cube1vel` by interpolating `cube2vel_smooth` onto `cube1vel`'s grid:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7655a346", + "metadata": {}, + "outputs": [], + "source": [ + "cube2vel_spectralresample = cube2vel_smooth.spectral_interpolate(cube1vel.spectral_axis,\n", + " suppress_smooth_warning=True)\n", + "cube2vel_spectralresample" + ] + }, + { + "cell_type": "markdown", + "id": "2eda47cf", + "metadata": {}, + "source": [ + "Note that we included the `suppress_smooth_warning=True` argument. That is to hide this warning:\n", + "```\n", + "WARNING: SmoothingWarning: Input grid has too small a spacing. The data should be smoothed prior to resampling. [spectral_cube.spectral_cube]\n", + "```\n", + "which will tell you if the operation will under-sample the original data. The smoothing work we did above is specifically to make sure we are properly sampling, so this warning does not apply." + ] + }, + { + "cell_type": "markdown", + "id": "00926359", + "metadata": {}, + "source": [ + "# Step 5. Spatial Smoothing" + ] + }, + { + "cell_type": "markdown", + "id": "14e99bbb", + "metadata": {}, + "source": [ + "Now that we've done spectral smoothing, we also need to follow a similar procedure of smoothing then resampling for the spatial axes. \n", + "\n", + "The `beam` is the resolution element of our cubes:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4f0512f4", + "metadata": {}, + "outputs": [], + "source": [ + "cube1vel.beam, cube2vel_spectralresample.beam" + ] + }, + { + "cell_type": "markdown", + "id": "c93654ca", + "metadata": {}, + "source": [ + "`cube1` again hase the larger beam, so we'll smooth `cube2` to its resolution" + ] + }, + { + "cell_type": "markdown", + "id": "e6c4d626", + "metadata": {}, + "source": [ + "#### Aside: mixed beams \n", + "\n", + "If cube1 and cube2 had different sized beams, but neither was clearly larger, we would have to convolve _both_ to a [common beam](https://radio-beam.readthedocs.io/en/latest/commonbeam.html#finding-the-smallest-common-beam).\n", + "\n", + "In this case, it's redundant and we could have just used `cube1`'s beam, but this is the more general approach:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b3cd31a7", + "metadata": {}, + "outputs": [], + "source": [ + "import radio_beam\n", + "common_beam = radio_beam.commonbeam.common_2beams(radio_beam.Beams(beams=[cube1vel.beam, cube2vel.beam]))\n", + "common_beam" + ] + }, + { + "cell_type": "markdown", + "id": "74ff5976", + "metadata": {}, + "source": [ + "We then convolve:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "eaa9f644", + "metadata": {}, + "outputs": [], + "source": [ + "# for v<0.6, we convert to Kelvin to ensure the units are preserved:\n", + "# cube2vel_spatialspectralsmooth = cube2vel_spectralresample.to(u.K).convolve_to(common_beam)\n", + "# in more recent versions, the unit conversion is handled appropriately,\n", + "# so unit conversion isn't needed\n", + "cube2vel_spatialspectralsmooth = cube2vel_spectralresample.convolve_to(common_beam)\n", + "cube2vel_spatialspectralsmooth" + ] + }, + { + "cell_type": "markdown", + "id": "cdc73378", + "metadata": {}, + "source": [ + "# Step 6. Reprojection\n", + "\n", + "Now we can do the spatial resampling as the final step for producing two cubes matched to the same spatial and spectral pixel grid:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "de8b8a48", + "metadata": {}, + "outputs": [], + "source": [ + "cube2vel_reproj = cube2vel_spatialspectralsmooth.reproject(cube1vel.header)\n", + "cube2vel_reproj" + ] + }, + { + "cell_type": "markdown", + "id": "69d8a90b", + "metadata": {}, + "source": [ + "These two cubes are now on an identical grid, and can be directly compared:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5bd3de8d", + "metadata": {}, + "outputs": [], + "source": [ + "cube2vel_reproj, cube1vel" + ] + }, + { + "cell_type": "markdown", + "id": "8cf176f3", + "metadata": {}, + "source": [ + "These spectra can now be overplotted as they are in the same unit with the same beam." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a196eb8b", + "metadata": {}, + "outputs": [], + "source": [ + "cube1vel[:,125,125].quicklook()\n", + "cube2vel_reproj[:,125,125].quicklook()" + ] + }, + { + "cell_type": "markdown", + "id": "13afd6f7", + "metadata": {}, + "source": [ + "# Dask\n", + "\n", + "All of the above can be done using `dask` as the underlying framework to parallelize the operations.\n", + "See [the spectral-cube documentation on dask integration](https://spectral-cube.readthedocs.io/en/latest/dask.html) or the [dask documentation](https://dask.org/) for further details.\n", + "\n", + "The dask approach can be made more memory-efficient (avoid using too much RAM) by writing intermediate steps to disk. The non-dask approach used above will generally need to read the whole cube into memory. Depending on the situation, either approach may be faster, but `dask` may be needed if the cube is larger than memory." + ] + }, + { + "cell_type": "markdown", + "id": "0934aa4b", + "metadata": {}, + "source": [ + "We repeat all the operations above using dask. We use a `ProgressBar` so you can see how long it takes. We also suppress warnings to make the output look cleaner (we already saw all the important warnings above)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "03e1a267", + "metadata": {}, + "outputs": [], + "source": [ + "from dask.diagnostics import ProgressBar\n", + "import warnings" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14d47820", + "metadata": {}, + "outputs": [], + "source": [ + "with warnings.catch_warnings():\n", + " warnings.simplefilter('ignore')\n", + " with ProgressBar():\n", + " cube2dask = SpectralCube.read(filename_2, use_dask=True)\n", + " cube2daskvel = cube2dask.with_spectral_unit(u.km/u.s,\n", + " velocity_convention='radio', rest_value=217.104984*u.GHz)\n", + " cube2daskvel_cutout = cube2daskvel.spectral_slab(cube1vel.spectral_axis.min() - velocity_res_2,\n", + " cube1vel.spectral_axis.max())\n", + " cube2daskvel_smooth = cube2daskvel_cutout.spectral_smooth(spectral_smoothing_kernel)\n", + " cube2daskvel_spectralresample = cube2daskvel_smooth.spectral_interpolate(cube1vel.spectral_axis,\n", + " suppress_smooth_warning=True)\n", + " cube2daskvel_spatialspectralsmooth = cube2daskvel_spectralresample.convolve_to(common_beam)\n", + " cube2daskvel_reproj = cube2daskvel_spatialspectralsmooth.reproject(cube1vel.header)\n", + "cube2daskvel_reproj" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4af66ead", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}