import numpy as np
import galsim
import matplotlib.pyplot as plt
from shrbk.data import get_galaxy_image
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
/tmp/ipykernel_4318/947104330.py in <module>
      2 import galsim
      3 import matplotlib.pyplot as plt
----> 4 from shrbk.data import get_galaxy_image

ModuleNotFoundError: No module named 'shrbk'

12.3. Blending bias

12.3.1. Definition

Blending bias is present when working with crowded fields. It is defined as the bias that arise due to the presence of other galaxies near the one which is being measured. It can have two origins : either the neighbour galaxies directly impact the measured shape, or if blended galaxies are removed from the sample it may give rise to a selection bias.

12.3.2. Shape measurement impact through examples

Here we will show the impact of unrecognized blends on shape measurement. We will reuse here the examples from the Blending bias section ; as a reminder, we used an exponential profile for galaxies, and applied shear to each of the galaxies. We also convolved the galaxies with a simple Moffat PSF profile.

12.3.2.1. Unblended case

First we generate a central galaxy, whose shape will be measured using the KSB method. It will be used as a reference to compare with the cases where there is actual blending.

image_central, image_psf = get_galaxy_image(gal_flux=1.e5,gal_r0=1.4,g1=0.1,g2=0.2,psf_beta=5,psf_re=1.0,
                                            pixel_scale=0.2)

plt.imshow(image_central.array)
plt.title("Central galaxy")
plt.show()

measured_shear = galsim.hsm.EstimateShear(image_central,image_psf,shear_est='KSB')
print("Shear measurement results : g1 = ",round(measured_shear.corrected_g1,3),", g2 = ",round(measured_shear.corrected_g2,3))
print("Relative error : \u03B4g1 = ",round((measured_shear.corrected_g1 - 0.1)/0.1 * 100,3),"%, \u03B4g2 = ",round((measured_shear.corrected_g2 - 0.2)/0.2*100,3),"%")
../_images/blending-bias_4_0.png
Shear measurement results : g1 =  0.103 , g2 =  0.206
Relative error : δg1 =  2.824 %, δg2 =  2.811 %

Then we generate a companion galaxy with which it will be blended.

image_companion, _ = get_galaxy_image(gal_flux=3.e4,gal_r0=0.6,g1=0.1,g2=0.2,psf_beta=5,psf_re=1.0,
                                            pixel_scale=0.2,shift_x=15,shift_y=15)

plt.imshow(image_companion.array)
plt.title("Companion galaxy")
plt.show()
../_images/blending-bias_6_0.png

We will first generate a scene where the two galaxies are far apart from each other. Observe that the shape measurement is unaffected by the presence of the companion galaxy.

blended_scene = galsim.ImageD(256,256)
blended_scene += image_central
blended_scene += image_companion

plt.imshow(blended_scene.array)
plt.title("Scene with separated galaxies")
plt.show()

measured_shear = galsim.hsm.EstimateShear(blended_scene,image_psf,shear_est='KSB')
print("Shear measurement results : g1 = ",round(measured_shear.corrected_g1,3),", g2 = ",round(measured_shear.corrected_g2,3))
print("Relative error : \u03B4g1 = ",round((measured_shear.corrected_g1 - 0.1)/0.1 * 100,3),"%, \u03B4g2 = ",round((measured_shear.corrected_g2 - 0.2)/0.2*100,3),"%")
../_images/blending-bias_8_0.png
Shear measurement results : g1 =  0.103 , g2 =  0.206
Relative error : δg1 =  2.824 %, δg2 =  2.811 %

If we plot the light profile of the image along the x axis, we can see that the two profiles are completely separated here.

fig, (ax1,ax2) = plt.subplots(1,2,figsize=(20,10))
profile_central = np.sum(image_central.array, axis=0)
profile_companion = np.sum(image_companion.array, axis=0)
profile_blended = np.sum(blended_scene.array, axis=0)
ax1.plot(list(range(256)),profile_blended,color="blue",label="Blended scene",linewidth=5,alpha=0.3)
ax1.plot(list(range(256)),profile_central,'--',color="red",label="Central galaxy",linewidth=2)
ax1.plot(list(range(256)),profile_companion,'--',color="green",label="Companion galaxy",linewidth=2)
ax1.set_title("x axis profile")
profile_central = np.sum(image_central.array, axis=1)
profile_companion = np.sum(image_companion.array, axis=1)
profile_blended = np.sum(blended_scene.array, axis=1)
ax2.plot(list(range(256)),profile_blended,color="blue",label="Blended scene",linewidth=5,alpha=0.3)
ax2.plot(list(range(256)),profile_central,'--',color="red",label="Central galaxy",linewidth=2)
ax2.plot(list(range(256)),profile_companion,'--',color="green",label="Companion galaxy",linewidth=2)
ax2.set_title("y axis profile")
plt.legend(fontsize="x-large")
plt.show()
../_images/blending-bias_10_0.png

12.3.2.2. Blended case

Now, we will generate a scene where the galaxies truly overlap each other. We regenerate the image of the companion galaxy closer to the central one.

image_companion, _ = get_galaxy_image(gal_flux=3.e4,gal_r0=0.6,g1=0.1,g2=0.2,psf_beta=5,psf_re=1.0,
                                            pixel_scale=0.2,shift_x=5,shift_y=6)

We can now carry the measurement on the blended image. Notice how the measurement is largely affected by the presence of the companion galaxy near the central one.

blended_scene = galsim.ImageD(256,256)
blended_scene += image_central
blended_scene += image_companion

plt.imshow(blended_scene.array)
plt.title("Blended scene")
plt.show()

measured_shear = galsim.hsm.EstimateShear(blended_scene,image_psf,shear_est='KSB')
print("Shear measurement results : g1 = ",round(measured_shear.corrected_g1,3),", g2 = ",round(measured_shear.corrected_g2,3))
print("Relative error : \u03B4g1 = ",round((measured_shear.corrected_g1 - 0.1)/0.1 * 100,3),"%, \u03B4g2 = ",round((measured_shear.corrected_g2 - 0.2)/0.2*100,3),"%")
../_images/blending-bias_14_0.png
Shear measurement results : g1 =  0.07 , g2 =  0.296
Relative error : δg1 =  -30.452 %, δg2 =  47.837 %

We can plot again the light profile and see how the profiles overlap : there are no longer two separate peaks but instead only one with two “subpeaks”.

fig, (ax1,ax2) = plt.subplots(1,2,figsize=(20,10))
profile_central = np.sum(image_central.array, axis=0)
profile_companion = np.sum(image_companion.array, axis=0)
profile_blended = np.sum(blended_scene.array, axis=0)
ax1.plot(list(range(256)),profile_blended,color="blue",label="Blended scene",linewidth=5,alpha=0.3)
ax1.plot(list(range(256)),profile_central,'--',color="red",label="Central galaxy",linewidth=2)
ax1.plot(list(range(256)),profile_companion,'--',color="green",label="Companion galaxy",linewidth=2)
ax1.set_title("x axis profile")
profile_central = np.sum(image_central.array, axis=1)
profile_companion = np.sum(image_companion.array, axis=1)
profile_blended = np.sum(blended_scene.array, axis=1)
ax2.plot(list(range(256)),profile_blended,color="blue",label="Blended scene",linewidth=5,alpha=0.3)
ax2.plot(list(range(256)),profile_central,'--',color="red",label="Central galaxy",linewidth=2)
ax2.plot(list(range(256)),profile_companion,'--',color="green",label="Companion galaxy",linewidth=2)
ax2.set_title("y axis profile")
plt.legend(fontsize="x-large")
plt.show()
../_images/blending-bias_16_0.png

12.3.2.3. Blended case with different shear

Bias due to blended galaxies is hard to estimate, as it not only depends on the presence of the neighbour but also on its shear, which may be different if the two galaxies are situated at different redshifts. We can see it by modifying the shear of the companion galaxy in the previous scene.

image_companion, _ = get_galaxy_image(gal_flux=3.e4,gal_r0=0.6,g1=-0.5,g2=-0.3,psf_beta=5,psf_re=1.0,
                                            pixel_scale=0.2,shift_x=5,shift_y=6)
blended_scene = galsim.ImageD(256,256)
blended_scene += image_central
blended_scene += image_companion

plt.imshow(blended_scene.array)
plt.title("Blended scene with two different shears")
plt.show()

measured_shear = galsim.hsm.EstimateShear(blended_scene,image_psf,shear_est='KSB')
print("Shear measurement results : g1 = ",round(measured_shear.corrected_g1,3),", g2 = ",round(measured_shear.corrected_g2,3))
print("Relative error : \u03B4g1 = ",round((measured_shear.corrected_g1 - 0.1)/0.1 * 100,3),"%, \u03B4g2 = ",round((measured_shear.corrected_g2 - 0.2)/0.2*100,3),"%")
../_images/blending-bias_19_0.png
Shear measurement results : g1 =  0.11 , g2 =  0.28
Relative error : δg1 =  9.585 %, δg2 =  39.905 %

12.3.2.4. Deblending

One of the ways of dealing with the problem is called “deblending”, ie trying to separate the images of the different galaxies and perform shape measurement on those images. We explore how deblending can be carried out in Deblending.

12.3.3. Quantifying blending bias

To date, many weak lensing analysis have been carried using only very basic deblending solutions, eg Source Extractor, or rejecting altogether blended samples. However, this approach will become impractical with next generation surveys, as they are deeper than the previous one, resulting in more galaxies being detected leading to more blending ; according to [BoschArmstrongBickerton+18], 58% of the HSC detections correspond to blended objects. A key recurrent issue with blending is the fact that it is impossible to detect all the blends ; this means that one should calibrate not only for detected blends but also for unrecognized blends, adding to the difficulty of the problem ; there is also evidence of bias induced by blended objects which have higher magnitude than the threshold and hence are not detected.

Some papers have analyzed the bias due to blending :

  • [DawsonSchneiderTysonJee16] has tried to estimate the impact of blending in real data from HSC, using HST higher resolution data to identify ambiguous blends (when only one object is detected). They estimate, using conservative assumptions, that those undetected blends increase the measured shear variance by (1) increasing the variance of the galaxy ellipticities and (2) decreasing the total number of observed objects. They find an increase in the shear noise of 14% when considering galaxies up to magnitude 27, and by 7% when limiting the magnitude to 25.3 (corresponding to the LSST “Gold” Sample).

  • [MLL+18] presents calibration experiments for HSC, using HST data from the COSMOS catalog for simulating HSC images. Using two different samples, one removing the neighbours and the other keeping them, they identify a multiplicative bias due to blending ranging from 4% to 10%. They also show how cutting at a higher magnitude increases the multiplicative bias due to blending.

  • [SBMJ20] present a modification of Metacalibration to calibrate specifically for blending bias, called Metadetection ; in the paper they obtain a 3.5% for LSST Y10 simulations using the standard Metacalibration, and only 0.15% using Metadetection . They also measured the multiplicative bias as a function of the distance separating the blended galaxies ; they report that the main point of failure both for Metacalibration and Metadetection is at intermediate distances (1.5 arcseconds), when the blending is ambiguous, meaning that the objects are sometimes identified as two objects and sometimes only one by the detection algorithm (Metadetection is still much better than Metacalibration at this scale). Please note that the simulations from this paper use the same shear for all the sources and thus do not account for the possible bias due to different shears for different redshifts.

  • [MBM+20] establishes a formal framework for dealing with blended galaxies at different redshifts, assuming that they may have different shear. In constant shear simulation, using Metacalibration for shape measurement, they report a multiplicative bias of around 2% in their fiducial simulation, dropping to 0.4% in a simulation where no objects are blended. They also show that a having blended objects with a different shear (corresponding to the case where the objects are at different redshifts) induce a bias in the Metacalibration measurement.

12.3.4. Selection bias due to blending

One of the ways to avoid dealing with blended galaxies is to reject them before carrying shape measurement. While we have seen in the previous section that this is not a good idea for new surveys due to the very large proportion of blended galaxies, it may also be a source of bias in previous weak lensing analysis, inducing a selection bias.

Namely, [HartlapJHilbertSSchneiderPHildebrandtH11] describe two effects of rejecting blended galaxies. First, the redshift distribution is modified : galaxies at low redshift tend to be larger and are thus more prone to blending, meaning that rejecting them leads to a shift in the distribution towards higher redshifts. The second effect is more severe : the rejection of galaxy is density dependent, as an higher density of matter tends to lead to more galaxies, and thus more blending and more rejected galaxies. The opposite is true in underdense regions. This can induce a bias as high as a few percents on shear correlation statistics, potentially having large consequences on the results of weak lensing analysis.