<jupyter_start><jupyter_text>DC2 Object Catalog Run2.2i GCR tutorial -- Part III: Guided ChallengesOwners: **<NAME> [@fjaviersanchez](https://github.com/LSSTDESC/DC2-analysis/issues/new?body=@fjaviersanchez), <NAME> [@EiffL](https://github.com/LSSTDESC/DC2-analysis/issues/new?body=@EiffL)**  Last Run: **2020-07-02** (by @jrbogart)This notebook is the third of four in the Run2.2i object catalog with GCR series ([Part I](object_gcr_1_intro.ipynb), [Part II](object_gcr_2_lensing_cuts.ipynb)). Here, we propose some challenges for you as science use cases of the object catalogs. We will provide a solution here but you are encouraged to create your own!__Logistics__: This notebook is intended to be run through the JupyterHub NERSC interface available here: https://jupyter-dev.nersc.gov. To setup your NERSC environment, please follow the instructions available here: https://confluence.slac.stanford.edu/display/LSSTDESC/Using+Jupyter-dev+at+NERSC<jupyter_code>%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Challenge 1: Galaxy counts-in-cells Blending affects the accuracy of centroid and flux measurements. It can potentially generate a systematic effect in different measurements (for example 2-point statistics). The stack returns a very useful value to check (partially) for the presence of these kind of systematics, which is the `blendedness` parameter (more details on Section 4.9.11 of [Bosch et al. 2017](https://arxiv.org/pdf/1705.06766.pdf)* Q: Why only partially? A very simple tool to measure the different statistical moments of galaxies is Counts-in-cells (CiC) [Peebles et al. 1980](https://press.princeton.edu/titles/724.html). Here, we are going to use a simplified version of CiC to check the possible systematic effects due to differences in the `blendedness` measurements. So what's CiC?1) Count the number of galaxies in a cell of a given scale.2) Measure the density contrast distribution and its moments.3) Change the scale and repeat. Note: We are focusing on the galaxy number CiC but it is possible to do with shapes<jupyter_code>from utils.cic import cic_analysis
import healpy as hp
import GCRCatalogs
from GCR import GCRQuery
# Load the object catalog
run = 'run2.2i_dr3'
tract = 3081    # has good depth
tract_s = str(tract)
catalog = GCRCatalogs.load_catalog(f'dc2_object_{run}')
#catalog = GCRCatalogs.load_catalog('dc2_object_run1.2i_all_columns')
# Let's use almost the same cuts as in the WL sample

cic_cuts_base = [
    GCRQuery((np.isfinite, 'mag_i_cModel')), # (from this and below) remove nan entries
    GCRQuery((np.isfinite, 'ext_shapeHSM_HsmShapeRegauss_resolution')),
    GCRQuery((np.isfinite, 'ext_shapeHSM_HsmShapeRegauss_e1')),
    GCRQuery((np.isfinite, 'ext_shapeHSM_HsmShapeRegauss_e2')),
    GCRQuery('good'), 
    GCRQuery('snr_i_cModel >= 10'), # (from this and below) cut on object properties
    GCRQuery('ext_shapeHSM_HsmShapeRegauss_sigma <= 0.4'),
]

is_blended = GCRQuery('blendedness >= 10**(-0.375)')

cic_cuts_nb = cic_cuts_base + [~is_blended]
cic_cuts_b = cic_cuts_base + [is_blended]

quantities = ['ra','dec']

# was tract 4849.  3081 has greater depth for Run2.2i_dr3
d_nb = catalog.get_quantities(quantities, 
                              filters=cic_cuts_nb, 
                              native_filters=[f'tract == {tract_s}'])
d_b = catalog.get_quantities(quantities, 
                             filters=cic_cuts_b, 
                             native_filters=[f'tract == {tract_s}'])
# Read in the depth map created in the previous tutorial 
import pickle
import gzip
f = gzip.open(f'assets/m10map_{run}_{tract_s}.pklz', 'rb')
m10map = pickle.load(f)
f.close()
mask = np.zeros_like(m10map)
mask[m10map > 23.0] = 1.0
hp.gnomview(mask, rot=(d_nb['ra'].mean(), d_nb['dec'].mean()), title='Run 2.2i Depth', reso=0.5, unit='10-$\sigma$ i-band depth')
sigma_b, sigma_err_b, skw_b, skw_err_b, kurtosis_b, kurtosis_err_b, pixel_scale = cic_analysis(d_b, mask, nboot=100)
sigma_nb, sigma_err_nb, skw_nb, skw_err_nb, kurtosis_nb, kurtosis_err_nb, _ = cic_analysis(d_nb, mask, nboot=100)
f, ax = plt.subplots(1, 1, figsize=(6,6))
ax.errorbar(pixel_scale[2:], sigma_b[2:], sigma_err_b[2:], fmt='o', linestyle='none', label='High blendedness')
ax.errorbar(pixel_scale[2:], sigma_nb[2:], sigma_err_nb[2:], fmt='x', linestyle='none', label='Low blendedness')
ax.legend()
ax.set_xlabel('Pixel scale [deg]')
ax.set_ylabel(r'$\sigma$')
f.draw(f.canvas.get_renderer())<jupyter_output><empty_output><jupyter_text>There's definitely something going on with the high blendedness sources!  Challenge 2: Check if PSF residuals are within requirementsIn this section, we will try to apply all the tools we have covered during this tutorial to test the quality of the DM stack PSF model on run 2.2i_dr3.The challenge will be to select a clean sample of stars, compute their size and ellipticity using second moments, and compare those to the PSF model predicted by the DM stack. We will test the one point and two point fuctions of these residuals to make diagnostic plots that would directly go into a weak lensing shape catalog paper.See for example Section 3. of [(Zuntz et al. 2017)](https://arxiv.org/pdf/1708.01533.pdf), (and in particular section 3.2) for an introduction to, and real life example of, PSF systematics tests.  Step 1: Select a clean sample of point sourcesWe want to restrict our sample to objects that follow these constraints: - Point sources - Not corrupted or with any defects - Successful second moment measurements - Sufficiently high signal to noise in the i-band, above 50 Remember to use [SCHEMA.md](https://github.com/LSSTDESC/gcr-catalogs/blob/master/GCRCatalogs/SCHEMA.mdschema-for-dc2-coadd-catalogs) as a reference to build your cuts<jupyter_code>import GCRCatalogs
from GCR import GCRQuery

catalog = GCRCatalogs.load_catalog(f'dc2_object_{run}')

filters = [
    # Add the required filters
    # ...
]<jupyter_output><empty_output><jupyter_text>Step 2: Computes size and ellipticity from second momentsWe will use the following definitions:  $g_1 = \frac{I_{xx} - I_{yy}}{I_{xx} + I_{yy}}$  $g_2 = \frac{2 I_{xy}}{I_{xx} + I_{yy}}$  $\sigma = ( I_{xx} I_{yy} - I_{xy}^2)^{1/4}$ Using the `add_derived_quantity` of the GCR (documented [here](https://yymao.github.io/generic-catalog-reader/index.htmlGCR.BaseGenericCatalog.add_derived_quantity)), add modifiers to compute these quantities for the sources and the PSF model evaluated at the position of the sources.Again the schema is your friend ;-)<jupyter_code># Here is a template to fill for g1
g1_modif = lambda ixx,iyy,ixy: (ixx-iyy)/(ixx+iyy)
catalog.add_derived_quantity('g1', g1_modif, 'Ixx', 'Iyy', 'Ixy')

# Define in the same way g2, sigma, psf_g1, psf_g2 psf_sigma
#...<jupyter_output><empty_output><jupyter_text>Step 3: Extract sampleNow that we have all the pieces, let's extract the quantities specified below from the catalog:<jupyter_code>quantities = ['ra', 'dec', 
              'mag_i', 'snr_i_cModel', 'psf_fwhm_i',
              'g1', 'g2', 'sigma',
              'psf_g1', 'psf_g2', 'psf_sigma']

# Extract these quantities from the catalog into data
# data = ...<jupyter_output><empty_output><jupyter_text>Step 4: Size and ellipticity residuals as a function of magnitude and seeingTry to reproduce plots similar to the ones you can find in the PSF section of your favorite experiment's shape catalog paper (for instance Section 4. of [Mandelbaum et al. 2017](https://arxiv.org/pdf/1705.06745.pdf)).For instance you can look at the fractional difference in size, as a function of magnitude, or seeing. You can also look at the distribution of ellipticity residuals, make sure they are centered on 0 , and again see if you can spot a dependence on seeing or magnitude.<jupyter_code># plt.figure(figsize=(10,5))

# plt.subplot(121)
# plt.hist2d(data['mag_i'], (data['sigma'] - data['psf_sigma'])/data['psf_sigma'], 100, range=[[15,23],[-0.02,0.02]]);
# plt.xlabel('i mag')
# plt.ylabel('$f \delta_\sigma$')
# plt.colorbar(label='Number of objects')
# plt.subplot(122)
# plt.hist2d(data['psf_fwhm_i'], (data['sigma'] - data['psf_sigma'])/data['psf_sigma'], 100, range=[[0.4,1.0],[-0.02,0.02]]);
# plt.xlabel('seeing FWHM (arcsec)')
# plt.colorbar(label='Number of objects');
# plt.figure(figsize=(15,10))

# plt.subplot(221)
# plt.hist2d(data['mag_i'], (data['g1'] - data['psf_g1']), 100, range=[[15,23],[-0.02,0.02]]);
# plt.xlabel('i mag')
# plt.ylabel('$g_1 - g_1^{PSF}$')
# plt.colorbar(label='Number of objects')
# plt.subplot(222)
# plt.hist2d(data['psf_fwhm_i'], (data['g1'] - data['psf_g1']), 100, range=[[0.4,1.0],[-0.02,0.02]]);
# plt.xlabel('seeing FWHM (arcsec)')
# plt.ylabel('$g_1 - g_1^{PSF}$')
# plt.colorbar(label='Number of objects')
# plt.subplot(223)
# plt.hist((data['g1'] - data['psf_g1']), 100, range=[-0.04,0.04]);
# plt.xlabel('$g_1 - g_1^{PSF}$')
# plt.axvline(0)
# plt.subplot(224)
# plt.hist((data['g2'] - data['psf_g2']), 100, range=[-0.04,0.04]);
# plt.xlabel('$g_2 - g_2^{PSF}$')
# plt.axvline(0)<jupyter_output><empty_output><jupyter_text>Step 5: Compute $\rho$-statistics in StileNo shear catalog paper would be complete without the so-called $\rho$-statistics ([Rowe 2010](https://academic.oup.com/mnras/article/404/1/350/3101560), [Jarvis et al. 2016](https://academic.oup.com/mnras/article-abstract/460/2/2245/2609178?redirectedFrom=fulltext)), which check the two-point correlations of the PSF residuals. See section 3.4 of Jarvis et al. 2016 for instance to see the definition of these statistics.We are going to use the [Stile](https://github.com/msimet/Stile) package developed by <NAME> [@msimet](https://github.com/msimet), <NAME> [@HironaoMiyatake](https://github.com/HironaoMiyatake), <NAME> [@rmandelb](https://github.com/rmandelb), and <NAME> [@dr-guangtou](https://github.com/dr-guangtou).This is an incredibly useful package which already implements a range of WL related systematics tests, including the $\rho$ statistics. Checkout the documentation for Stile [here](http://stile.readthedocs.io/en/latest/index.html).Uncomment the following cells to try them out.<jupyter_code># import pandas
# import stile

# # Stile expects numpy structured arrays, here is an easy way to do that:
# d = pandas.DataFrame(data)
# # We also add a weight column, set to 1
# d['w'] =1
# d = d.to_records(index=False)
# # Here is an example of how to compute the rho1 statistics
# stile_args = {'ra_units': 'degrees', 'dec_units': 'degrees',
#               'min_sep': 0.05, 'max_sep': 1, 'sep_units': 'degrees', 'nbins': 20}

# rho1 = stile.CorrelationFunctionSysTest('Rho1')

# r1 = rho1(d, config=stile_args)
# # Have a look at the content of the result array
# r1
# # Stile offers utility functions to generate 
# f = rho1.plot(r1)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1312096029547553, 'fraction_numerical': 0.03250230840258541, 'mean_word_length': 5.044084821428571, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 38, 'https://': 17, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '208814', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3961, 'n_tokens_neox': 3715, 'n_words': 1091, 'path': 'tutorials/object_gcr_3_challenges.ipynb', 'repo_name': 'LSSTDESC/DC2-analysis', 'size': '10880'}