---
title: ASTR 597 W23 HW 3
author: Andy Tzanidakis and friends
tags: LSST Class
---

[TOC]

# Homework 3 - ASTR 597 W23
- Author: Anastasios (Andy) Tzanidakis
- Course ID: ASTR 597 W23

# Problem 1

## Problem 1a - DP0.2 Tracts
How many DP0.2 tracts have r-band `'deepCoadd'`s?

We can query Butler and count the number of unique tracts for LSST-r in `deepCoadd`:

butler = dafButler.Butler('dp02', collections='2.2i/runs/DP0.2') # call butler
registry = butler.registry # fetch registry
datasetRefs = registry.queryDatasets(datasetType='deepCoadd', band='r') # basic butler query 

tracts = [j.dataId['tract'] for j in datasetRefs]

print (f"Number of unique tracts in deepCoadd LSST-r: {len(set(tracts))}")

>> Number of unique tracts in deepCoadd LSST-r: 157

There are 157 unique LSST-r band tracts in deepCoadd.

## Problem 1b - deepCoadd Plots
Find the tract and patch corresponding to (ra, dec) 62.0, -37.0. If there is more than one r-band `deepCoad` at that position, make a plot that explains why.

# target coordinates
ra, dec = 62.0, -37.0
target_band = 'r'

level = 10 # the resolution of the HTM grid
pixelization = lsst.sphgeom.HtmPixelization(level)

htm_id = pixelization.index(
    lsst.sphgeom.UnitVector3d(
        lsst.sphgeom.LonLat.fromDegrees(ra, dec)
    )
)

datasetRefs = registry.queryDatasets(datasetType='deepCoadd', htm20=htm_id, band=target_band)

Now we can query and append the unique `tract` and `patch` for each ID. For this approach I used the corners of the image using the WCS to determine if the region of interest was within this box:

_tract, _patch = [], []
for j in tqdm(datasetRefs):
    # Query deepCoadd box & WCS
    bbox = butler.get('deepCoadd.bbox', dataId=j.dataId, band=target_band)
    wcs = butler.get('deepCoadd.wcs', dataId=j.dataId, band=target_band)
    
    corners = [wcs.pixelToSky(bbox.beginX, bbox.beginY),
               wcs.pixelToSky(bbox.beginX, bbox.endY),
               wcs.pixelToSky(bbox.endX, bbox.endY),
               wcs.pixelToSky(bbox.endX, bbox.beginY)]
    
    # extrapolate coordinate box
    ra_list, dec_list = [], []
    for i in range(4):
        ra_list.append(round(corners[i].getRa().asDegrees(), 3))
        dec_list.append(round(corners[i].getDec().asDegrees(), 3))
        
    # Seek condition: is the target coordinate within the box?
    condition = (min(ra_list)<ra) & (max(ra_list)>ra) & (max(dec_list)>dec) & (min(dec_list)<dec)
    
    if condition:
        tract, patch = j.dataId['tract'], j.dataId['patch']
        
        _tract.append(tract)
        _patch.append(patch)

From our above code we have shown that there's more than r-band `deepCoadd` at this position. Here are the overlapping tract and patch list:
- Tract ID number: 3831
- Patch ID number: 10, 3

We can also visually show this by plotting the edges of the CCD and visualizing the central target coordinate we're querying:

|<img src="https://i.imgur.com/Ddw0KOh.png" alt="Image" width="450" height="400" style="display: block; margin: 0 auto" />|
|:--:|
| Figure 1. -- Overlapping CCD area of patch 10 and 3 for patch number 3831. |

## Problem 1c - PSF FWHM Distribution
Using the (tract, patch) with the lowest patch number, plot histograms of the seeing values of the images that were coadded in the r-band `deepCoadd` and `goodSeeingCoadd` for that (tract, patch). Tutorial notebook 09a (Custom Coadd)will give you some ideas on how to get the `coaddInputs`. The `coaddInputs.ccds` table will be helpful.

For each input `calexp`, you will need to retrieve the PSF and then use `psf.computeShape(psf.getAveragePosition()).getDeterminantRadius()` to get the PSF width in sigma. Convert these to FWHM in arcseconds for your plots.

my_username = getpass.getuser()

my_collection_identifier = 'custom_coadd_1'

my_outputCollection = 'u/'+my_username+'/'+my_collection_identifier
print('Name of new butler collection for my output: ', my_outputCollection)

simpleButler = SimplePipelineExecutor.prep_butler('dp02',
                                                  inputs=['2.2i/runs/DP0.2'],
                                                  output=my_outputCollection)

Now, we can query for each coadd input the PSF FHWM in arcseconds. We repeat this process for both `goodSeeingCoadd` and `deepCoadd`:

doa = simpleButler.get("deepCoadd.coaddInputs", tract=3831, patch=3, band='r')
ccd_table = doa.ccds.asAstropy() # contains visit ID 

seeing_values = []

datasetType = 'calexp'
j = 0 

for vis in tqdm(ccd_table['visit']):

    dataId = {'visit': vis, 'detector': 175}
    calexp = butler.get(datasetType, dataId=dataId, band='r')
    
    psf = calexp.getPsf()
    sigma = psf.computeShape(psf.getAveragePosition())
    seeing = sigma.getDeterminantRadius() * calexp.wcs.getPixelScale().asArcseconds() * 2.355
    
    seeing_values.append(seeing)

# FWHM 
seeing_values = np.array(seeing_values)

|<img src="https://i.imgur.com/V3nckUa.png" alt="Image" width="680" height="280" style="display: block; margin: 0 auto" />|
|:--:|
| Figure N. -- Distribution of the PSF FWHM sigma in arcseconds for `deepCoadd`(purple) and `goodSeeingCoadd` (pink).|

# Problem 2 - DIY Good-seeing Coadd

## Problem 2a
Complete the `09a_Custom_Coadd` tutorial notebook, but rather than coadding all of the visits within a temporal window, instead coadd the three first three visits that have seeing values in the range included in the `goodSeeingCoadd` inputs. Display the result and compare to the standard `deepCoadd` and `goodSeeingCoadd`.

Coadded template images are a necessary input to alert production. What are the implications of this exercise for alert generation early in the LSST survey? How many images were taken before three images had good enough seeing? How long did the survey run before the images were acquired?

We will begin as suggested from the`09_Custom_Coadd` tutorial:

# coordinates & filter
ra_targ = 55.745834
dec_targ = -32.269167
flt = 'i'

my_spherePoint = lsst.geom.SpherePoint(ra_targ*lsst.geom.degrees, dec_targ*lsst.geom.degrees)

skymap = butler.get('skyMap')
tract = skymap.findTract(my_spherePoint)
my_tract = tract.tract_id
my_patch = tract.findPatch(my_spherePoint).getSequentialIndex()
print('My tract and patch: ', my_tract, my_patch)

_dataId = {'band': flt, 'tract': my_tract, 'patch': my_patch}
good_seeing_coadd = butler.get('goodSeeingCoadd', dataId=_dataId)

>> My tract and patch:  4431, 17

We're interested in tract 4431 and patch 17. Now we can continue our query:

input_good_s = good_seeing_coadd.getInfo().getCoaddInputs()
visit_good_info = input_good_s.ccds.asAstropy()
visitTableRef = list(butler.registry.queryDatasets('visitTable'))
visitTable = butler.get(visitTableRef[0])
my_Goodcoadd_visits_mjds = visitTable.loc[input_good_s.visits['id']]['expMidptMJD'] #MJD times

We can also write a quick helper function to fetch the PSF FWHM in arcseconds like we did in the previous excercise:

def fetch_goodSeeingCoadd_vals(tract, patch, band, visit_id):
    """ Return the PSF FWHM in arcseconds for a given tract, patch, band and visitID. 
    
        Parameters
        ----------
        tract (int): tract number
        patch (int): patch number
        band (str): photometric bandpass
        visit_id (int): visit ID
    
        Returns
        -------
        PSF FHWM (float): Full-width half maximum of the PSF
    
    """
    query = simpleButler.get("goodSeeingCoadd.coaddInputs", tract=tract, patch=patch, band=band)
    ccd_table = query.ccds.asAstropy() # contains visit ID 
    
    dataId = {'visit': visit_id, 'detector': 175}
    calexp = butler.get('calexp', dataId=dataId, band=band)
    
    psf = calexp.getPsf()
    sigma = psf.computeShape(psf.getAveragePosition())
    seeing = sigma.getDeterminantRadius() * calexp.wcs.getPixelScale().asArcseconds() * 2.355
    
    return seeing

# Apply helper function and print the PSF FWHM in arcseconds including the ID's of the first 3  
vis_ids = visitTable.loc[input_good_s.visits['id']].visit.values
seeing_goodS_ca = [fetch_goodSeeingCoadd_vals(4431, 17, 'i', j) for j in tqdm(vis_ids)]

>> IDs of visit: [174602 177422 192350]
>> PSF FWHM [arcsec]: [0.7476371542388519, 0.7325203429002716, 0.7566516395201798]

As we see in our above code, we will be coadding visits [174602 177422 192350] that span a PSF FWHM 0.73 - 0.76 arcseconds. Now we can continue with the coadding sequence:

ids_good = tuple(vis_ids[0:3]) # ids of the best images

yaml_file = '$DRP_PIPE_DIR/pipelines/LSSTCam-imSim/DRP-DP0.2.yaml'
steps = 'makeWarp,assembleCoadd'
my_uri = yaml_file + '#' + steps

assembleCoaddPipeline = Pipeline.from_uri(my_uri)

assembleCoaddPipeline.addConfigOverride('makeWarp', 'doApplyFinalizedPsf', False)

queryString = f"tract = 4431 AND patch = 17 AND " + \
              f"visit in {ids_good} AND skymap = 'DC2'"

print(queryString)
>> tract = 4431 AND patch = 17 AND visit in (174602, 177422, 192350) AND skymap = 'DC2'

spe = SimplePipelineExecutor.from_pipeline(assembleCoaddPipeline, where=queryString, butler=simpleButler)
quanta = spe.run()
my_new_deepCoadd = simpleButler.get(quanta[-1].outputs['deepCoadd'][0])

my_new_deepCoadd_inputs = simpleButler.get("deepCoadd.coaddInputs",
                                           {'band': 'i', 'tract': 4431, 'patch': 17})

|<img src="https://i.imgur.com/8HQJJXM.jpg" alt="Image" width="950" height="280" style="display: block; margin: 0 auto" />|
|:--:|
| Figure N. -- `goodSeeingCoadd` three image coadded image (left), `deepCoadd` coadded image using the full coadded visit history >500 images stacked. |

- We find that in our custom `goodSeeingCoadd` using the first three `calexp` images that the depth not similar to that of the `deepCoadd` that uses the full stack of images. The DIY `goodSeeingCoadd` coadded image spanned approximately 27 days.

- Deep coadded template images will be vital for the alert production of LSST. Coadded images play a vital role in both the photometric depth and decreasing the point spread function width ([Zackay & Ofek 2015](https://arxiv.org/pdf/1512.06872.pdf)). We also discover from this exercise that it will take a while (i.e at least >27 days) to build good enough `deepCoaddd` templates that can be used for the LSST alert production.

# Acknowledgements
Many concepts and ideas for the workflow of this project were heavily inspired by conversations I had with Tom Wagg, David Wang, Tobin Wainer, and Jake Kurlander.