Accessing the Map

Before diving into how you can access the 3D dust map, we'd like to make you aware of a few important points about the dataset. For a detailed discussion of the map, please read Green et al. (2015).

We use the same definition of E(B-V ) as the Schlegel, Finkbeiner & Davis (1998) dust map. Although this was originally supposed to be the excess B-V in the Landolt filter system, Schlafly & Finkbeiner (2011) found that it differs somewhat from the true stellar B-V excess. Therefore, in order to convert our values of E(B-V ) to extinction in various filter systems, consult Table 6 of Schlafly & Finkbeiner (2011) (use the values in the RV = 3.1 column).

For each sightline, we provide multiple estimates of the distance vs. reddening relationship. Alongside the maximum-probability density estimate (essentially, the best-fit) distance-reddening curve, we also provide samples of the distance-reddening curve, which are representative of the formal uncertainty in the relation. Most statistics you may wish to derive, like the median reddening to a given distance, are best determined by using the representative samples, rather than the best-fit relation.

We include a number of pieces of information on the reliability of each pixel. A convergence flag marks whether our fit to the line-of-sight reddening curve converged. This is a formal indicator, meaning that we correctly sampled the spread of possible distance-reddening relations, given our model assumptions. It does not, however, indicate that our model assumptions were correct for that pixel. This convergence flag is based on the Gelman-Rubin diagnostic, a method of flagging Markov Chain Monte Carlo non-convergence.

Additionally, minimum and maximum reliable distances are provided for each pixel, based on the distribution of stars along the sightline. Because we determine the presence of dust by observing differential reddening of foreground and background stars, we cannot trace dust beyond the farthest stars in a given pixel. Our estimates of dust reddening closer than the nearest observed stars in a pixel are similarly uncertain. We therefore urge caution in using our maps outside of the distance range indicated for each pixel.

You can query the Argonaut web server directly from your code, without having to download the full map to your local disk. If you're not interested in playing around with the entire map, this is the way we recommend you access it.

In , the following function will call the Argonaut server and return the line-of-sight reddening for any coordinates, or set of coordinates, on the sky:
import json, requests

def query(lon, lat, coordsys='gal', mode='full'):
    '''
    Send a line-of-sight reddening query to the Argonaut web server.
    
    Inputs:
      lon, lat: longitude and latitude, in degrees.
      coordsys: 'gal' for Galactic, 'equ' for Equatorial (J2000).
      mode: 'full', 'lite' or 'sfd'
    
    In 'full' mode, outputs a dictionary containing, among other things:
      'distmod':    The distance moduli that define the distance bins.
      'best':       The best-fit (maximum proability density)
                    line-of-sight reddening, in units of SFD-equivalent
                    E(B-V), to each distance modulus in 'distmod.' See
                    Schlafly & Finkbeiner (2011) for a definition of the
                    reddening vector (use R_V = 3.1).
      'samples':    Samples of the line-of-sight reddening, drawn from
                    the probability density on reddening profiles.
      'success':    1 if the query succeeded, and 0 otherwise.
      'converged':  1 if the line-of-sight reddening fit converged, and
                    0 otherwise.
      'n_stars':    # of stars used to fit the line-of-sight reddening.
      'DM_reliable_min':  Minimum reliable distance modulus in pixel.
      'DM_reliable_max':  Maximum reliable distance modulus in pixel.
    
    Less information is returned in 'lite' mode, while in 'sfd' mode,
    the Schlegel, Finkbeiner & Davis (1998) E(B-V) is returned.
    '''
    
    url = 'http://argonaut.skymaps.info/gal-lb-query-light'
    
    payload = {'mode': mode}
    
    if coordsys.lower() in ['gal', 'g']:
        payload['l'] = lon
        payload['b'] = lat
    elif coordsys.lower() in ['equ', 'e']:
        payload['ra'] = lon
        payload['dec'] = lat
    else:
        raise ValueError("coordsys '{0}' not understood.".format(coordsys))
    
    headers = {'content-type': 'application/json'}
    
    r = requests.post(url, data=json.dumps(payload), headers=headers)
    
    try:
        r.raise_for_status()
    except requests.exceptions.HTTPError as e:
        print('Response received from Argonaut:')
        print(r.text)
        raise e
    
    return json.loads(r.text)

This code can be adapted for any programming language that can issue HTTP POST requests. The code can also be found on GitHub.

To query one sightline, say, Galactic coordinates (ℓ, b) = (90°, 10°), you would call

>>> # Query the Galactic coordinates (l, b) = (90, 10):
>>> qresult = query(90, 10, coordsys='gal')
>>> 
>>> # See what information is returned for each pixel:
>>> qresult.keys()
[u'b', u'GR', u'distmod', u'l', u'DM_reliable_max', u'ra', u'samples',
u'n_stars', u'converged', u'success', u'dec', u'DM_reliable_min', u'best']
>>> 
>>> qresult['n_stars']
750
>>> qresult['converged']
1
>>> # Get the best-fit E(B-V) in each distance slice
>>> qresult['best']
[0.00426, 0.00678, 0.0074, 0.00948, 0.01202, 0.01623, 0.01815, 0.0245,
0.0887, 0.09576, 0.10139, 0.12954, 0.1328, 0.21297, 0.23867, 0.24461,
0.37452, 0.37671, 0.37684, 0.37693, 0.37695, 0.37695, 0.37696, 0.37698,
0.37698, 0.37699, 0.37699, 0.377, 0.37705, 0.37708, 0.37711]
>>> 
>>> # See the distance modulus of each distance slice
>>> qresult['distmod']
[4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0, 10.5,
11.0, 11.5, 12.0, 12.5, 13.0, 13.5, 14.0, 14.5, 15.0, 15.5, 16.0, 16.5,
17.0, 17.5, 18.0, 18.5, 19.0]

You can also query multiple sightlines simultaneously, simply by passing longitude and latitude as lists:

>>> qresult = query([45, 170, 250], [0, -20, 40])
>>> 
>>> qresult['n_stars']
[352, 162, 254]
>>> qresult['converged']
[1, 1, 1]
>>> # Look at the best fit for the first pixel:
>>> qresult['best'][0]
[0.00545, 0.00742, 0.00805, 0.01069, 0.02103, 0.02718, 0.02955, 0.03305,
0.36131, 0.37278, 0.38425, 0.41758, 1.53727, 1.55566, 1.65976, 1.67286,
1.78662, 1.79262, 1.88519, 1.94605, 1.95938, 2.0443, 2.39438, 2.43858,
2.49927, 2.54787, 2.58704, 2.58738, 2.58754, 2.58754, 2.58755]

If you are going to be querying large numbers of sightlines at once, we kindly request that you use this batch syntax, rather than calling query() in a loop. It will be faster, because you only have to contact the server once, and it will reduce the load on the Argonaut server.

Two additional query modes are provided, beyond the default 'full' mode. If mode = 'lite' is passed to the query() function, then less information is returned per sightline:

>>> qresult = query(180, 0, coordsys='gal', mode='lite')
>>> qresult.keys()
[u'b', u'success', u'distmod', u'sigma', u'median', u'l',
u'DM_reliable_max', u'ra', u'n_stars', u'converged', u'dec',
u'DM_reliable_min', u'best']
>>> 
>>> # Get the median E(B-V) to each distance slice:
>>> qresult['median']
[0.0204, 0.02747, 0.03027, 0.03036, 0.03047, 0.05214, 0.05523, 0.0748,
0.07807, 0.10002, 0.13699, 0.2013, 0.20158, 0.20734, 0.23129, 0.73734,
0.76125, 0.83905, 0.90236, 1.05944, 1.08085, 1.11408, 1.11925, 1.12212,
1.12285, 1.12289, 1.12297, 1.12306, 1.12308, 1.12309, 1.12312]
>>> 
>>> # Get the standard deviation of E(B-V) in each slice
>>> # (actually, half the difference between the 84th and 16th percentiles):
>>> qresult['sigma']
[0.03226, 0.03476, 0.03452, 0.03442, 0.03439, 0.03567, 0.03625, 0.0317,
0.03238, 0.03326, 0.05249, 0.0401, 0.03919, 0.03278, 0.08339, 0.05099,
0.03615, 0.04552, 0.05177, 0.03678, 0.03552, 0.05246, 0.05055, 0.05361,
0.05422, 0.0538, 0.05381, 0.05381, 0.0538, 0.0538, 0.05379]

Finally, for the convenience of many users who also want to query the two-dimensional Schlegel, Finkbeiner & Davis (1998) map of dust reddening, the option mode = 'sfd' is also provided:

>>> qresult = query([0, 10, 15], [75, 80, 85], coordsys='gal', mode='sfd')
>>> 
>>> qresult.keys()
[u'EBV_SFD', u'b', u'dec', u'l', u'ra']
>>> 
>>> # E(B-V), in magnitudes:
>>> qresult['EBV_SFD']
[0.02119, 0.01813, 0.01352]

The entire 3D map can be downloaded in either HDF5 or FITS format from the Harvard Dataverse.

The full map comes to over 4.5 GB in compressed HDF5 format, so if you're only interested in individual sightlines, we strongly recommend you use our web API, which you can include directly as a function call in a number of programming languages.

The HDF5 format is a self-documenting, highly flexible format for scientific data. It has a number of powerful features, such as internal compression and compound datatypes (similar to numpy structured arrays), and has bindings in many different programming languages, including C, Python, Fortran and IDL.

The HDF5 file we provide has four datasets:

  • /pixel_info : pixel locations and metadata.
  • /samples : samples of distance vs. reddening profile in each pixel.
  • /best_fit : best-fit distance vs. reddening profile in each pixel.
  • /GRDiagnostic : Gelman-Rubin convergence diagnostic in each pixel.

All four datasets are ordered in the same way, so that the nth element of the /samples dataset corresponds to the same pixel as described by the nth entry in /pixel_info. As our 3D dust map contains pixels of different sizes, /pixel_info specifies each pixel by a HEALPix nside and nested pixel index.

An example in Python will help illustrate the structure of the file:

>>> import numpy as np
>>> import h5py
>>> 
>>> f = h5py.File('dust-map-3d.h5', 'r')
>>> pix_info = f['/pixel_info'][:]
>>> samples = f['/samples'][:]
>>> best_fit = f['/best_fit'][:]
>>> GR = f['/GRDiagnostic'][:]
>>> f.close()
>>> 
>>> print(pix_info['nside'])
[512 512 512 ..., 1024 1024 1024]
>>> 
>>> print(pix_info['healpix_index'])
[1461557 1461559 1461602 ..., 6062092 6062096 6062112]
>>> 
>>> print(pix_info['n_stars'])
[628 622 688 ..., 322 370 272]
>>> 
>>> # Best-fit E(B-V) in each pixel
>>> best_fit.shape  # (# of pixels, # of distance bins)
(2437292, 31)
>>> 
>>> # Get the best-fit E(B-V) in each distance bin for the first pixel
>>> best_fit[0]
array([ 0.00401   ,  0.00554   ,  0.012     ,  0.01245   ,  0.01769   ,
        0.02089   ,  0.02355   ,  0.03183   ,  0.04297   ,  0.08127   ,
        0.11928   ,  0.1384    ,  0.95464998,  0.9813    ,  1.50296998,
        1.55045998,  1.81668997,  1.86567998,  1.9109    ,  2.00281   ,
        2.01739001,  2.02519011,  2.02575994,  2.03046989,  2.03072   ,
        2.03102994,  2.03109002,  2.03109002,  2.03110003,  2.03110003,
        2.03111005], dtype=float32)
>>> 
>>> # Samples of E(B-V) from the Markov Chain
>>> samples.shape  # (# of pixels, # of samples, # of distance bins)
(2437292, 20, 31)
>>> 
>>> # The Gelman-Rubin convergence diagnostic in the first pixel.
>>> # Each distance bin has a separate value.
>>> # Typically, GR > 1.1 indicates non-convergence.
>>> GR[0]
array([ 1.01499999,  1.01999998,  1.01900005,  1.01699996,  1.01999998,
        1.01999998,  1.02400005,  1.01600003,  1.00800002,  1.00600004,
        1.00100005,  1.00199997,  1.00300002,  1.02499998,  1.01699996,
        1.00300002,  1.01300001,  1.00300002,  1.00199997,  1.00199997,
        1.00199997,  1.00199997,  1.00100005,  1.00100005,  1.00100005,
        1.00100005,  1.00100005,  1.00100005,  1.00100005,  1.00100005,
        1.00100005], dtype=float32)

As a simple example of how to work with the full data cube, we will plot the median reddening in the farthest distance bin. We begin by opening the HDF5 file and extracting the information we need:

>>> import numpy as np
>>> import h5py
>>> import healpy as hp
>>> import matplotlib.pyplot as plt
>>> 
>>> # Open the file and extract pixel information and median reddening in the far limit
>>> f = h5py.File('dust-map-3d.h5', 'r')
>>> pix_info = f['/pixel_info'][:]
>>> EBV_far_median = np.median(f['/samples'][:,:,-1], axis=1)
>>> f.close()

The variable pix_info specifies the location of each pixel (by nside and healpix_index), while EBV_far_median contains the median reddening in each pixel in the farthest distance bin. We want to construct a single-resolution HEALPix map, which we can use standard library routines to plot.

We find the maximum nside present in the map, and create an empty array, pix_val, to house the upsampled map:

>>> # Construct an empty map at the highest HEALPix resolution present in the map
>>> nside_max = np.max(pix_info['nside'])
>>> n_pix = hp.pixelfunc.nside2npix(nside_max)
>>> pix_val = np.empty(n_pix, dtype='f8')
>>> pix_val[:] = np.nan

Now, we have to fill the upsampled map, by putting every pixel in the original map into the correct location(s) in the upsampled map. Because our original map has multiple resolutions, pixels that are below the maximum resolution correspond to multiple pixels in the upsampled map. We loop through the nside resolutions present in the original map, placing all the pixels of the same resolution into the upsampled map at once:

>>> # Fill the upsampled map
>>> for nside in np.unique(pix_info['nside']):
...     # Get indices of all pixels at current nside level
...     idx = pix_info['nside'] == nside
...     
...     # Extract E(B-V) of each selected pixel
...     pix_val_n = EBV_far_median[idx]
...     
...     # Determine nested index of each selected pixel in upsampled map
...     mult_factor = (nside_max/nside)**2
...     pix_idx_n = pix_info['healpix_index'][idx] * mult_factor
...     
...     # Write the selected pixels into the upsampled map
...     for offset in range(mult_factor):
...         pix_val[pix_idx_n+offset] = pix_val_n[:]

Now we have an array, pix_val, that represents reddening in the farthest distance bin. We can use one of healpy's built-in visualization functions to plot the map:

>>> # Plot the results using healpy's matplotlib routines
>>> hp.visufunc.mollview(pix_val, nest=True, xsize=4000,
                         min=0., max=4., rot=(130., 0.),
                         format=r'$%g$',
                         title=r'$\mathrm{E} ( B-V \, )$',
                         unit='$\mathrm{mags}$')
>>> plt.show()

Here's the resulting map. Note that it's centered on (ℓ, b) = (130°, 0°):

All code snippets included on this page are covered by the MIT License. In other words, feel free to use the code provided here in your own work.