This notebook provides a user's guide to functions in desiutil.plots for producing all-sky plots. Refer to the docstring for each function for details on its parameters and usage. Operations involving healpix maps require that the healpy package is installed (not available via conda, so use pip instead).
In [1]:
%matplotlib inline
import numpy as np
from matplotlib.projections import get_projection_names
from astropy.coordinates import ICRS
In [2]:
from desiutil.plots import prepare_data, init_sky, plot_grid_map, plot_healpix_map, plot_sky_circles, plot_sky_binned
In [3]:
stars = ICRS(['05h14m32.27s', '05h55m10.29s', '05h25m07.87s', '05h36m12.81s', '05h40m45.52s', '05h47m45.39s', '05h32m00.40s',
'06h45m09.25s', '06h58m37.55s', '07h08m23.49s', '06h22m41.99s', '07h24m05.71s',
'02h31m47.08s', '14h50m42.40s', '15h20m43.75s', '16h45m58.16s', '15h44m03.46s', '17h32m12.90s', '16h17m30.50s',
'12h54m01.63s', '11h03m43.84s', '13h47m32.55s', '13h23m55.54s', '11h01m50.39s', '11h53m49.74s', '12h15m25.45s',],
['−08d12m05.9s', '+07d24m25.3s', '+06d20m59.0s', '−01d12m06.9s', '−01d56m33.3s', '−09d40m10.6s', '−00d17m56.7s',
'-16d42m47.3s', '-28d58m19.5s', '−26d23m35.5s', '−17d57m21.3s', '−29d18m11.2s',
'+89d15m50.9s', '+74d09m19.7s', '+71d50m02.3s', '+82d02m14.1s', '+77d47m40.2s', '+86d35m10.8s', '+75d45m16.9s',
'+55d57m35.4s', '+61d45m04.0s', '+49d18m47.9s', '+54d55m31.3s', '+56d22m56.4s', '+53d41m41.0s', '+57d01m57.4s'])
Use init_sky() to prepare the default Mollweide all-sky projection, which can then be used for other operations like any other axes object:
In [4]:
ax = init_sky()
p = ax.scatter(ax.projection_ra(stars.ra.degree), ax.projection_dec(stars.dec.degree), marker='.', color='black')
By default, the map is oriented so the DESI footprint does not wrap around in RA and shows the galactic plane in red, but you can change most of this with options:
In [5]:
ax = init_sky(projection='hammer', ra_center=0, galactic_plane_color=None)
p = ax.scatter(ax.projection_ra(stars.ra.degree), ax.projection_dec(stars.dec.degree), marker='.', color='black')
To see the available projections, use matplotlib.projections.get_projection_names(). The code in this package will probably work best with 'mollweide', 'hammer' or 'aitoff', where 'mollweide' is the default.
In [6]:
get_projection_names()
Out[6]:
This section describes two methods that interpret an array of values as a map on the sphere and create a corresponding plot. If you have a list of (RA, Dec) coordinates that you wish to plot, instead of a map, refer to the next section on catalog plotting.
Use plot_grid_map() to map a 2D array of values tabulated on a grid of RA and DEC values, with RA indexed by the first index. In addition to the grid values, the grid edges must also be specified. Since values are associated with the cells between grid lines, a grid with dimensions $N_{\mathrm{RA}}, N_{\mathrm{Dec}}$ has corresponding edge arrays of length $N_{\mathrm{RA}} + 1$ and $N_{\mathrm{Dec}} + 1$. The grid must wrap around the sphere in RA but does not need to cover the full range of Dec. The grid can be chosen independently of the central RA of the map projection and any wrap-around grid column will be split automatically.
In [7]:
np.random.seed(20200415) # For reproducible maps.
n_ra, n_dec = 15, 10
ra_edges = np.linspace(0., 360., n_ra + 1)
dec_edges = np.linspace(-60., 60., n_dec + 1)
data = np.random.uniform(size=(n_dec, n_ra))
ax = plot_grid_map(data, ra_edges, dec_edges, label='Value',
galactic_plane_color=None, ecliptic_plane_color=None)
Use plot_map_healpix() to map a 1D array of values associated with a healpix map. The array size must exactly match the expected size for some value of the healpix NSIDE parameter. Each pixel is rendered as a matplotlib polygon, so this function is somewhat slower than plot_grid_map() and also does not automatically wrap around pixels at the RA edges (which are not drawn at all, but this is less visible at larger values of NSIDE).
In [8]:
import healpy as hp
In [9]:
nside, nest = 4, False
data = np.random.uniform(size=hp.nside2npix(nside))
ax = plot_healpix_map(data, nest, galactic_plane_color=None, ecliptic_plane_color=None)
Use the colorbar option to suppress the default colorbar, or the label option to add a colorbar label.
This section describes two methods for plotting a catalog of objects indexed by (RA, DEC).
Use plot_sky_circles() to draw a circular field of view centered on each (RA, DEC). This function is useful in DESI for displaying observing tiles and the default field of view is set accordingly.
In [10]:
ax = plot_sky_circles(ra_center=[0, 30, 60, 90, 120, 150, 180, 210, 240],
dec_center=[-60, -45, -30, -15, 0, 15, 30, 45, 60],
field_of_view=10)
Note that the projection of a circular field of view is generally not circular. Face colors can be specified for each circle independently.
In [11]:
ax = plot_sky_circles(ra_center=[60, 90, 120, 150],
dec_center=[-30, -15, 0, 15],
field_of_view=30,
facecolors=['r', 'g', 'b', 'k'])
Each (RA, Dec) can optionally be associated with some value that is used to determine the color of each circle.
In [12]:
ax = plot_sky_circles(ra_center=[0, 60, 120, 180],
dec_center=[-30, 0, 30, 60],
field_of_view=30,
data=[0, 1, 2, 3], cmap='gray', edgecolor='r');
Use plot_sky_binned() to automatically histogram object locations and display the resulting pixel counts per unit sky area as a map, using either the default plot_grid_map() or plot_healpix_map() (which requires that healpy is installed).
In [13]:
n = 50000
ra = 360 * np.random.uniform(size=n)
dec = np.clip(20 * np.random.normal(size=n), -90, +90)
ax = plot_sky_binned(ra, dec)
In [14]:
ax = plot_sky_binned(ra, dec, plot_type='healpix')
The automatic histogramming is controlled by the value (in sq.deg.) of the max_bin_area parameter. The default grid map is binned in cos(Dec) and RA and aims for roughly square pixels on the sky. Use the verbose option to print details of the automatic histogramming. Note that a healpix map is more constrained so will generally result in smaller pixels for the same value of max_bin_area (and corrrespondingly larger shot noise).
In [15]:
ax = plot_sky_binned(ra, dec, max_bin_area=50, verbose=True)
In [16]:
ax = plot_sky_binned(ra, dec, plot_type='healpix', max_bin_area=50, verbose=True)
Instead of plotting object densities, you can provide an array of data associated with each object and the resulting map will display the mean data value per pixel. Note that empty bins, where the mean data is undefined, are automatically masked.
In [17]:
ax = plot_sky_binned(ra, dec, data=ra)
Here's a more useful example of a "source" on the sky. Note that we have to prepare the RA values to deal with wrap-around.
In [18]:
ra_source = np.random.normal(loc=0.0, scale=5.0, size=(10000,))
dec_source = np.random.normal(loc=30.0, scale=5.0, size=(10000,))
ra_source[ra_source < 0] += 360.0
ax = plot_sky_binned(ra_source, dec_source, max_bin_area=2.0, verbose=True)
These plotting routines all pass the data to be displayed through prepare_data() in order to provide a uniform API for how plot values are masked and mapped to display colors.
In [19]:
n_ra, n_dec = 15, 10
ra_edges = np.linspace(0., 360., n_ra + 1)
dec_edges = np.linspace(-60., 60., n_dec + 1)
data = np.random.uniform(size=(n_dec, n_ra))
ax = plot_grid_map(data, ra_edges, dec_edges)
For more control, you can call prepare_data() yourself before calling the map functions. For example, you can explicitly mask some pixels in a map:
In [20]:
data2 = prepare_data(data, mask=(0.4 < data) & (data < 0.6))
ax = plot_grid_map(data2, ra_edges, dec_edges)
You can also specify the clipping to apply at both ends (lo/hi) using a mixture of percentiles or absolute values.
In [21]:
data2 = prepare_data(data, clip_lo='20%', clip_hi=0.8)
ax = plot_grid_map(data2, ra_edges, dec_edges)
Values outside the clip range are normally displayed (after clipping) but can instead be masked by adding '!' to the clip specification. Note that quotes are optional for numeric clip values (0.8) but required when using '%' or '!'.
In [22]:
data2 = prepare_data(data, clip_lo='!20%', clip_hi='!0.8')
ax = plot_grid_map(data2, ra_edges, dec_edges)
If the range of clipped data is smaller than the clip limits, the tighter limits will be used for the color scale by default:
In [23]:
data2 = prepare_data(data, clip_lo=-0.5, clip_hi=1.5)
ax = plot_grid_map(data2, ra_edges, dec_edges)
However, you can force the color scale to use either clip limit by passing save_limits=True. This is useful when making a series of plots with a consistent colorscale.
In [24]:
data2 = prepare_data(data, clip_lo=-0.5, clip_hi=1.5, save_limits=True)
ax = plot_grid_map(data2, ra_edges, dec_edges)
The same techniques can be used with the data passed to plot_sky_circles().
In [25]:
data = prepare_data([0, 1, 2, 3], mask=[False, True, False, False])
ax = plot_sky_circles(ra_center=[0, 60, 120, 180],
dec_center=[-30, 0, 30, 60],
field_of_view=30, data=data, edgecolor='k')
The plot_sky_binned() function calculates the data to plot internally, but has parameters clip_lo and clip_hi that are used to prepare this data.
In [26]:
ax = plot_sky_binned(ra, dec, data=ra, clip_hi='!180')
In [ ]: