Working with FITS-cubes

Working with FITS-cubes#

Authors#

Dhanesh Krishnarao (DK), Shravan Shetty, Diego Gonzalez-Casanova, Audra Hernandez, Kris Stern, Kelle Cruz, Stephanie Douglas

Learning Goals#

  • Find and download data using astroquery

  • Read and plot slices across different dimensions of a data cube

  • Compare different data sets (2D and 3D) by overploting contours

  • Transform coordinate projections and match data resolutions with reproject

  • Create intensity moment maps / velocity maps with spectral_cube

Keywords#

FITS, image manipulation, data cubes, radio astronomy, WCS, astroquery, reproject, spectral cube, matplotlib, contour plots, colorbar

Summary#

In this tutorial we will visualize 2D and 3D data sets in Galactic and equatorial coordinates.

The tutorial will walk you though a visual analysis of the Small Magellanic Cloud (SMC) using HI 21cm emission and a Herschel 250 micron map. We will learn how to read in data from a file, query and download matching data from Herschel using astroquery, and plot the resulting images in a multitude of ways.

The primary libraries we’ll be using are: astroquery, spectral_cube, reproject, matplotlib)

They can be installed using conda:

conda install -c conda-forge astroquery
conda install -c conda-forge spectral-cube
conda install -c conda-forge reproject

Alternatively, if you don’t use conda, you can use pip.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

import astropy.units as u
from astropy.utils.data import download_file
from astropy.io import fits  # We use fits to open the actual data file

from astropy.utils import data

from spectral_cube import SpectralCube

from astroquery.esasky import ESASky
from astroquery.utils import TableList
from astropy.wcs import WCS
from reproject import reproject_interp

data.conf.remote_timeout = 60
/opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Download the HI Data#

We’ll be using HI 21 cm emission data from the HI4Pi survey. We want to look at neutral gas emission from the Magellanic Clouds and learn about the kinematics of the system and column densities. Using the VizieR catalog, we’ve found a relevant data cube to use that covers this region of the sky. You can also download an allsky data cube, but this is a very large file, so picking out sub-sections can be useful!

For us, the relevant file is available via ftp from CDS Strasbourg. We have a reduced version of it which will be a FITS data cube in Galactic coordinates using the tangential sky projection.

Sure, we could download this file directly, but why do that when we can load it up via one line of code and have it ready to use in our cache?

Download the HI Fits Cube#

# Downloads the HI data in a fits file format
hi_datafile = download_file(
    "http://data.astropy.org/tutorials/FITS-cubes/reduced_TAN_C14.fits",
    cache=True,
    show_progress=True,
)

Awesome, so now we have a copy of the data file (a FITS file). So how do we do anything with it?

Luckily for us, the spectral_cube package does a lot of the nitty gritty work for us to manipulate this data and even quickly look through it. So let’s open up our data file and read in the data as a SpectralCube!

The variable cube has the data using SpectralCube and hi_data is the data cube from the FITS file without the special formating from SpectralCube.

hi_data = fits.open(hi_datafile)  # Open the FITS file for reading
cube = SpectralCube.read(hi_data)  # Initiate a SpectralCube
hi_data.close()  # Close the FITS file - we already read it in and don't need it anymore!
If you happen to already have the FITS file on your system, you can also skip the fits.open step and just directly read a FITS file with SpectralCube like this:

#cube = SpectralCube.read('path_to_data_file/TAN_C14.fits')

So what does this SpectralCube object actually look like? Let’s find out! The first check is to print out the cube.

print(cube)
SpectralCube with shape=(450, 150, 150) and unit=K:
 n_x:    150  type_x: GLON-TAN  unit_x: deg    range:   286.727203 deg:  320.797623 deg
 n_y:    150  type_y: GLAT-TAN  unit_y: deg    range:   -50.336450 deg:  -28.401234 deg
 n_s:    450  type_s: VRAD      unit_s: m / s  range:  -598824.534 m / s:  600409.133 m / s

Some things to pay attention to here:#

A data cube has three axes. In this case, there is Galactic Longitude (x), Galactic Latitude (y), and a spectral axis in terms of a LSR Velocity (z - listed as s with spectral_cube).

The data hidden in the cube lives as an ndarray with shape (n_s, n_y, n_x) so that axis 0 corresponds with the Spectral Axis, axis 1 corresponds with the Galactic Latitude Axis, and axis 2 corresponds with the Galactic Longitude Axis.

When we print(cube) we can see the shape, size, and units of all axes as well as the data stored in the cube. With this cube, the units of the data in the cube are temperatures (K). The spatial axes are in degrees and the Spectral Axis is in (meters / second).

The cube also contains information about the coordinates corresponding to the data in the form of a WCS (World Coordinate System) object.

SpectralCube is clever and keeps all the data masked until you really need it so that you can work with large sets of data. So let’s see what our data actually looks like!

SpectralCube has a quicklook() method which can give a handy sneak-peek preview of the data. It’s useful when you just need to glance at a slice or spectrum without knowing any other information (say, to make sure the data isn’t corrupted or is looking at the right region.)

To do this, we index our cube along one axis (for a slice) or two axes (for a spectrum):

cube[
    300, :, :
].quicklook()  # Slice the cube along the spectral axis, and display a quick image
_images/9975fe7596148d50097c14981a812a829aedb2254de77f112d76982489f11327.png
cube[:, 75, 75].quicklook()  # Extract a single spectrum through the data cube
_images/51c97ca462963e6399451d6411b343a844dc2dd3d68ce38af859c6a1f7fcc456.png

Try messing around with slicing the cube along different axes, or picking out different spectra#

Make a smaller cube, focusing on the Magellanic Clouds#

The HI data cube we downloaded is bigger than we actually need it to be. Let’s try zooming in on just the part we need and make a new sub_cube.

The easiest way to do this is to cut out part of the cube with indices or coordinates.

We can extract the world coordinates from the cube using the .world() method.

Warning: using .world() will extract coordinates from every position you ask for. This can be a TON of data if you don't slice through the cube. One work around is to slice along two axes and extract coordinates just along a single dimension.

The output of .world() is an Astropy Quantity representing the pixel coordinates, which includes units. You can extract these Astropy Quantity objects by slicing the data.

_, lat, _ = cube.world[0, :, 0]  # extract latitude world coordinates from cube
_, _, lon = cube.world[0, 0, :]  # extract longitude world coordinates from cube

You can then extract a sub_cube in the spatial coordinates of the cube

# Define desired latitude and longitude range
lat_range = [-46, -40] * u.deg
lon_range = [306, 295] * u.deg

# Create a sub_cube cut to these coordinates
sub_cube = cube.subcube(
    xlo=lon_range[0], xhi=lon_range[1], ylo=lat_range[0], yhi=lat_range[1]
)

print(sub_cube)
SpectralCube with shape=(450, 38, 57) and unit=K:
 n_x:     57  type_x: GLON-TAN  unit_x: deg    range:   294.113498 deg:  306.009028 deg
 n_y:     38  type_y: GLAT-TAN  unit_y: deg    range:   -46.014280 deg:  -40.027398 deg
 n_s:    450  type_s: VRAD      unit_s: m / s  range:  -598824.534 m / s:  600409.133 m / s

Cut along the Spectral Axis:#

We don’t really need data from such a large velocity range so let’s just extract a little slab. We can do this in any units that we want using the .spectral_slab() method.

sub_cube_slab = sub_cube.spectral_slab(-300.0 * u.km / u.s, 300.0 * u.km / u.s)

print(sub_cube_slab)
SpectralCube with shape=(226, 38, 57) and unit=K:
 n_x:     57  type_x: GLON-TAN  unit_x: deg    range:   294.113498 deg:  306.009028 deg
 n_y:     38  type_y: GLAT-TAN  unit_y: deg    range:   -46.014280 deg:  -40.027398 deg
 n_s:    226  type_s: VRAD      unit_s: m / s  range:  -299683.842 m / s:  301268.441 m / s

Moment Maps#

Moment maps are a useful analysis tool to study data cubes. In short, a moment is a weighted integral along an axis (typically the Spectral Axis) that can give information about the total Intensity (or column density), mean velocity, or velocity dispersion along lines of sight.

SpectralCube makes this very simple with the .moment() method. We can convert to friendlier spectral units of km/s and these new 2D projections can be saved as new FITS files, complete with modified WCS information as well.

moment_0 = sub_cube_slab.with_spectral_unit(u.km / u.s).moment(
    order=0
)  # Zero-th moment
moment_1 = sub_cube_slab.with_spectral_unit(u.km / u.s).moment(order=1)  # First moment

# Write the moments as a FITS image
# moment_0.write('hi_moment_0.fits')
# moment_1.write('hi_moment_1.fits')

print("Moment_0 has units of: ", moment_0.unit)
print("Moment_1 has units of: ", moment_1.unit)

# Convert Moment_0 to a Column Density assuming optically thin media
hi_column_density = moment_0 * 1.82 * 10**18 / (u.cm * u.cm) * u.s / u.K / u.km
Moment_0 has units of:  K km / s
Moment_1 has units of:  km / s

Display the Moment Maps#

The WCSAxes framework in Astropy allows us to display images with different coordinate axes and projections.

As long as we have a WCS object associated with the data, we can transfer that projection to a matplotlib axis. SpectralCube makes it possible to access just the WCS object associated with a cube object.

print(moment_1.wcs)  # Examine the WCS object associated with the moment map
WCS Keywords

Number of WCS axes: 2
CTYPE : 'GLON-TAN' 'GLAT-TAN' 
CRVAL : 303.75 -40.0 
CRPIX : 11.557620817843869 40.3723930634757 
PC1_1 PC1_2  : 1.0 0.0 
PC2_1 PC2_2  : 0.0 1.0 
CDELT : -0.1494444443846667 0.1538888888273333 
NAXIS : 0  0

As expected, the first moment image we created only has two axes (Galactic Longitude and Galactic Latitude). We can pass in this WCS object directly into a matplotlib axis instance.

# Initiate a figure and axis object with WCS projection information
fig = plt.figure(figsize=(18, 12))
ax = fig.add_subplot(111, projection=moment_1.wcs)

# Display the moment map image
im = ax.imshow(moment_1.hdu.data, cmap="RdBu_r", vmin=0, vmax=200)
ax.invert_yaxis()  # Flips the Y axis

# Add axes labels
ax.set_xlabel("Galactic Longitude (degrees)", fontsize=16)
ax.set_ylabel("Galactic Latitude (degrees)", fontsize=16)

# Add a colorbar
cbar = plt.colorbar(im, pad=0.07)
cbar.set_label("Velocity (km/s)", size=16)

# Overlay set of RA/Dec Axes
overlay = ax.get_coords_overlay("fk5")
overlay.grid(color="white", ls="dotted", lw=2)
overlay[0].set_axislabel("Right Ascension (J2000)", fontsize=16)
overlay[1].set_axislabel("Declination (J2000)", fontsize=16)

# Overplot column density contours
levels = (1e20, 5e20, 1e21, 3e21, 5e21, 7e21, 1e22)  # Define contour levels to use
ax.contour(hi_column_density.hdu.data, cmap="Greys_r", alpha=0.5, levels=levels)
<matplotlib.contour.QuadContourSet at 0x7ff8a0cdfe60>
_images/7ea6258532f3b8288f0a99a7b00b5be61254d83625099953158e04a32eefee2e.png

As you can see, the WCSAxes framework is very powerful and similiar to making any matplotlib style plot.

Display a Longitude-Velocity Slice#

The WCSAxes framework in Astropy also lets us slice the data accross different dimensions. It is often useful to slice along a single latitude and display an image showing longtitude and velocity information only (position-velocity or longitude-velocity diagram).

This can be done by specifying the slices keyword and selecting the appropriate slice through the data.

slices requires a 3D tuple containing the index to be sliced along and where we want the two axes to be displayed. This should be specified in the same order as the WCS object (longitude, latitude, velocity) as opposed to the order of numpy array holding the data (velocity, latitude, longitude).

We then select the appropriate data by indexing along the numpy array.

lat_slice = 18  # Index of latitude dimension to slice along

# Initiate a figure and axis object with WCS projection information
fig = plt.figure(figsize=(18, 12))
ax = fig.add_subplot(111, projection=sub_cube_slab.wcs, slices=("y", lat_slice, "x"))
# Above, we have specified to plot the longitude along the y axis, pick only the lat_slice
# indicated, and plot the velocity along the x axis

# Display the slice
im = ax.imshow(
    sub_cube_slab[:, lat_slice, :].transpose().data
)  # Display the image slice
ax.invert_yaxis()  # Flips the Y axis

# Add axes labels
ax.set_xlabel("LSR Velocity (m/s)", fontsize=16)
ax.set_ylabel("Galactic Longitude (degrees)", fontsize=16)

# Add a colorbar
cbar = plt.colorbar(im, pad=0.07, orientation="horizontal")
cbar.set_label("Temperature (K)", size=16)
WARNING: WCSWarning: Slicing across a celestial axis results in an invalid WCS, so the celestial projection (TAN) is being removed.  The WCS indices being kept were [0 2]. [spectral_cube.wcs_utils]
WARNING: WCSWarning: Slicing across a celestial axis results in an invalid WCS, so the celestial projection (TAN) is being removed.  The view used was (slice(None, None, None), 18, slice(None, None, None)). [spectral_cube.wcs_utils]
_images/3cf78b5ca13f10f3aa2179da4f7aed610b78a12038de74e295efddaf64fa357c.png

As we can see, the SMC seems to be only along positive velocities.

Try:#

Create a new spectral slab isolating just the SMC and slice along a different dimension to create a latitude-velocity diagram#

Find and Download a Herschel Image#

This is great, but we want to compare the HI emission data with Herschel 350 micron emission to trace some dust. This can be done with astroquery. We can query for the data by mission, take a quick look at the table of results, and download data after selecting a specific wavelength or filter.

Since we are looking for Herschel data from an ESA mission, we will use the astroquery.ESASky class.

Specifically, the ESASKY.query_region_maps() method allows us to search for a specific region of the sky either using an Astropy SkyCoord object or a string specifying an object name. In this case, we can just search for the SMC. A radius to search around the object can also be specified.

# Query for Herschel data in a 1 degree radius around the SMC
result = ESASky.query_region_maps("SMC", radius=1 * u.deg, missions="Herschel")

print(result)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/coordinates/angles/formats.py:280, in _AngleParser.parse(self, angle, unit, debug)
    279 try:
--> 280     found_angle, found_unit = self._thread_local._parser.parse(
    281         angle, lexer=self._thread_local._lexer, debug=debug
    282     )
    283 except ValueError as e:

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/utils/parsing.py:114, in ThreadSafeParser.parse(self, *args, **kwargs)
    113 with self._lock:
--> 114     return self.parser.parse(*args, **kwargs)

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/extern/ply/yacc.py:333, in LRParser.parse(self, input, lexer, debug, tracking, tokenfunc)
    332 else:
--> 333     return self.parseopt_notrack(input, lexer, debug, tracking, tokenfunc)

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/extern/ply/yacc.py:1063, in LRParser.parseopt_notrack(self, input, lexer, debug, tracking, tokenfunc)
   1062 if not lookaheadstack:
-> 1063     lookahead = get_token()     # Get the next token
   1064 else:

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/extern/ply/lex.py:386, in Lexer.token(self)
    385 self.lexpos = lexpos
--> 386 newtok = self.lexerrorf(tok)
    387 if lexpos == self.lexpos:
    388     # Error method didn't change text position at all. This is an error.

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/coordinates/angles/formats.py:164, in _AngleParser._make_parser.<locals>.t_error(t)
    163 def t_error(t):
--> 164     raise ValueError(f"Invalid character at col {t.lexpos}")

ValueError: Invalid character at col 0

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/coordinates/sky_coordinate_parsers.py:532, in _parse_coordinate_arg(coords, frame, units)
    529     for frame_attr_name, repr_attr_class, value, unit in zip(
    530         frame_attr_names, repr_attr_classes, values, units
    531     ):
--> 532         components[frame_attr_name] = repr_attr_class(
    533             value, unit=unit, copy=COPY_IF_NEEDED
    534         )
    535 except Exception as err:

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/coordinates/angles/core.py:724, in Longitude.__new__(cls, angle, unit, wrap_angle, **kwargs)
    721     raise TypeError(
    722         "A Longitude angle cannot be created from a Latitude angle."
    723     )
--> 724 self = super().__new__(cls, angle, unit=unit, **kwargs)
    725 if wrap_angle is None:

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/coordinates/angles/core.py:163, in Angle.__new__(cls, angle, unit, dtype, copy, **kwargs)
    162 if isinstance(angle, str):
--> 163     angle, angle_unit = formats.parse_angle(angle, unit)
    164     if angle_unit is None:

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/coordinates/angles/formats.py:358, in parse_angle(angle, unit, debug)
    328 """
    329 Parses an input string value into an angle value.
    330 
   (...)    356     string.
    357 """
--> 358 return _AngleParser().parse(angle, unit, debug=debug)

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/coordinates/angles/formats.py:284, in _AngleParser.parse(self, angle, unit, debug)
    283 except ValueError as e:
--> 284     raise ValueError(
    285         f"{str(e) or 'syntax error'} parsing angle {angle!r}"
    286     ) from e
    288 if unit is None and found_unit is None:

ValueError: Invalid character at col 0 parsing angle 'SMC'

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astroquery/utils/commons.py:67, in parse_coordinates(coordinates)
     66 try:
---> 67     c = SkyCoord(coordinates, frame="icrs")
     68     warnings.warn("Coordinate string is being interpreted as an "
     69                   "ICRS coordinate.", InputWarning)

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/coordinates/sky_coordinate.py:231, in SkyCoord.__init__(self, copy, *args, **kwargs)
    230 args = list(args)  # Make it mutable
--> 231 skycoord_kwargs, components, info = _parse_coordinate_data(
    232     frame_cls(**frame_kwargs), args, kwargs
    233 )
    235 # In the above two parsing functions, these kwargs were identified
    236 # as valid frame attributes for *some* frame, but not the frame that
    237 # this SkyCoord will have. We keep these attributes as special
    238 # skycoord frame attributes:

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/coordinates/sky_coordinate_parsers.py:279, in _parse_coordinate_data(frame, args, kwargs)
    274 elif len(args) == 1:
    275     # One arg which must be a coordinate.  In this case coord_kwargs
    276     # will contain keys like 'ra', 'dec', 'distance' along with any
    277     # frame attributes like equinox or obstime which were explicitly
    278     # specified in the coordinate object (i.e. non-default).
--> 279     _skycoord_kwargs, _components = _parse_coordinate_arg(args[0], frame, units)
    281     # Copy other 'info' attr only if it has actually been defined.

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/coordinates/sky_coordinate_parsers.py:536, in _parse_coordinate_arg(coords, frame, units)
    535 except Exception as err:
--> 536     raise ValueError(
    537         f'Cannot parse first argument data "{value}" for attribute'
    538         f" {frame_attr_name}"
    539     ) from err
    540 return skycoord_kwargs, components

ValueError: Cannot parse first argument data "SMC" for attribute ra

During handling of the above exception, another exception occurred:

KeyboardInterrupt                         Traceback (most recent call last)
Cell In[14], line 2
      1 # Query for Herschel data in a 1 degree radius around the SMC
----> 2 result = ESASky.query_region_maps("SMC", radius=1 * u.deg, missions="Herschel")
      4 print(result)

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astroquery/esasky/core.py:703, in ESASkyClass.query_region_maps(self, position, radius, missions, get_query_payload, cache, row_limit, verbose)
    700 query_result = {}
    702 sesame_database.set('simbad')
--> 703 coordinates = commons.parse_coordinates(position)
    705 self._store_query_result(query_result=query_result, names=sanitized_missions,
    706                          descriptors=self._get_observation_info(), coordinates=coordinates,
    707                          radius=sanitized_radius, row_limit=sanitized_row_limit,
    708                          get_query_payload=get_query_payload, cache=cache, verbose=verbose)
    709 if get_query_payload:

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astroquery/utils/commons.py:86, in parse_coordinates(coordinates)
     84                 c = SkyCoord.from_name(coordinates, frame="icrs")
     85         else:
---> 86             c = SkyCoord.from_name(coordinates, frame="icrs")
     88 elif isinstance(coordinates, CoordClasses):
     89     if hasattr(coordinates, 'frame'):

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/coordinates/sky_coordinate.py:1909, in SkyCoord.from_name(cls, name, frame, parse, cache)
   1875 """
   1876 Given a name, query the CDS name resolver to attempt to retrieve
   1877 coordinate information for that object. The search database, sesame
   (...)   1905     Instance of the SkyCoord class.
   1906 """
   1907 from .name_resolve import get_icrs_coordinates
-> 1909 icrs_coord = get_icrs_coordinates(name, parse, cache=cache)
   1910 icrs_sky_coord = cls(icrs_coord)
   1911 if frame in ("icrs", icrs_coord.__class__):

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/coordinates/name_resolve.py:176, in get_icrs_coordinates(name, parse, cache)
    173 for url in urls:
    174     try:
    175         resp_data = get_file_contents(
--> 176             download_file(url, cache=cache, show_progress=False)
    177         )
    178         break
    179     except urllib.error.URLError as e:

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/utils/data.py:1548, in download_file(remote_url, cache, show_progress, timeout, sources, pkgname, http_headers, ssl_context, allow_insecure)
   1546 for source_url in sources:
   1547     try:
-> 1548         f_name = _download_file_from_source(
   1549             source_url,
   1550             timeout=timeout,
   1551             show_progress=show_progress,
   1552             cache=cache,
   1553             remote_url=remote_url,
   1554             pkgname=pkgname,
   1555             http_headers=http_headers,
   1556             ssl_context=ssl_context,
   1557             allow_insecure=allow_insecure,
   1558         )
   1559         # Success!
   1560         break

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/utils/data.py:1332, in _download_file_from_source(source_url, show_progress, timeout, remote_url, cache, pkgname, http_headers, ftp_tls, ssl_context, allow_insecure)
   1329         else:
   1330             raise
-> 1332 with _try_url_open(
   1333     source_url,
   1334     timeout=timeout,
   1335     http_headers=http_headers,
   1336     ftp_tls=ftp_tls,
   1337     ssl_context=ssl_context,
   1338     allow_insecure=allow_insecure,
   1339 ) as remote:
   1340     info = remote.info()
   1341     try:

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/site-packages/astropy/utils/data.py:1246, in _try_url_open(source_url, timeout, http_headers, ftp_tls, ssl_context, allow_insecure)
   1243 req = urllib.request.Request(source_url, headers=http_headers)
   1245 try:
-> 1246     return urlopener.open(req, timeout=timeout)
   1247 except urllib.error.URLError as exc:
   1248     reason = exc.reason

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/urllib/request.py:515, in OpenerDirector.open(self, fullurl, data, timeout)
    512     req = meth(req)
    514 sys.audit('urllib.Request', req.full_url, req.data, req.headers, req.get_method())
--> 515 response = self._open(req, data)
    517 # post-process response
    518 meth_name = protocol+"_response"

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/urllib/request.py:532, in OpenerDirector._open(self, req, data)
    529     return result
    531 protocol = req.type
--> 532 result = self._call_chain(self.handle_open, protocol, protocol +
    533                           '_open', req)
    534 if result:
    535     return result

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/urllib/request.py:492, in OpenerDirector._call_chain(self, chain, kind, meth_name, *args)
    490 for handler in handlers:
    491     func = getattr(handler, meth_name)
--> 492     result = func(*args)
    493     if result is not None:
    494         return result

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/urllib/request.py:1392, in HTTPSHandler.https_open(self, req)
   1391 def https_open(self, req):
-> 1392     return self.do_open(http.client.HTTPSConnection, req,
   1393                         context=self._context)

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/urllib/request.py:1344, in AbstractHTTPHandler.do_open(self, http_class, req, **http_conn_args)
   1342 try:
   1343     try:
-> 1344         h.request(req.get_method(), req.selector, req.data, headers,
   1345                   encode_chunked=req.has_header('Transfer-encoding'))
   1346     except OSError as err: # timeout error
   1347         raise URLError(err)

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/http/client.py:1338, in HTTPConnection.request(self, method, url, body, headers, encode_chunked)
   1335 def request(self, method, url, body=None, headers={}, *,
   1336             encode_chunked=False):
   1337     """Send a complete request to the server."""
-> 1338     self._send_request(method, url, body, headers, encode_chunked)

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/http/client.py:1384, in HTTPConnection._send_request(self, method, url, body, headers, encode_chunked)
   1380 if isinstance(body, str):
   1381     # RFC 2616 Section 3.7.1 says that text default has a
   1382     # default charset of iso-8859-1.
   1383     body = _encode(body, 'body')
-> 1384 self.endheaders(body, encode_chunked=encode_chunked)

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/http/client.py:1333, in HTTPConnection.endheaders(self, message_body, encode_chunked)
   1331 else:
   1332     raise CannotSendHeader()
-> 1333 self._send_output(message_body, encode_chunked=encode_chunked)

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/http/client.py:1093, in HTTPConnection._send_output(self, message_body, encode_chunked)
   1091 msg = b"\r\n".join(self._buffer)
   1092 del self._buffer[:]
-> 1093 self.send(msg)
   1095 if message_body is not None:
   1096 
   1097     # create a consistent interface to message_body
   1098     if hasattr(message_body, 'read'):
   1099         # Let file-like take precedence over byte-like.  This
   1100         # is needed to allow the current position of mmap'ed
   1101         # files to be taken into account.

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/http/client.py:1037, in HTTPConnection.send(self, data)
   1035 if self.sock is None:
   1036     if self.auto_open:
-> 1037         self.connect()
   1038     else:
   1039         raise NotConnected()

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/http/client.py:1472, in HTTPSConnection.connect(self)
   1469 def connect(self):
   1470     "Connect to a host on a given (SSL) port."
-> 1472     super().connect()
   1474     if self._tunnel_host:
   1475         server_hostname = self._tunnel_host

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/http/client.py:1003, in HTTPConnection.connect(self)
   1001 """Connect to the host and port specified in __init__."""
   1002 sys.audit("http.client.connect", self, self.host, self.port)
-> 1003 self.sock = self._create_connection(
   1004     (self.host,self.port), self.timeout, self.source_address)
   1005 # Might fail in OSs that don't implement TCP_NODELAY
   1006 try:

File /opt/hostedtoolcache/Python/3.12.11/x64/lib/python3.12/socket.py:850, in create_connection(address, timeout, source_address, all_errors)
    848 if source_address:
    849     sock.bind(source_address)
--> 850 sock.connect(sa)
    851 # Break explicitly a reference cycle
    852 exceptions.clear()

KeyboardInterrupt: 

Here, the result is a TableList which contains 24 Herschel data products that can be downloaded. We can see what information is available in this TableList by examining the keys in the Herschel Table.

result["HERSCHEL"].keys()

We want to find a 350 micron image, so we need to look closer at the filters used for these observations.

result["HERSCHEL"]["filter"]

Luckily for us, there is an observation made with three filters: 250, 350, and 500 microns. This is the object we will want to download. One way to do this is by making a boolean mask to select out the Table entry corresponding with the desired filter. Then, the ESASky.get_maps() method will download our data provided a TableList argument.

Note that the below command requires an internet connection to download five quite large files. It could take several minutes to complete.
filters = result["HERSCHEL"]["filter"].astype(
    str
)  # Convert the list of filters from the query to a string

# Construct a boolean mask, searching for only the desired filters
mask = np.array(["250, 350, 500" == s for s in filters], dtype="bool")

# Re-construct a new TableList object containing only our desired query entry
target_obs = TableList(
    {"HERSCHEL": result["HERSCHEL"][mask]}
)  # This will be passed into ESASky.get_maps()

IR_images = ESASky.get_maps(target_obs)  # Download the images
IR_images["HERSCHEL"][0][
    "350"
].info()  # Display some information about the 350 micron image

Since we are just doing some qualitative analysis, we only need the image, but you can also access lots of other information from our downloaded object, such as errors.

Let’s go ahead and extract just the WCS information and image data from the 350 micron image.

herschel_header = IR_images["HERSCHEL"][0]["350"]["image"].header
herschel_wcs = WCS(IR_images["HERSCHEL"][0]["350"]["image"])  # Extract WCS information
herschel_imagehdu = IR_images["HERSCHEL"][0]["350"]["image"]  # Extract Image data
print(herschel_wcs)

With this, we can display this image using matplotlib with WCSAxes and the LogNorm() object so we can log scale our image.

# Set Nans to zero
himage_nan_locs = np.isnan(herschel_imagehdu.data)
herschel_data_nonans = herschel_imagehdu.data
herschel_data_nonans[himage_nan_locs] = 0

# Initiate a figure and axis object with WCS projection information
fig = plt.figure(figsize=(18, 12))
ax = fig.add_subplot(111, projection=herschel_wcs)

# Display the moment map image
im = ax.imshow(herschel_data_nonans, cmap="viridis", norm=LogNorm(vmin=2, vmax=50))
# ax.invert_yaxis() # Flips the Y axis

# Add axes labels
ax.set_xlabel("Right Ascension", fontsize=16)
ax.set_ylabel("Declination", fontsize=16)
ax.grid(color="white", ls="dotted", lw=2)

# Add a colorbar
cbar = plt.colorbar(im, pad=0.07)
cbar.set_label(
    "".join(["Herschel 350" r"$\mu$m ", "(", herschel_header["BUNIT"], ")"]), size=16
)

# Overlay set of Galactic Coordinate Axes
overlay = ax.get_coords_overlay("galactic")
overlay.grid(color="black", ls="dotted", lw=1)
overlay[0].set_axislabel("Galactic Longitude", fontsize=14)
overlay[1].set_axislabel("Galactic Latitude", fontsize=14)

Overlay HI 21 cm Contours on the IR 30 micron Image#

To visually compare the neutral gas and dust as traced by HI 21 cm emission and IR 30 micron emission, we can use contours and colorscale images produced using the WCSAxes framework and the .get_transform() method.

The WCSAxes.get_transform() method returns a transformation from a specified frame to the pixel/data coordinates. It accepts a string specifying the frame or a WCS object.

# Initiate a figure and axis object with WCS projection information
fig = plt.figure(figsize=(18, 12))
ax = fig.add_subplot(111, projection=herschel_wcs)

# Display the moment map image
im = ax.imshow(
    herschel_data_nonans, cmap="viridis", norm=LogNorm(vmin=5, vmax=50), alpha=0.8
)
# ax.invert_yaxis() # Flips the Y axis

# Add axes labels
ax.set_xlabel("Right Ascension", fontsize=16)
ax.set_ylabel("Declination", fontsize=16)
ax.grid(color="white", ls="dotted", lw=2)

# Extract x and y coordinate limits
x_lim = ax.get_xlim()
y_lim = ax.get_ylim()

# Add a colorbar
cbar = plt.colorbar(im, fraction=0.046, pad=-0.1)
cbar.set_label(
    "".join(["Herschel 350" r"$\mu$m ", "(", herschel_header["BUNIT"], ")"]), size=16
)

# Overlay set of RA/Dec Axes
overlay = ax.get_coords_overlay("galactic")
overlay.grid(color="black", ls="dotted", lw=1)
overlay[0].set_axislabel("Galactic Longitude", fontsize=14)
overlay[1].set_axislabel("Galactic Latitude", fontsize=14)

hi_transform = ax.get_transform(
    hi_column_density.wcs
)  # extract axes Transform information for the HI data

# Overplot column density contours
levels = (2e21, 3e21, 5e21, 7e21, 8e21, 1e22)  # Define contour levels to use
ax.contour(
    hi_column_density.hdu.data,
    cmap="Greys_r",
    alpha=0.8,
    levels=levels,
    transform=hi_transform,
)  # include the transform information with the keyword "transform"

# Overplot velocity image so we can also see the Gas velocities
im_hi = ax.imshow(
    moment_1.hdu.data,
    cmap="RdBu_r",
    vmin=0,
    vmax=200,
    alpha=0.5,
    transform=hi_transform,
)

# Add a second colorbar for the HI Velocity information
cbar_hi = plt.colorbar(im_hi, orientation="horizontal", fraction=0.046, pad=0.07)
cbar_hi.set_label("HI " r"$21$cm Mean Velocity (km/s)", size=16)

# Apply original image x and y coordinate limits
ax.set_xlim(x_lim)
ax.set_ylim(y_lim)

Using reproject to match image resolutions#

The reproject package is a powerful tool allowing for image data to be transformed into a variety of projections and resolutions. It’s most powerful use is in fact to transform data from one map projection to another without losing any information and still properly conserving flux values within the data. It even has a method to perform a fast reprojection if you are not too concerned with the absolute accuracy of the data values.

A simple use of the reproject package is to scale down (or up) resolutions of an image artificially. This could be a useful step if you are trying to get emission line ratios or directly compare the Intensity or Flux from a tracer to that of another tracer in the same physical point of the sky.

From our previously made images, we can see that the IR Herschel Image has a higher spatial resolution than that of the HI data cube. We can look into this more by taking a better look at both header objects and using reproject to downscale the Herschel Image.

print("IR Resolution (dx,dy) = ", herschel_header["cdelt1"], herschel_header["cdelt2"])

print(
    "HI Resolution (dx,dy) = ",
    hi_column_density.hdu.header["cdelt1"],
    hi_column_density.hdu.header["cdelt1"],
)
Note: Different ways of accessing the header are shown above corresponding to the different object types (coming from SpectralCube vs astropy.io.fits)

As we can see, the IR data has over 10 times higher spatial resolution. In order to create a new projection of an image, all we need to specifiy is a new header containing WCS information to transform into. These can be created manually if you wanted to completely change something about the projection type (i.e. going from a Mercator map projection to a Tangential map projection). For us, since we want to match our resolutions, we can just “steal” the WCS object from the HI data. Specifically, we will be using the reproject_interp() function. This takes two arguments: an HDU object that you want to reproject, and a header containing WCS information to reproject onto.

rescaled_herschel_data, _ = reproject_interp(
    herschel_imagehdu,
    # reproject the Herschal image to match the HI data
    hi_column_density.hdu.header,
)

rescaled_herschel_imagehdu = fits.PrimaryHDU(
    data=rescaled_herschel_data,
    # wrap up our reprojection as a new fits HDU object
    header=hi_column_density.hdu.header,
)

rescaled_herschel_imagehdu will now behave just like the other FITS images we have been working with, but now with a degraded resolution matching the HI data. This includes having its native coordinates in Galactic rather than RA and Dec.

# Set Nans to zero
image_nan_locs = np.isnan(rescaled_herschel_imagehdu.data)
rescaled_herschel_data_nonans = rescaled_herschel_imagehdu.data
rescaled_herschel_data_nonans[image_nan_locs] = 0

# Initiate a figure and axis object with WCS projection information
fig = plt.figure(figsize=(18, 12))
ax = fig.add_subplot(111, projection=WCS(rescaled_herschel_imagehdu))

# Display the moment map image
im = ax.imshow(
    rescaled_herschel_data_nonans,
    cmap="viridis",
    norm=LogNorm(vmin=5, vmax=50),
    alpha=0.8,
)
# im = ax.imshow(rescaled_herschel_imagehdu.data, cmap = 'viridis',
#               norm = LogNorm(), vmin = 5, vmax = 50, alpha = .8)
ax.invert_yaxis()  # Flips the Y axis

# Add axes labels
ax.set_xlabel("Galactic Longitude", fontsize=16)
ax.set_ylabel("Galactic Latitude", fontsize=16)
ax.grid(color="white", ls="dotted", lw=2)

# Extract x and y coordinate limits
x_lim = ax.get_xlim()
y_lim = ax.get_ylim()

# Add a colorbar
cbar = plt.colorbar(im, fraction=0.046, pad=-0.1)
cbar.set_label(
    "".join(["Herschel 350" r"$\mu$m ", "(", herschel_header["BUNIT"], ")"]), size=16
)

# Overlay set of RA/Dec Axes
overlay = ax.get_coords_overlay("fk5")
overlay.grid(color="black", ls="dotted", lw=1)
overlay[0].set_axislabel("Right Ascension", fontsize=14)
overlay[1].set_axislabel("Declination", fontsize=14)

hi_transform = ax.get_transform(
    hi_column_density.wcs
)  # extract axes Transform information for the HI data

# Overplot column density contours
levels = (2e21, 3e21, 5e21, 7e21, 8e21, 1e22)  # Define contour levels to use
ax.contour(
    hi_column_density.hdu.data,
    cmap="Greys_r",
    alpha=0.8,
    levels=levels,
    transform=hi_transform,
)  # include the transform information with the keyword "transform"

# Overplot velocity image so we can also see the Gas velocities
im_hi = ax.imshow(
    moment_1.hdu.data,
    cmap="RdBu_r",
    vmin=0,
    vmax=200,
    alpha=0.5,
    transform=hi_transform,
)

# Add a second colorbar for the HI Velocity information
cbar_hi = plt.colorbar(im_hi, orientation="horizontal", fraction=0.046, pad=0.07)
cbar_hi.set_label("HI " r"$21$cm Mean Velocity (km/s)", size=16)

# Apply original image x and y coordinate limits
ax.set_xlim(x_lim)
ax.set_ylim(y_lim)

The real power of reproject is in actually changing the map projection used to display the data. This is done by creating a WCS object that contains a different projection type such as CTYPE : 'RA---CAR'  'DEC--CAR' as opposed to CTYPE : 'RA---TAN'  'DEC--TAN'.

Challenge:#

Use reproject and WCS to create a new WCS object in a different map projection and see how distortions in the image can change.