Surveying dust structure via GALEX MIS

Surveying dust structure via GALEX MIS#

The objective of this tutorial is to create a UV mosaic of a high-latitute cloud, using data from GALEX. We’ll query the MAST Archive for Observations, then use the data to create a composite image.

Learning goals#

By the end of this tutorial you will be able to:

  • Understand UV images and how they are useful to study dust

  • Access MAST cloud data

  • Create masks for circular images

  • Create a mosaic from several images using WCS

Introduction#

GALEX background: The Galaxy Evolution Explorer (GALEX) was a satellite designed to investigate star formation. It observed the sky in two different bands: in the Near UV (NUV) (\(1750-27504\) Å) and in the Far UV (FUV) (\(1350-1750\) Å). The GALEX database contains over 600 million source measurements in the ultraviolet domain, with some sources having more than one measurement, making them useful for studying their variable sources.

GALEX Medium Imaging survey (MIS) background: Single orbit exposures (1,500s) of 1000 square degrees in positions that match the Sloan Digital Sky Survey (SDSS). These images present a higher resolution in comparison with AIS, since their exposure time was longer.

Clouds can be visible in UV when they are found close to hot stars. The objective of this tutorial is to extract and display an intensity map image of a high-latitude cloud retrieved from GALEX MIS. High-latitude clouds (Galactic latitude: \(|b| > 20-30º\)) are interesting because they are considered ideal candidates to study triggered star formation.

Table of Contents#

  • Imports

  • MBM 15 Search in GALEX MIS

  • Product Selection and Download

  • Plotting a Test Image

  • Data Downsizing

  • Creating a Circular Mask

  • Assembling the Mosaic

  • Plotting the Mosaic

Imports#

  • astropy to use tools needed for performing astronomy and astrophysics with Python, including handling fits files, defining coordinates or choosing the right limits for image visualization.

  • Observations from astroquery.mast to query the Barbara A. Mikulski Archive for Space Telescopes (MAST).

  • matplotlib to visualize images.

  • numpy to manipulate arrays.

  • reproject to select the correct WCS for our images

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np

from astropy.io import fits
from astroquery.mast import Observations
from astropy.visualization import ZScaleInterval
from photutils.aperture import CircularAperture
from reproject import reproject_interp
from reproject.mosaicking import find_optimal_celestial_wcs

# turn on access to the cloud dataset
Observations.enable_cloud_dataset()

INFO: Using the S3 STScI public dataset [astroquery.mast.cloud]

MBM 15 Search in GALEX MIS#

We’ll access GALEX MIS products using astroquery.mast. Let’s search for any Observations within a degree of our target, the molecular cloud MBM 15. We should also specify mission with the filter obs_collection = GALEX.

Note that we add two additional filters to our query below. First, we use project="MIS" to specify the medium imaging survey. There are a few imaging surveys we could choose from, but MIS will give us good coverage of the cloud with a longer exposure time than the all-sky survey. GALEX has two detectors, so let’s also choose the Near UV with filters="NUV". Since we’re just looking to create a mosaic, it doesn’t matter too much which filter we pick; you might have a valid scientific reason to prefer one over the other.

obs = Observations.query_criteria(objectname="MBM 15", radius="1 deg", obs_collection="GALEX", project="MIS", filters="NUV")
len(obs)

9

Product Selection and Download#

With nine Observations, we should be able to a mosaic with a good view of the nebula. Now we need to get the data products: the actual files associated with these Observations.

prod = Observations.get_product_list(obs)
len(prod)

2592

Woah! Over 2000 files from just 9 Observations.

However, many of these files are not useful for our analysis. We’ll filter for files that are:

  • Marked as "SCIENCE" type: that is, not calibration or data validation files

  • Imaging only data: no catalog files.

  • We’ll also keep only files that have been flagged as a “minimum recommended products”. These files are selected by MAST archive scientists as the most relevant for scientific analysis. They should work well for our mosaic.

# filter files on the above specifications
filt_prod = Observations.filter_products(prod, productType="SCIENCE", productSubGroupDescription="Imaging Only", mrp_only=True)

# get the filenames
cloud_uris = Observations.get_cloud_uris(filt_prod)
cloud_uris

['s3://stpubdata/galex/GR6/pipe/01-vsn/03329-MISDR1_18916_0459/d/01-main/0001-img/07-try/MISDR1_18916_0459-fd-int.fits.gz',
 's3://stpubdata/galex/GR6/pipe/01-vsn/03329-MISDR1_18916_0459/d/01-main/0001-img/07-try/MISDR1_18916_0459-nd-int.fits.gz',
 's3://stpubdata/galex/GR6/pipe/01-vsn/17464-MISWZS03_18917_0284/d/01-main/0001-img/07-try/MISWZS03_18917_0284-fd-int.fits.gz',
 's3://stpubdata/galex/GR6/pipe/01-vsn/17464-MISWZS03_18917_0284/d/01-main/0001-img/07-try/MISWZS03_18917_0284-nd-int.fits.gz',
 's3://stpubdata/galex/GR6/pipe/01-vsn/17477-MISWZS03_27307_0183/d/01-main/0001-img/07-try/MISWZS03_27307_0183-fd-int.fits.gz',
 's3://stpubdata/galex/GR6/pipe/01-vsn/17477-MISWZS03_27307_0183/d/01-main/0001-img/07-try/MISWZS03_27307_0183-nd-int.fits.gz',
 's3://stpubdata/galex/GR6/pipe/01-vsn/17542-MISWZS03_28512_0284/d/01-main/0001-img/07-try/MISWZS03_28512_0284-fd-int.fits.gz',
 's3://stpubdata/galex/GR6/pipe/01-vsn/17542-MISWZS03_28512_0284/d/01-main/0001-img/07-try/MISWZS03_28512_0284-nd-int.fits.gz',
 's3://stpubdata/galex/GR6/pipe/01-vsn/17543-MISWZS03_28513_0284/d/01-main/0001-img/07-try/MISWZS03_28513_0284-fd-int.fits.gz',
 's3://stpubdata/galex/GR6/pipe/01-vsn/17543-MISWZS03_28513_0284/d/01-main/0001-img/07-try/MISWZS03_28513_0284-nd-int.fits.gz',
 's3://stpubdata/galex/GR6/pipe/01-vsn/17573-MISWZS03_27260_0183o/d/01-main/0001-img/07-try/MISWZS03_27260_0183o-fd-int.fits.gz',
 's3://stpubdata/galex/GR6/pipe/01-vsn/17573-MISWZS03_27260_0183o/d/01-main/0001-img/07-try/MISWZS03_27260_0183o-nd-int.fits.gz',
 's3://stpubdata/galex/GR6/pipe/01-vsn/17574-MISWZS03_28514_0459o/d/01-main/0001-img/07-try/MISWZS03_28514_0459o-fd-int.fits.gz',
 's3://stpubdata/galex/GR6/pipe/01-vsn/17574-MISWZS03_28514_0459o/d/01-main/0001-img/07-try/MISWZS03_28514_0459o-nd-int.fits.gz',
 's3://stpubdata/galex/GR7/pipe/02-vsn/53017-MISWZS03_18917_0284_css29469/d/01-main/0001-img/07-try/MISWZS03_18917_0284_css29469-fd-int.fits.gz',
 's3://stpubdata/galex/GR7/pipe/02-vsn/53017-MISWZS03_18917_0284_css29469/d/01-main/0001-img/07-try/MISWZS03_18917_0284_css29469-nd-int.fits.gz',
 's3://stpubdata/galex/GR7/pipe/02-vsn/53022-MISWZS03_27307_0183_css24257/d/01-main/0001-img/07-try/MISWZS03_27307_0183_css24257-fd-int.fits.gz',
 's3://stpubdata/galex/GR7/pipe/02-vsn/53022-MISWZS03_27307_0183_css24257/d/01-main/0001-img/07-try/MISWZS03_27307_0183_css24257-nd-int.fits.gz']

Great! Now, you might notice that the files end with either fd-int.fits.gz or nd-int.fits.gz; due to the way GALEX Observations are handled, we’re actually getting both the fd (far UV) and nd (near UV) data. We’ll handle this gracefully a few cells from now.

In the meantime, let’s open the files into a list of fits objects:

# need "anon":true to anonymously access the cloud files
hdus = [fits.open(filename, fsspec_kwargs={"anon": True}) for filename in cloud_uris]

Since that cell doesn’t have any output, we should run a test to make sure everything is working.

Plotting a Test Image#

Let’s run our test by getting the data from the first file in the list and plotting it. We won’t worry (yet) about plotting the correct coordinates, although we will need to appropriately scale the data based on brightness so it’s visible in our plot.

# Create a figure on which to plot our data
fig = plt.figure(0, [9, 9])
ax = fig.add_subplot(111)

# Get the primary data from the first fits file
test_data = hdus[0][0].data

# Automatically scale the brightness based on the data
interval = ZScaleInterval(contrast=0.2)
lims = interval.get_limits(test_data)

# Show our data with the scaling from above
ax.imshow(test_data, vmin=lims[0], vmax=lims[1], cmap='BuPu_r')

<matplotlib.image.AxesImage at 0x7f008c80dc90>
../../../_images/8419deae84c55d8931a4fedb6046ad6c4a19ea837120c5ca660d5ec2585f9d89.png

Excellent. We can aleady see some of the “wispy” or cloud-like structure of our target.

We should take note of the shape of the data: namely, it is circular. We need to keep this in mind as we construct our mosaic.

Data Downsizing#

These .fits files that we’ve downloaded are actually quite large. To keep Notebook runtime short, we’re going to scale down the resolution of the image. The “downgrade factor” is, in effect, making the pixels larger and thus have lower memory/processing requirements. With a default downgrade factor of five, our most complicated cell (coming later) will take about 30 seconds to run. To see the original quality image, you can set downgrade_factor = 1 (i.e. original quality).

At this step, we also need to find an optimal WCS that is big enough to display our combined mosaic. In essence, we know how big each image is and where it was taken on the sky; find_optimal_celestial_wcs will do the hard work for us and select an appropriate image size and coordinates.

# determine the angular resolution of the image
res = np.abs(hdus[0][0].header['CDELT1'])

# set up a "downgrade factor" to decrease the resolution of the final combined mosaic, because these files are BIG
downgrade_factor = 5     # = 1 for original quality

# use astropy's find_optimal_celestial_wcs to determine the output WCS 
# for the combined image based on the information in the headers
wcs_out, shape_out = find_optimal_celestial_wcs(hdus, resolution=res*downgrade_factor*u.deg)

Creating a Circular Mask#

Since the GALEX data is circular, we need to create a mask for the data. What we want is a set of valid values within the \(0.6^\circ\) observation, surrounded by NaNs.

Our function below takes map_in, keeping only valid values within the circular aperture. All other values are set to NaN; this allows us to compute a meaningful “average” value where the observations overlap.

def galex_mask(map_in, header):
    # define the radius of the circle to be 0.6 degrees: use header information to convert this to pixels
    dr = np.rint(0.6 / header['CDELT2'])
    app = CircularAperture((header['CRPIX1'], header['CRPIX2']), r=dr)
    
    # convert aperture to mask
    amask = app.to_mask(method='center')
    # pad the mask
    amask = np.pad(np.array(amask.data), 480)
    # convert mask to boolean values
    amask = np.array(amask, dtype=bool)
    
    # mask all pixels outside of aperture to "NaN"
    map_in[~amask] = np.nan
    
    return map_in

We’ll use this function in our “assembly” below.

Assembling the Mosaic#

It’s time to put all of the pieces together! Let’s combine these images, following the steps below:

  1. Get the data from each HDU.

  2. Only keep the NUV data.

  3. Set all data outside the aperture to NaN.

  4. Reproject the data onto the mosaic WCS to correctly orient all of the images.

  5. Take the average of the all the maps to gracefully handle any overlapping regions.

Note: this cell will take roughly 30 seconds to run with a downscaling factor of five. If you changed downgrade_factor to something lower than five, it will likely take longer.

# Create an empty "maps" list to hold our outputs
maps = []

# Loop through the hdus, getting the data and header from each
for hdu in hdus:
    fmap = hdu[0].data
    fhead = hdu[0].header
    
    # Only keep the NUV bands. (FUV = 1)
    if fhead['band'] == 2:
        # for each image, use our function to set all pixels outside of the GALEX footprint to NaN
        fmap = galex_mask(fmap, fhead)

        # use reproject to reproject each image to the "combined" output WCS defined earlier
        regmap, foot = reproject_interp((fmap, fhead), wcs_out, shape_out)

        # append the reprojected map to a list
        maps.append(regmap)
        
# average all of the reprojected maps together
combined = np.nanmean(np.dstack((maps)), 2)

/tmp/ipykernel_2056/1879320563.py:21: RuntimeWarning: Mean of empty slice
  combined = np.nanmean(np.dstack((maps)), 2)

Excellent. Now it’s time for the moment of truth.

Plotting the Mosaic#

Just like with our test image, we’ll need to create a figure. In this instance, we want to specify the projection to be our mosaic WCS so that the coordinates are labeled.

Otherwise, we proceed as normal.

# Create the figure for the new map
fig = plt.figure(0, [9, 9])
ax = fig.add_subplot(111, projection=wcs_out)

# Automatically scale the brightness
interval = ZScaleInterval(contrast=0.4)
lims = interval.get_limits(combined)

# Plot the new map
ax.imshow(combined, vmin=lims[0], vmax=lims[1], origin='lower', cmap='BuPu_r')
ax.set_xlabel('RA')
ax.set_ylabel('Dec')

# Un-comment the below to save the figure as a jpg
#plt.savefig('mbm15.jpg', dpi=600, bbox_inches='tight')
../../../_images/579196ad729325361d99a2373fbe026a97e3ca976523be4c6442d89013568c09.png

Beautiful! Our mosaic is complete. You may notice some artifacts where the data was incomplete; you can also see an arc in the lower left corner where our circular NaN aperture didn’t quite match the data. How might you improve upon this?

This image corresponds to the Molecular cloud MBM 15. This molecular cloud belongs to the Orion-Eridanus superbubble, west of the Orion Nebula. More information about this cloud can be found in Joubaud et al. (2019).

Additional Resources#

For more information about the MAST archive and details about the tutorial:

About this Notebook#

Authors: Clara Puerto Sánchez, Claire Murray, Thomas Dutkiewicz
Keywords: mosaic, dust, GALEX
Last Updated: October 2024

For support, please contact the Archive HelpDesk at archive@stsci.edu.

Top of Page Space Telescope Logo

Contents