<jupyter_start><jupyter_text>Synthetic strong gravitational lens injection Demo This notebook shows how to insert synthetic lensed sources behind DC2 and Rubin images. Owner: **<NAME>** ([@herjy](https://github.com/herjy/DESC-Lamp))Last Verified to Run: **2021-11-22** (by @herjy)This early version of the demo aims at generating a lensed source galaxy image with an arbitrary lens profile and insert it in the background of an existing DC2 image of a galaxy at the coadd level. The content of this notebook builds up on the sythetic sourrce injection (example notebook){https://github.com/LSSTDESC/ssi-tools/blob/main/examples/DC2_calexp_injection.ipynb} by <NAME> inn the ssi tools packa LogisticsThis is intended to be runnable at NERSC through the https://jupyter.nersc.gov interface from a local git clone of https://github.com/herjy/DESC-Lamp in your NERSC directory.  But you can also run it wherever, with appropriate adjustment of the 'repo' location to point to a place where you have a Butler repo will all of the images. This notebook uses the `desc-stack-weekly-latest` kernel. Instructions for setting up the proper DESC python kernel can be found here: https://confluence.slac.stanford.edu/x/o5DVE<jupyter_code># A few common packages
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.pyplot as plt
from matplotlib.transforms import Affine2D
import pickle
import pandas as pd
%matplotlib inline

# And also DESC packages to get the data path
import GCRCatalogs
import galsim

from desclamp import train_set<jupyter_output><empty_output><jupyter_text>Setup Fetchinng Butler for deep coadd of run 2.2i dr6, wide field deep. and extracting cutout information for a given location<jupyter_code>butler = Butler(REPOS['2.2i_dr6_wfd'])

skymap = butler.get("deepCoadd_skyMap")

# Near the center of DC2
ra = 65.0  # degrees
dec = -35.0 # degrees
point = geom.SpherePoint(ra, dec, geom.degrees)
cutoutSize = geom.ExtentI(201, 201)

tractInfo = skymap.findTract(point)
patchInfo = tractInfo.findPatch(point)
xy = geom.PointI(tractInfo.getWcs().skyToPixel(point))
bbox = geom.BoxI(xy - cutoutSize//2, cutoutSize)
coaddId = {
    'tract':tractInfo.getId(), 
    'patch':f"%d,%d"%patchInfo.getIndex(),
    'filter':'r'
}
coadd = butler.get("deepCoadd_sub", bbox=bbox, immediate=True, dataId=coaddId)


# Get visits contributing to this cutout
visitIds = []
for ccd in tqdm(coadd.getInfo().getCoaddInputs().ccds):
    dataId = {
        'visit':ccd['visit'], 
        'detector':ccd['ccd'], 
        'filter':ccd['filter'],
        'tract':tractInfo.getId()
    }
    wcs = butler.get('calexp_wcs', dataId=dataId)   # Note 99% of cell execution time is this line and the next
    bbox = butler.get('calexp_bbox', dataId=dataId)
    poly = makeSkyPolygonFromBBox(bbox, wcs)
    if poly.contains(point.getVector()):
        visitIds.append(dataId)


plt.close('all')

# Get image cutouts
n_max_cutout = 4
 
#file=  open('lensed_source','rb')
profile = galsim.Sersic(n, half_light_radius=0.2, flux=1.0)

lensed_source = pickle.load(file, encoding='bytes')
gsobj = galsim.InterpolatedImage(galsim.Image(lensed_source), scale = 0.1, flux = 100000)

for dataId in visitIds[:n_max_cutout]:
    wcs = butler.get('calexp_wcs', dataId=dataId)
    bbox = butler.get('calexp_bbox', dataId=dataId)
    xy = geom.PointI(wcs.skyToPixel(point))
    cutout = geom.BoxI(xy - cutoutSize//2, cutoutSize)
    bbox.clip(cutout)
    image = butler.get("calexp_sub", bbox=bbox, immediate=True, dataId=dataId)

    mat = np.eye(3)
    mat[:2,:2] = wcs.getCdMatrix()
    transform = Affine2D(mat)
    
    #Image to inject
    point = SpherePoint(ra, dec, geom.degrees)
    
    fig = plt.figure()
    arr = np.copy(image.maskedImage.image.array)
    plot_extents = 0, bbox.width, 0, bbox.height
    helper = floating_axes.GridHelperCurveLinear(
        transform, plot_extents, 
        tick_formatter1=DictFormatter({}),
        tick_formatter2=DictFormatter({}),
        grid_locator1=MaxNLocator(nbins=1),
        grid_locator2=MaxNLocator(nbins=1),
    )
    ax = floating_axes.FloatingSubplot(fig, 111, grid_helper=helper)
    ax.imshow(arr, vmin=0, vmax=300, transform=transform+ax.transData, cmap = 'gist_stern')
    ax.set_title(repr(dataId))
    ax.scatter(
        xy.x - bbox.minX, 
        xy.y - bbox.minY, 
        c='r', marker='+', transform=transform+ax.transData
    )
    fig.add_subplot(ax)
    fig.show()
    
    print("Injected images")
    _add_fake_sources(image, [(point, gsobj)])
    
    fig2 = plt.figure()
    inj_arr = image.maskedImage.image.array
    plot_extents = 0, bbox.width, 0, bbox.height
    helper = floating_axes.GridHelperCurveLinear(
        transform, plot_extents, 
        tick_formatter1=DictFormatter({}),
        tick_formatter2=DictFormatter({}),
        grid_locator1=MaxNLocator(nbins=1),
        grid_locator2=MaxNLocator(nbins=1),
    )
    ax = floating_axes.FloatingSubplot(fig2, 111, grid_helper=helper)
    ax.imshow(inj_arr, vmin=0, vmax=300, transform=transform+ax.transData, cmap = 'gist_stern')
    ax.set_title(repr(dataId))
    ax.scatter(
        xy.x - bbox.minX, 
        xy.y - bbox.minY, 
        c='r', marker='+', transform=transform+ax.transData
    )
    fig2.add_subplot(ax)
    fig2.show()
    
    fig3 = plt.figure()
    ax = floating_axes.FloatingSubplot(fig3, 111, grid_helper=helper)
    ax.imshow(inj_arr-arr, vmin=0, transform=transform+ax.transData, cmap = 'gist_stern')
    ax.set_title(repr(dataId))
    ax.scatter(
        xy.x - bbox.minX, 
        xy.y - bbox.minY, 
        c='r', marker='+', transform=transform+ax.transData
    )
    fig3.add_subplot(ax)
    fig3.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11917735981719108, 'fraction_numerical': 0.01669889259975391, 'mean_word_length': 3.9178910976663786, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 11, 'https://': 5, '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': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '329085', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2099, 'n_tokens_neox': 2020, 'n_words': 463, 'path': 'notebooks/Lens_injection.ipynb', 'repo_name': 'herjy/DESCLamp', 'size': '5699'}