301.4. Extended Chandra Deep Field South (ECDFS)

301.4. Extended Chandra Deep Field South (ECDFS)¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: DP1
Container Size: large
LSST Science Pipelines version: r29.1.1
Last verified to run: 06-24-2025
Repository: github.com/lsst/tutorial-notebooks

Learning objective: To explore the Extended Chandra Deep Field South (ECDFS) target field from commissioning with LSSTComCam, using Data Preview 1 (DP1) data.

LSST data products: deep_coadd, CcdVisit, Visit, Object

Packages: lsst.rsp, lsst.daf.butler, lsst.afw.display

Credit: Originally developed by the Rubin Community Science team. Please consider acknowledging them if this notebook is used for the preparation of journal articles, software releases, or other notebooks.

Get Support: Everyone is encouraged to ask questions or raise issues in the Support Category of the Rubin Community Forum. Rubin staff will respond to all questions posted there.

1. Introduction¶

This notebook examines the Data Preview 1 (DP1) data products in the Extended Chandra Deep Field South (ECDFS), including magnitude limits, visit distribution with time, data quality, and the distributions of stars and galaxies in color-magnitude and color-color diagrams.

The Extended Chandra Deep Field South (ECDFS) is a popular extragalatic deep field that has attracted a wealth of panchromatic imaging and spectroscopy. The Rubin map of ECDFS is centered at (RA, Dec) = (53.16 -28.10) The region covered spans a diameter of about 1 degree.

Related tutorials: The 100-level tutorials demonstrate how to use the butler, the Table Access Protocol (TAP) service, and the Firefly image display. The 200-level tutorials introduce the types of image and catalog data.

1.1. Import packages¶

Import numpy, a fundamental package for scientific computing with arrays in Python (numpy.org), and matplotlib, a comprehensive library for data visualization (matplotlib.org; matplotlib gallery), including custom shapes (Polygon) and lines (mlines). itertools supports efficient iteration and combinatorics.

From the lsst package, import modules for accessing the Table Access Protocol (TAP) service, for retrieving datasets from the butler, and for image displaying from the LSST Science Pipelines (pipelines.lsst.io). Additional modules support spherical geometry (sphgeom), 2D geometry (geom), and standardized multiband plotting (lsst.utils.plotting) for LSST data analysis and visualization.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.lines as mlines
from matplotlib.patches import Polygon
import itertools

from lsst.rsp import get_tap_service
from lsst.daf.butler import Butler
import lsst.afw.display as afw_display
import lsst.sphgeom as sphgeom
import lsst.geom as geom

from lsst.utils.plotting import (get_multiband_plot_colors,
                                 get_multiband_plot_symbols,
                                 get_multiband_plot_linestyles)

1.2. Define parameters and functions¶

Create an instance of the TAP service, and assert that it exists.

In [2]:

service = get_tap_service("tap")
assert service is not None

Create an instance of the butler, and assert that it exists.

In [3]:

butler = Butler('dp1', collections="LSSTComCam/DP1")
assert butler is not None

Define the approximate central coordinates of the ECDFS in degrees, and the region to use when querying for overlapping images.

In [4]:

target_ra = 53.2
target_dec = -28.1

ra_cen = target_ra
dec_cen = target_dec

radius = 1.0
region = sphgeom.Region.from_ivoa_pos(f"CIRCLE {target_ra} {target_dec} {radius}")

Define parameters to use colorblind-friendly colors with matplotlib.

In [5]:

plt.style.use('seaborn-v0_8-colorblind')
prop_cycle = plt.rcParams['axes.prop_cycle']
colors = prop_cycle.by_key()['color']

Define colors, symbols, and linestyles to represent the six LSST filters, $ugrizy$.

In [6]:

filter_colors = get_multiband_plot_colors()
filter_symbols = get_multiband_plot_symbols()
filter_linestyles = get_multiband_plot_linestyles()

Set afwDisplay to use Firefly to display images, and open Firefly frame 1.

In [7]:

afw_display.setDefaultBackend('firefly')
display = afw_display.Display(frame=1)

2. Coadded coverage¶

The individual LSSTLSSTComCam images have been combined (stacked) into deep_coadd images, and sky maps of the accumulated exposure time and image depth have been created.

2.1. Display a deep coadd¶

Query the butler for all $r$-band deep_coadd images that overlap the defined region.

In [8]:

coadd_datasetrefs = butler.query_datasets("deep_coadd",
                                          where="patch.region OVERLAPS region AND band='r'",
                                          bind={"region": region},
                                          with_dimension_records=True,
                                          order_by=["patch.tract"])

In [9]:

print('There are %s r-band deep_coadd images.' % len(coadd_datasetrefs))

There are 82 r-band deep_coadd images.

Retrieve the deep_coadd image associated with one of the dataset references returned by the butler query.

In [10]:

coadd = butler.get(coadd_datasetrefs[10])

Display the image, and set the mask plane to be fully transparent.

In [11]:

display.mtv(coadd)
display.setMaskTransparency(100)

2.2. Exposure time and depth¶

Retrieve the entire survey property map of the $r$-band magnitude limit as hspmap_rmaglim, and of the $r$-band exposure time as hspmap_rexptime.

In [12]:

hspmap_rmaglim = butler.get('deepCoadd_psf_maglim_consolidated_map_weighted_mean',
                            skymap='lsst_cells_v1', band='r')
hspmap_rexptime = butler.get('deepCoadd_exposure_time_consolidated_map_sum',
                             skymap='lsst_cells_v1', band='r')

Extract the healsparse map values within the radius of the field center and divide the region into a 250×250 grid. Evaluate each survey property map at these grid points. Replace all negative values (no-data placeholders) with NaN to enable proper map display.

In [13]:

dec_size = radius
ra_size = dec_size / np.cos(np.radians(dec_cen))
ra_min, ra_max = ra_cen - ra_size, ra_cen + ra_size
dec_min, dec_max = dec_cen - dec_size, dec_cen + dec_size

ra = np.linspace(ra_min, ra_max, 250)
dec = np.linspace(dec_min, dec_max, 250)
x, y = np.meshgrid(ra, dec)

values_rmaglim = hspmap_rmaglim.get_values_pos(x, y)
values_rmaglim = np.where(values_rmaglim < 0.0, np.nan, values_rmaglim)

values_rexptime = hspmap_rexptime.get_values_pos(x, y)
values_rexptime = np.where(values_rexptime < 0.0, np.nan, values_rexptime)

Clean up.

In [14]:

del ra, dec

Print the mean, median, and min/max values for the $r$-band magnitude limit and exposure time. Exposure times are accumulated over multiple visits, not single exposures.

In [15]:

print('r-band magnitude limit mean/median: %5.2f %5.2f' % (np.nanmean(values_rmaglim),
                                                           np.nanmedian(values_rmaglim)))
print('r-band magnitude min/max: %5.2f %5.2f' % (np.nanmin(values_rmaglim),
                                                 np.nanmax(values_rmaglim)))
print('\n')
print('r-band exposure time (s) mean/median: %5.1f %5.1f' % (np.nanmean(values_rexptime),
                                                             np.nanmedian(values_rexptime)))
print('r-band exposure time min/max (s): %5.1f %5.1f' % (np.nanmin(values_rexptime),
                                                         np.nanmax(values_rexptime)))

r-band magnitude limit mean/median: 26.07 26.34
r-band magnitude min/max: 21.46 27.12


r-band exposure time (s) mean/median: 2433.5 1650.0
r-band exposure time min/max (s):  30.0 6870.0

2.3. Coverage map¶

Define unique colors and linestyles for plotting each of the tracts in the dataset.

In [16]:

alltracts = [rec.dataId['tract'] for rec in coadd_datasetrefs]
unique_tracts = np.unique(alltracts)

colors = plt.cm.tab10(np.linspace(0, 1, len(unique_tracts)))
linestyles = ['-', '--', '-.', ':'] * ((len(unique_tracts) // 4) + 1)

style_dict = {
    tract: {'color': colors[i], 'linestyle': linestyles[i]}
    for i, tract in enumerate(unique_tracts)
}

The coadd_datasetrefs returned from the butler include region information for each dataset, stored as a ConvexPolygon3D with four vertices defining the sky footprint of each overlapping tract and patch.

Convert these vertices to 2D sky coordinates (RA, Dec) using geom.SpherePoint, and use matplotlib.patches.Polygon along with ax.add_patch() to draw each patch outline. Each tract is plotted with a distinct color and linestyle for visual clarity.

In [17]:

fig, ax = plt.subplots(figsize=(6, 5))

mesh = ax.pcolormesh(x, y, values_rmaglim, cmap='Greys_r', shading='auto')
fig.colorbar(mesh, ax=ax, label="r-band limiting magnitude (mag)")

for rec in coadd_datasetrefs:
    vertices = rec.dataId.patch.region.getVertices()
    vertices_deg = []
    for vertex in vertices:
        vertices_deg.append([geom.SpherePoint(vertex).getRa().asDegrees(),
                             geom.SpherePoint(vertex).getDec().asDegrees()])
    polygon = Polygon(vertices_deg, closed=True, facecolor='None',
                      edgecolor=style_dict[rec.dataId['tract']]['color'],
                      linestyle=style_dict[rec.dataId['tract']]['linestyle'],
                      linewidth=2)
    ax.add_patch(polygon)

ax.set_xlim(ra_max, ra_min)
ax.set_ylim(dec_min, dec_max)
ax.set_xlabel('RA (deg)')
ax.set_ylabel('Dec (deg)')

handles = [
    mlines.Line2D([], [], color=style['color'], linestyle=style['linestyle'],
                  label=f"Tract {tract}")
    for tract, style in style_dict.items()
]
ax.legend(handles=handles, loc='upper left', ncol=2)

ax.minorticks_on()
ax.set_title('Tracts and patches overlapping ECDFS')

plt.tight_layout()
plt.show()

Figure 1: For the region of the ECDFS, the $r$-band magnitude limit survey property map in greyscale, overplotted with the patch boundaries colored by their tract number (as in the legend).

Clean up.

In [18]:

del coadd_datasetrefs
del hspmap_rmaglim, hspmap_rexptime
del x, y, values_rmaglim, values_rexptime
del style_dict

3. Visits (observations)¶

Retrieve all visits (individual observations of the sky) that fall within a circular region centered at (RA, Dec) = (ra_cen, dec_cen) with a radius of 1 degree. Return the visit ID, band, observation midpoint time in both MJD and calendar date.

In [19]:

query = """
        SELECT visit, band, expMidptMJD, expMidpt
        FROM dp1.Visit
        WHERE CONTAINS(POINT('ICRS', ra, dec),CIRCLE('ICRS', {}, {}, {}))=1
        ORDER BY expMidptMJD
        """.format(ra_cen, dec_cen, radius)

In [20]:

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)
if job.phase == 'ERROR':
    job.raise_if_error()

Job phase is COMPLETED

In [21]:

assert job.phase == 'COMPLETED'
visits = job.fetch_result().to_table()
print(f"Total number of visits: {len(visits)}")

Total number of visits: 855

Option to display the table of results.

In [22]:

# visits

3.1. Filter distribution¶

Count and list the number of visits per filter.

In [23]:

all_bands = np.array([visit['band'] for visit in visits], dtype=str)
unique_bands, counts = np.unique(all_bands, return_counts=True)

for band, count in zip(unique_bands, counts):
    print(band, count)

g 230
i 162
r 237
u 43
y 30
z 153

In [24]:

field_filters = ['u', 'g', 'r', 'i', 'z', 'y']

Clean up.

In [25]:

del all_bands, unique_bands, counts

3.2. Visit dates¶

Plot two characteristics of the visit dates:

1. Plot the cumulative distribution of visit dates for each of the filters observed for this field.
2. Plot the distribution of the time separation between visits for each band and for all bands combined.

In [26]:

t_start = visits[0]['expMidptMJD'] - 1.0
t_end = visits[-1]['expMidptMJD'] + 1.0
time_bins = np.arange(t_end - t_start, dtype='int') + t_start
diff_bins = np.arange(0, np.max(time_bins) - np.min(time_bins) + 3, 1)
all_time_diffs = []

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4))

for band in field_filters:
    mjds = sorted([
        visit['expMidptMJD']
        for visit in visits
        if visit['band'] == band
    ])
    mjds = np.array(mjds, dtype=float)
    label = f"{band} ({len(mjds)})"

    counts, _, patches = ax1.hist(
        mjds, bins=time_bins, cumulative=True, histtype='step',
        linewidth=2, color=filter_colors[band], label=label
    )
    for patch in patches:
        patch.set_linestyle(filter_linestyles[band])

    tdiffs = []
    for comb in itertools.combinations(mjds, 2):
        tdiffs.append(comb[1] - comb[0])
    time_diffs = np.array(tdiffs)
    all_time_diffs.append(time_diffs)

    counts_diff, _, patches_diff = ax2.hist(
        time_diffs, bins=diff_bins, histtype='step', linewidth=2,
        color=filter_colors[band], label=band
    )
    for patch in patches_diff:
        patch.set_linestyle(filter_linestyles[band])

all_time_diffs = np.concatenate(all_time_diffs)
ax2.hist(
    all_time_diffs, bins=diff_bins, histtype='step', linewidth=2,
    color='black', label='all filters'
)

ax1.set_xlabel('Modified Julian Date')
ax1.set_ylabel('Cumulative number of visits')
ax1.minorticks_on()

ax2.set_xlabel('Separation between observations (d)')
ax2.set_ylabel('Number of observations')
ax2.set_yscale('log')
ax2.legend(loc='upper right', ncol=2)
ax2.minorticks_on()

plt.show()

Figure 2: The left panel shows the cumulative number of visits per filter as a function of time. The right panel shows histograms of visit-to-visit time separations for each filter individually (colors), and all together (black).

Clean up.

In [27]:

del time_bins, all_time_diffs

3.3. Visit image quality¶

Retrieve statistics about the image quality of the individual visit_images from the CcdVisit table.

In [28]:

query = """
        SELECT visitId, ra, dec, band, magLim, seeing
        FROM dp1.CcdVisit
        WHERE CONTAINS(POINT('ICRS', ra, dec),CIRCLE('ICRS', {}, {}, {}))=1
        ORDER BY visitId
        """.format(ra_cen, dec_cen, radius)

In [29]:

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)
if job.phase == 'ERROR':
    job.raise_if_error()

Job phase is COMPLETED

In [30]:

assert job.phase == 'COMPLETED'
ccd_visits = job.fetch_result().to_table()

Plot the distribution of seeing, which is the mean full width at half maximum of the PSF, for all visits and detector combinations.

In [31]:

min_seeing = np.min(ccd_visits['seeing'])
max_seeing = np.max(ccd_visits['seeing'])
seeing_bins = np.arange(min_seeing, max_seeing, 0.05)

fig = plt.figure(figsize=(7, 5))
for band in field_filters:
    print(f"Median in {band}-band: \
            {np.ma.median(ccd_visits[ccd_visits['band'] == band]['seeing']): .2f}")
    n, bins, patches = plt.hist(ccd_visits[ccd_visits['band'] == band]['seeing'],
                                seeing_bins, histtype='step',
                                linewidth=2,
                                color=filter_colors[band],
                                label=band)
    for patch in patches:
        patch.set_linestyle(filter_linestyles[band])

plt.legend(loc='upper right')
plt.xlabel('PSF FWHM (arcsec)')
plt.ylabel('Number of visits')
plt.minorticks_on()
plt.show()

Median in u-band:              1.40
Median in g-band:              1.14
Median in r-band:              1.08
Median in i-band:              1.00
Median in z-band:              1.00
Median in y-band:              1.07

Figure 3: Histogram of delivered image quality of detector images, per filter.

Clean up.

In [32]:

del seeing_bins

Plot the distribution of the magnitude limit, magLim.

In [33]:

min_mag = np.min(ccd_visits['magLim'])
max_mag = np.max(ccd_visits['magLim'])
maglim_bins = np.arange(min_mag, max_mag, 0.1)

fig = plt.figure(figsize=(7, 5))

for band in field_filters:
    print(f"Median in {band}-band: \
            {np.ma.median(ccd_visits[ccd_visits['band'] == band]['magLim']): .2f}")
    n, bins, patches = plt.hist(ccd_visits[ccd_visits['band'] == band]['magLim'],
                                maglim_bins, histtype='step',
                                linewidth=2,
                                color=filter_colors[band],
                                label=band)
    for patch in patches:
        patch.set_linestyle(filter_linestyles[band])

plt.legend(loc='upper left')
plt.xlabel('limiting magnitude')
plt.ylabel('Number of visits')
plt.minorticks_on()
plt.show()

Median in u-band:              23.48
Median in g-band:              24.68
Median in r-band:              24.38
Median in i-band:              23.97
Median in z-band:              23.15
Median in y-band:              21.93

Figure 4: Histogram of limiting magnitude at signal-to-noise=5 on detector images, per filter.

Clean up.

In [34]:

del maglim_bins

4. Object (detections)¶

The Object table contains forced measurements on the deep coadded images at the locations of all objects detected with signal-to-noise ratio > 5 in a deep_coadd of any filter.

Query the Object table for objects in the field, and return their PSF and cModel AB magnitudes.

In [35]:

query = """SELECT objectId, coord_ra, coord_dec, ebv,
        {}_psfMag, {}_cModelMag, {}_psfMag, {}_cModelMag,
        {}_psfMag, {}_cModelMag, {}_psfMag, {}_cModelMag,
        {}_psfMag, {}_cModelMag, {}_psfMag, {}_cModelMag,
        refExtendedness
        FROM dp1.Object
        WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec),
              CIRCLE('ICRS', {}, {}, {})) = 1
        ORDER BY objectId ASC
        """.format(field_filters[0], field_filters[0], field_filters[1], field_filters[1],
                   field_filters[2], field_filters[2], field_filters[3], field_filters[3],
                   field_filters[4], field_filters[4], field_filters[5], field_filters[5],
                   ra_cen, dec_cen, radius)
print(query)

SELECT objectId, coord_ra, coord_dec, ebv,
        u_psfMag, u_cModelMag, g_psfMag, g_cModelMag,
        r_psfMag, r_cModelMag, i_psfMag, i_cModelMag,
        z_psfMag, z_cModelMag, y_psfMag, y_cModelMag,
        refExtendedness
        FROM dp1.Object
        WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec),
              CIRCLE('ICRS', 53.2, -28.1, 1.0)) = 1
        ORDER BY objectId ASC

In [36]:

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)
if job.phase == 'ERROR':
    job.raise_if_error()

Job phase is COMPLETED

In [37]:

assert job.phase == 'COMPLETED'
objtab = job.fetch_result().to_table()

4.1. Magnitude distributions¶

Plot histograms of object magnitudes, separating the samples into stars and galaxies.

Use the refExtendedness flag to distinguish them: treat objects with refExtendedness == 1 as likely galaxies (extended) and those with refExtendedness == 0 as likely stars (point sources).

Use cModelMag mags for galaxies and psfMag for stars.

In [38]:

ptsource = (objtab['refExtendedness'] == 0)

In [39]:

mag_bins = np.arange(15.4, 28.6, 0.2)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 4), sharey=True)
plt.subplots_adjust(hspace=0, wspace=0)

for band in field_filters:
    mag_psf = objtab[f"{band}_psfMag"]
    mag_cmodel = objtab[f"{band}_cModelMag"]

    valid_mags = (15 < mag_psf) & (mag_psf < 30) & (15 < mag_cmodel) & (mag_cmodel < 30)

    print(f"Number of objects in {band}-band: {np.sum(valid_mags)}")

    n_psf, bins_psf, patches_psf = ax1.hist(mag_psf[valid_mags & ptsource], bins=mag_bins,
                                            histtype='step', linewidth=2,
                                            color=filter_colors[band], label=band,)
    for patch in patches_psf:
        patch.set_linestyle(filter_linestyles[band])

    n_cmodel, bins_cmodel, patches_cmodel = ax2.hist(
        mag_cmodel[valid_mags & ~ptsource],
        bins=mag_bins,
        histtype='step',
        linewidth=2,
        color=filter_colors[band],
        label=band,
    )
    for patch in patches_cmodel:
        patch.set_linestyle(filter_linestyles[band])

ax1.set_xlim(mag_bins.min(), mag_bins.max())
ax1.set_yscale('log')
ax1.set_xlabel('Magnitude')
ax1.set_ylabel('Number of objects')
ax1.set_title('Point sources (PSF mags)')
ax1.legend(loc='upper left', ncols=2)
ax1.minorticks_on()

ax2.set_xlim(mag_bins.min(), mag_bins.max())
ax2.set_xlabel('Magnitude')
ax2.set_title('Extended sources (cModel mags)')
ax2.minorticks_on()

plt.show()

Number of objects in u-band: 368728
Number of objects in g-band: 442595
Number of objects in r-band: 437816
Number of objects in i-band: 440706
Number of objects in z-band: 424969
Number of objects in y-band: 341877

Figure 5: Histogram of object magnitudes separated by morphological type. The left panel shows likely stars (refExtendedness = 0) using PSF magnitudes, while the right panel shows likely galaxies cModel magnitudes. The bump near 19 mag in the stellar distribution corresponds to red clump stars in the SMC.

Clean up.

In [40]:

del mag_bins, mag_psf, mag_cmodel, valid_mags, n_psf, bins_psf, patches_psf
del n_cmodel, bins_cmodel, patches_cmodel

5. Analysis¶

A short demonstration of some standard scientific analyses.

5.1. Color magnitude diagram¶

Plot $r$ vs. $(g - r)$ CMDs and $(r - i)$ vs. $(g - r)$ CCDs. Create separate plots for stars and galaxies using PSF magnitudes for stars and cModel magnitudes for galaxies.

In [41]:

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(7, 6),
                                             height_ratios=[1, 2.5])
plt.subplots_adjust(hspace=0, wspace=0)

grmin, grmax = -0.9, 2.3
rimin, rimax = -1.3, 2.8
magmin, magmax = 15.8, 26.8

ax1.hexbin(objtab[ptsource]['g_psfMag']-objtab[ptsource]['r_psfMag'],
           objtab[ptsource]['r_psfMag']-objtab[ptsource]['i_psfMag'],
           gridsize=(150, 150),
           extent=(grmin, grmax, rimin, rimax), bins='log', cmap='Grays')
ax1.set_title('Point sources (PSF mags)')
ax1.set_xlim(grmin, grmax)
ax1.set_ylim(rimin, rimax)
ax1.set_ylabel(r'$(r-i)$')
ax1.set_xticklabels([])
ax1.minorticks_on()

ax2.hexbin(objtab[~ptsource]['g_cModelMag']-objtab[~ptsource]['r_cModelMag'],
           objtab[~ptsource]['r_cModelMag']-objtab[~ptsource]['i_cModelMag'],
           gridsize=(150, 150),
           extent=(grmin, grmax, rimin, rimax), bins='log', cmap='Grays')
ax2.set_title('Extended sources (cModel mags)')
ax2.set_xlim(grmin, grmax)
ax2.set_ylim(rimin, rimax)
ax2.set_xticklabels([])
ax2.set_yticklabels([])
ax2.minorticks_on()

ax3.hexbin(objtab[ptsource]['g_psfMag']-objtab[ptsource]['r_psfMag'],
           objtab[ptsource]['r_psfMag'], gridsize=(150, 200),
           extent=(grmin, grmax, magmin, magmax), bins='log', cmap='Grays')

ax3.set_xlim(grmin, grmax)
ax3.set_ylim(magmax, magmin)
ax3.set_xlabel(r'$(g-r)$')
ax3.set_ylabel(r'$r$ magnitude')
ax3.minorticks_on()

ax4.hexbin(objtab[~ptsource]['g_cModelMag']-objtab[~ptsource]['r_cModelMag'],
           objtab[~ptsource]['r_cModelMag'], gridsize=(150, 200),
           extent=(grmin, grmax, magmin, magmax), bins='log', cmap='Grays')
ax4.set_xlim(grmin, grmax)
ax4.set_ylim(magmax, magmin)
ax4.set_xlabel(r'$(g-r)$')
ax4.set_yticklabels([])
ax4.minorticks_on()

plt.show()

Figure 6: CMDs and CCDs for objects classified by morphology. The left column shows likely stars (refExtendedness = 0) using PSF magnitudes, and the right column shows likely galaxies using cModel magnitudes. Top panels display $(r - i)$ vs. $(g - r)$ CCDs; bottom panels show $r$ vs. $(g - r)$ CMDs.

5.2. Galaxy evolution analyses¶

This section demonstrates a few tasks more specific to galaxy evolution analyses since the ECDFS is a popular well-studied extragalactic field. Two aspects are performed below:

1. First zoom in on the hubble ultra-deep field as a visual comparison to complementary datasets.
2. Color selection of high-redshift galaxies at redshift z~4, using a standard Lyman-break selection technique (optionally compare to some literature coordinates?)

First, to enable easier comparison to other data set the afw backend to matplotlib.

In [42]:

afw_display.setDefaultBackend('matplotlib')

Below, define the central coordinates of the Hubble Ultra Deep Field South (e.g. Beckwith et al. 2006; Ellis et al. 2013; Koekemoer et al. 2013). Then, identify the data ID of ECDFS deep_coadd r-band image that contains this region of the sky.

In [43]:

hudf_ra = 53.159
hudf_dec = -27.7815854

center_coadd_datasetrefs = butler.query_datasets("deep_coadd",
                                                 where="band.name='r' AND\
                                                        patch.region OVERLAPS POINT(ra, dec)",
                                                 bind={"ra": hudf_ra, "dec": hudf_dec},
                                                 with_dimension_records=True,
                                                 order_by=["patch.tract"])


coadd = butler.get(center_coadd_datasetrefs[0])
dataId = center_coadd_datasetrefs[0].dataId

The cell below will ask the Butler to generate a cutout that is 2500 LSST pixels on a side (8.3 x 8.3 arcmin) centered on the HUDF. Compare to a cutout of the same part of the sky from Hubble Space Telescope F606W imaging (similar effective wavelength to LSST r-band).

In [44]:

wcs = butler.get("deep_coadd.wcs", dataId=dataId)
radec_hudf = geom.SpherePoint(hudf_ra, hudf_dec, geom.degrees)
xy_hudf = geom.PointI(wcs.skyToPixel(radec_hudf))

cutoutSize = geom.ExtentI(2500, 2500)
bbox = geom.BoxI(xy_hudf - cutoutSize // 2, cutoutSize)
parameters = {'bbox': bbox}
cutout_image = butler.get('deep_coadd', parameters=parameters, dataId=dataId)

fig, ax = plt.subplots()
display = afw_display.Display(frame=fig)
display.scale('asinh', 'zscale')
display.mtv(cutout_image.image)
plt.show()

Figure 6: Cutout of the region of ECDFS containing the Hubble Ultra Deep Field.

Figure 7: a screenshot of a Hubble Space Telescope V-band (F606W; 0.6 micron wavelength) image of the same region.

5.3. Select candidate high-redshift galaxies¶

The cell below demonstrates the color selection of high-redshift galaxies in the LSST ECFDS deep_coadds by identifying galaxies with a strong lyman-break at redshift~4 (g-band dropouts). The colors defined below reference the selection defined in Overzier et al. 2006 which uses similar filtersets to the LSST filters.

First, query the Object Table to return cModelFlux values u,g,r,i,z filters for galaxies (r_extendedness=1) that are well detected (with signal to noise more than 20 in the r and i bands).

In [45]:

query_lines = [
    "SELECT objectId, coord_ra, coord_dec, i_extendedness, "
]

for filt in field_filters:
    query_lines.append(f"       {filt}_cModelFlux, {filt}_cModelFluxErr, ")
    query_lines.append(f"       {filt}_cModelMag AS {filt}mag, "
                       f"{filt}_cModelMagErr AS {filt}magErr, ")

query_lines.append("       refExtendedness")
query_lines.append("FROM dp1.Object")
query_lines.append(
    f"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), "
    f"CIRCLE('ICRS', {ra_cen}, {dec_cen}, {radius})) = 1 "
    f"AND r_cModelFlux/r_cModelFluxErr > 20 "
    f"AND i_cModelFlux/i_cModelFluxErr > 20"
    f"AND r_extendedness = 1"
)
query_lines.append("ORDER BY objectId ASC")

query = "\n".join(query_lines)
print(query)

SELECT objectId, coord_ra, coord_dec, i_extendedness,
       u_cModelFlux, u_cModelFluxErr,
       u_cModelMag AS umag, u_cModelMagErr AS umagErr,
       g_cModelFlux, g_cModelFluxErr,
       g_cModelMag AS gmag, g_cModelMagErr AS gmagErr,
       r_cModelFlux, r_cModelFluxErr,
       r_cModelMag AS rmag, r_cModelMagErr AS rmagErr,
       i_cModelFlux, i_cModelFluxErr,
       i_cModelMag AS imag, i_cModelMagErr AS imagErr,
       z_cModelFlux, z_cModelFluxErr,
       z_cModelMag AS zmag, z_cModelMagErr AS zmagErr,
       y_cModelFlux, y_cModelFluxErr,
       y_cModelMag AS ymag, y_cModelMagErr AS ymagErr,
       refExtendedness
FROM dp1.Object
WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), CIRCLE('ICRS', 53.2, -28.1, 1.0)) = 1 AND r_cModelFlux/r_cModelFluxErr > 20 AND i_cModelFlux/i_cModelFluxErr > 20AND r_extendedness = 1
ORDER BY objectId ASC

In [46]:

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)
if job.phase == 'ERROR':
    job.raise_if_error()

Job phase is COMPLETED

In [47]:

assert job.phase == 'COMPLETED'
results = job.fetch_result()
results.to_table()

Out[47]:

Table length=74388

objectId	coord_ra	coord_dec	i_extendedness	u_cModelFlux	u_cModelFluxErr	umag	umagErr	g_cModelFlux	g_cModelFluxErr	gmag	gmagErr	r_cModelFlux	r_cModelFluxErr	rmag	rmagErr	i_cModelFlux	i_cModelFluxErr	imag	imagErr	z_cModelFlux	z_cModelFluxErr	zmag	zmagErr	y_cModelFlux	y_cModelFluxErr	ymag	ymagErr	refExtendedness
	deg	deg		nJy	nJy	mag	mag	nJy	nJy	mag	mag	nJy	nJy	mag	mag	nJy	nJy	mag	mag	nJy	nJy	mag	mag	nJy	nJy	mag	mag
int64	float64	float64	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32
609780077793644595	52.98446522868694	-28.69956236608048	1.0	2143.64	427.446	23.0721	0.216498	4107.5	74.7722	22.3661	0.0197645	11146.1	239.547	21.2822	0.0233342	20625.0	228.074	20.614	0.0120062	27612.3	370.524	20.2972	0.0145693	--	--	--	--	1.0
609780077793644596	52.98620407898802	-28.693789376693445	1.0	847.695	459.305	24.0794	0.588283	4767.55	60.1868	22.2043	0.0137066	5112.41	225.404	22.1284	0.0478697	7284.69	192.724	21.744	0.0287242	12603.6	338.831	21.1488	0.0291885	--	--	--	--	1.0
609780077793644613	52.985939853056564	-28.69679674577051	--	39538.8	6031.78	19.9074	0.165633	31775.8	409.134	20.1448	0.0139795	26683.0	1090.82	20.3344	0.0443856	27020.4	1214.81	20.3208	0.0488135	50730.3	2369.3	19.6368	0.050708	--	--	--	--	1.0
609780077793645026	52.972552675021674	-28.682123721400128	1.0	31473.7	793.915	20.1551	0.0273874	97788.1	126.326	18.9243	0.00140259	170707.0	480.003	18.3194	0.00305292	229063.0	470.132	18.0001	0.00222838	261805.0	667.193	17.8551	0.00276693	307996.0	3182.95	17.6786	0.0112204	1.0
609780077793645262	52.986892074360064	-28.673412724028246	1.0	2659.0	449.749	22.8382	0.183644	6934.43	46.943	21.7975	0.00734995	11716.3	165.379	21.228	0.0153254	14161.2	188.708	21.0223	0.0144682	15880.3	283.233	20.8979	0.0193646	19922.1	1723.72	20.6517	0.093941	1.0
609780077793645282	52.97437078645207	-28.672854232930916	1.0	1860.02	442.091	23.2262	0.258059	2466.5	49.1083	22.9198	0.0216171	4597.71	202.407	22.2436	0.0477979	7121.18	188.88	21.7686	0.0287978	8733.35	274.911	21.547	0.0341772	7645.4	1958.94	21.6915	0.278193	1.0
609780764988408641	52.98128653098051	-28.640623249005124	1.0	266.497	555.361	25.3358	2.2626	849.475	44.1727	24.0771	0.0564583	2814.36	124.134	22.7766	0.0478889	5978.17	100.207	21.9586	0.0181992	9091.2	203.75	21.5034	0.0243333	12425.0	1291.36	21.1643	0.112843	1.0
609780764988408857	52.97082685787927	-28.63114343953832	1.0	1465.22	671.764	23.4852	0.49778	2667.87	48.0124	22.8346	0.0195395	9084.21	132.068	21.5043	0.0157847	20945.5	130.077	20.5973	0.00674269	28051.8	253.622	20.2801	0.00981634	30360.7	1543.53	20.1942	0.0551985	1.0
609780764988408924	52.97280437648549	-28.628976583906766	1.0	-29.7438	568.233	--	--	816.562	34.9825	24.12	0.0465143	2664.67	93.0737	22.8359	0.0379234	5605.54	87.3585	22.0285	0.0169205	9384.87	204.029	21.4689	0.0236041	12297.9	1491.32	21.1754	0.131663	1.0
609780764988408986	52.95220290778823	-28.62600524625896	1.0	1609.53	727.228	23.3833	0.490565	1635.09	44.4143	23.3661	0.029492	4307.37	114.032	22.3145	0.0287435	7838.02	110.947	21.6645	0.0153685	9826.83	211.323	21.419	0.0233484	13425.1	1495.71	21.0802	0.120964	1.0
609780764988409299	52.98453257487234	-28.612603957304692	1.0	2753.56	509.272	22.8003	0.200807	3650.03	29.7631	22.4943	0.00885331	6216.85	76.1334	21.9161	0.0132962	6954.96	67.5167	21.7943	0.01054	9311.59	146.163	21.4774	0.0170427	9348.62	1241.41	21.4731	0.144176	1.0
609780764988409376	52.93172362649338	-28.610100726547877	1.0	3499.09	618.087	22.5401	0.191787	2942.46	44.2835	22.7282	0.0163402	3412.32	138.241	22.5674	0.0439857	4903.45	98.9794	22.1737	0.0219163	6114.97	222.904	21.934	0.0395774	7760.63	1296.2	21.6753	0.181343	1.0
609780764988409474	52.93697040322151	-28.60597284755832	1.0	3719.34	481.412	22.4738	0.140532	3508.42	31.5444	22.5372	0.00976191	5122.92	75.4172	22.1262	0.0159837	7559.09	70.5	21.7038	0.0101261	9980.12	166.807	21.4022	0.018147	8142.79	991.899	21.6231	0.132257	1.0
609780764988409708	52.942149530980316	-28.596835723591248	1.0	728.088	493.315	24.2445	0.73564	1264.75	34.3896	23.645	0.0295221	1432.1	65.3425	23.5101	0.049539	2608.26	84.747	22.8591	0.0352775	4624.55	163.971	22.2373	0.0384965	3409.78	1036.21	22.5682	0.329949	1.0
609780764988409790	52.93778614378815	-28.593488845270418	1.0	574.531	499.466	24.5017	0.943881	1064.63	31.6428	23.832	0.0322702	1512.36	66.7623	23.4509	0.0479293	2528.03	69.8977	22.893	0.0300196	3016.45	135.584	22.7013	0.0488021	2560.39	1205.44	22.8792	0.511167	1.0
609780764988409910	52.93402709778855	-28.58903688661977	1.0	2681.73	402.329	22.829	0.162889	2009.54	31.6292	23.1423	0.017089	3327.35	74.3417	22.5948	0.0242582	4282.19	72.9947	22.3208	0.0185076	5677.92	133.14	22.0145	0.0254591	6202.29	1257.01	21.9186	0.220044	1.0
609780764988409954	52.96627354323881	-28.58683645335011	1.0	989.835	338.122	23.9111	0.370882	1314.16	23.7982	23.6034	0.0196617	2553.17	54.8921	22.8823	0.0233429	3014.9	49.308	22.7018	0.017757	3212.02	115.556	22.6331	0.0390604	3091.53	921.718	22.6746	0.323705	1.0
609780764988409959	52.927485529152925	-28.58722235254101	1.0	-786.615	680.869	--	--	1635.9	39.3256	23.3656	0.0261002	2684.68	86.5621	22.8278	0.0350074	3687.74	111.373	22.4831	0.0327904	3884.96	201.033	22.4265	0.0561831	3791.1	1416.82	22.4531	0.405763	1.0
609780764988410214	52.92265370900696	-28.577915798953477	1.0	-358.72	221.132	--	--	554.734	24.7296	24.5398	0.0484012	1497.19	42.6026	23.4618	0.0308947	3841.19	54.5275	22.4388	0.0154126	5411.75	118.544	22.0667	0.0237831	6686.34	843.027	21.837	0.136892	1.0
609780764988410329	52.93702721258248	-28.57339422713925	1.0	596.245	167.165	24.4614	0.304401	642.408	25.1305	24.3805	0.0424731	943.172	45.6748	23.9635	0.0525787	1282.16	49.6353	23.6301	0.0420314	2006.34	123.696	23.144	0.066938	-482.103	894.35	--	--	1.0
609780764988410433	52.99051983137868	-28.569458200303867	1.0	1687.05	310.908	23.3322	0.200092	1425.8	43.3485	23.5149	0.0330095	1898.46	54.8521	23.204	0.0313702	3707.29	80.618	22.4774	0.0236102	3882.35	155.534	22.4273	0.0434967	2604.5	934.612	22.8607	0.389612	1.0
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
611257134226611400	53.10368543331558	-27.498471044536128	1.0	470.728	398.662	24.7181	0.919514	3349.04	67.4596	22.5877	0.0218699	13116.5	237.552	21.1055	0.0196636	19716.0	236.894	20.663	0.0130455	26695.7	284.672	20.3339	0.0115779	--	--	--	--	1.0
611257134226611431	53.15061588556635	-27.497303286465264	1.0	51.2462	397.528	27.1258	8.42228	2042.56	53.1588	23.1246	0.028257	4094.1	163.731	22.3696	0.0434206	6809.54	121.407	21.8172	0.0193575	8269.11	249.158	21.6064	0.0327146	--	--	--	--	1.0
611257134226611457	53.145312504650946	-27.496350536107865	1.0	3002.79	385.657	22.7062	0.139444	2442.62	48.5875	22.9304	0.021597	4470.21	154.952	22.2742	0.0376351	7775.3	116.206	21.6732	0.016227	9087.2	271.45	21.5039	0.0324328	--	--	--	--	1.0
611257134226611503	53.11245520260492	-27.494147627013835	1.0	100188.0	2202.27	18.898	0.023866	554607.0	569.222	17.04	0.00111435	1482770.0	1850.33	15.9723	0.00135488	2239130.0	1885.96	15.5248	0.000914488	2916520.0	1736.2	15.2378	0.000646336	--	--	--	--	1.0
611257134226611504	53.11698399795581	-27.49756534439678	1.0	3785.35	567.295	22.4547	0.162715	14924.7	97.5663	20.9652	0.00709773	34165.9	286.178	20.066	0.00909426	47867.8	301.012	19.6999	0.00682754	60091.6	344.054	19.453	0.00621637	--	--	--	--	1.0
611257134226611505	53.11110195438482	-27.492089478120832	1.0	953.199	594.453	23.952	0.677109	6049.73	115.426	21.9457	0.0207153	17948.1	298.027	20.765	0.0180286	40388.0	335.162	19.8844	0.00901004	57247.9	379.864	19.5056	0.00720432	--	--	--	--	1.0
611257134226611506	53.119532689130644	-27.49339209301625	1.0	269.911	355.708	25.3219	1.43086	3122.45	59.7073	22.6638	0.0207614	7730.43	174.388	21.6795	0.0244928	11385.5	181.352	21.2591	0.0172939	14840.2	253.721	20.9714	0.0185626	--	--	--	--	1.0
611257134226611507	53.10897681651336	-27.49129789967202	1.0	-710.264	323.819	--	--	1332.92	59.9009	23.588	0.0487925	3494.36	144.513	22.5416	0.0449019	9620.46	165.07	21.442	0.0186293	16480.0	197.354	20.8576	0.0130021	--	--	--	--	1.0
611257134226611510	53.11722736427128	-27.500433916342732	1.0	3423.3	827.527	22.5639	0.262459	4481.76	118.481	22.2714	0.0287029	7617.96	367.15	21.6954	0.0523273	8990.96	303.368	21.5155	0.0366344	8589.67	388.843	21.5651	0.0491499	--	--	--	--	1.0
611257134226611513	53.113803248234966	-27.494719701814915	1.0	2634.83	352.61	22.8481	0.1453	1015.24	57.636	23.8836	0.0616382	4233.7	161.912	22.3332	0.0415226	6759.51	174.401	21.8252	0.0280128	9859.19	200.786	21.4154	0.0221114	--	--	--	--	1.0
611257134226611515	53.11068108784129	-27.492868419028262	--	466.85	1338.82	24.7271	3.11365	10158.5	254.92	21.3829	0.0272457	20093.3	554.696	20.6424	0.0299728	35368.1	693.97	20.0285	0.0213036	36360.2	621.415	19.9984	0.0185558	--	--	--	--	1.0
611257134226611516	53.113295632661284	-27.493035554319103	1.0	1576.51	359.31	23.4058	0.247455	918.989	65.1544	23.9917	0.0769765	3212.89	160.582	22.6328	0.0542658	4899.15	173.755	22.1747	0.0385072	7329.94	204.394	21.7372	0.0302755	--	--	--	--	1.0
611257134226611518	53.10672167429923	-27.493985144292544	--	11625.7	8914.71	21.2365	0.832552	36487.4	1552.5	19.9946	0.0461968	58875.8	2917.63	19.4752	0.0538043	111904.0	4083.79	18.7779	0.0396225	156079.0	5485.51	18.4166	0.038159	--	--	--	--	1.0
611257134226611525	53.10936334993844	-27.492954031975184	--	4365.21	2851.31	22.3	0.709191	11264.8	401.841	21.2707	0.0387307	23832.1	831.452	20.4571	0.0378791	34764.2	1243.45	20.0472	0.0388348	36891.7	860.777	19.9827	0.0253329	--	--	--	--	1.0
611257134226611526	53.11296832534903	-27.492289651616343	--	6981.83	1737.21	21.7901	0.270152	7089.36	164.556	21.7735	0.0252018	15919.7	418.362	20.8952	0.0285326	21005.8	464.392	20.5942	0.0240033	32000.9	577.703	20.1371	0.0196005	--	--	--	--	1.0
611257134226611666	53.09592963830107	-27.490889464378586	1.0	5824.59	573.212	21.9868	0.10685	18229.3	131.409	20.7481	0.00782667	37692.2	288.151	19.9594	0.00830027	54771.0	305.786	19.5536	0.00606166	68918.3	465.231	19.3042	0.00732924	--	--	--	--	1.0
611257134226611701	53.10242838039273	-27.486976286763714	1.0	19316.3	661.687	20.6852	0.0371923	106484.0	217.811	18.8318	0.00222086	272736.0	595.379	17.8106	0.00237015	406893.0	605.147	17.3763	0.00161475	524674.0	756.458	17.1003	0.00156538	--	--	--	--	1.0
611257134226611846	53.09683357709753	-27.485504113414045	1.0	-122.658	343.009	--	--	1811.46	78.3514	23.2549	0.0469616	3606.62	159.18	22.5072	0.0479194	5203.85	174.124	22.1092	0.0363293	6668.54	289.43	21.8399	0.0471235	--	--	--	--	1.0
611257271665557609	52.772696201702765	-27.514444124290666	0.0	--	--	--	--	19339.2	96.0317	20.6839	0.00539139	60162.8	1760.88	19.4517	0.0317779	156753.0	444.269	18.412	0.0030772	247867.0	846.264	17.9145	0.00370691	301990.0	1730.72	17.7	0.0062224	0.0
611257271665559218	52.82742692885282	-27.518964788444382	1.0	--	--	--	--	7723.6	122.995	21.6805	0.0172899	13853.0	347.726	21.0461	0.0272532	20425.0	274.674	20.6246	0.0146009	--	--	--	--	29858.4	1887.73	20.2123	0.0686433	1.0
611257271665559232	52.818324739149794	-27.517739172656288	1.0	--	--	--	--	39332.4	149.843	19.9131	0.00413629	82729.0	469.273	19.1059	0.00615875	121532.0	386.399	18.6883	0.003452	--	--	--	--	174969.0	2392.65	18.2926	0.0148471	1.0
611257271665559258	52.813808783243864	-27.51163242507894	1.0	--	--	--	--	28732.0	147.199	20.2541	0.00556241	69100.8	401.565	19.3013	0.00630953	104170.0	498.945	18.8556	0.00520035	--	--	--	--	166596.0	2114.9	18.3458	0.0137832	1.0

The color selection requires a robust measurement of the g-r color, but g-band dropouts are not detected, so color cannot be robustly recovered from the flux measurement. Instead, get a lower limit on how red it is using the g-band flux error (i.e. use the upper limit to the g-band flux to measure the g-r color). To do this, define a new parameter called gmag to store the g-band magnitude error, for use in measuring colors. In this case, if sources have signal to noise less than 1 in g_band, gmag will instead hold the gmagErr from the object table.

In [48]:

gmag = np.zeros(len(results['gmag']), dtype='float')
S2N = np.asarray(results['g_cModelFlux']/results['g_cModelFluxErr'], dtype='float')

for i in range(len(gmag)):
    if S2N[i] < 1:
        gmag[i] = results['gmagErr'][i]
    else:
        gmag[i] = results['gmag'][i]

/tmp/ipykernel_62120/914783169.py:8: UserWarning: Warning: converting a masked element to nan.
  gmag[i] = results['gmag'][i]
/tmp/ipykernel_62120/914783169.py:6: UserWarning: Warning: converting a masked element to nan.
  gmag[i] = results['gmagErr'][i]

The cell below will identify Lyman-break galaxies (LBG) using the colors defined by Overzier et al. 2006.

In [49]:

GRlim = 1.5

whlbg = np.where(((gmag-results['rmag']) > 1.5)
                 & ((results['gmag'] - results['rmag'])
                    > (results['rmag'] - results['zmag'] + 1.1))
                 & ((results['rmag'] - results['zmag']) < 1.0))[0]

Plot the g-r vs r-z color color diagram. The plot shows the colors of the galaxies in the object table (blue) and those with colors consistent with z~4 galaxies (green).

In [50]:

RmZ = np.arange(0.4, 1.1, 0.1)

plt.plot(RmZ, RmZ + 1.1, '--', color='k')
plt.plot([1., 1.], [2.1, 3], '--', color='k')
plt.plot([-1, 0.4], [1.5, 1.5], '--', color='k')

RminusZ = results['rmag'] - results['zmag']
GminusR = gmag - results['rmag']

plt.plot(RminusZ, GminusR, '.', label='all galaxies', alpha=.1)
plt.plot(RminusZ[whlbg], GminusR[whlbg], 's', alpha=.5,
         label='color-selected z~4 candidates')

plt.xlabel('r - z')
plt.ylabel('g - r')
plt.xlim([-0.5, 3])
plt.ylim([-0.5, 3])
plt.legend()
plt.show()

Figure 9: g - r vs r - z color color diagram for galaxies with S/N > 20 in the ECDFS field (blue). Galaxies with blue r-z colors and red g-r colors consistent with with a lyman-break at redshift z=4 are in green.