Sky brightness inaccurate due to units error

[ShrihanSolo]

This is for current version PyPI 0.0.2.3, on macOS Sonoma 14.5. The sky brightness input to deeplenstronomy is shown in the examples (https://deepskies.github.io/deeplenstronomy/Notebooks/ConfigFiles.html) to be:

sky_brightness: measured light contamination from the atmosphere in magnitude per square arcsecond

followed by,

SURVEY:
    PARAMETERS:
        BANDS: g,r,i,z,Y
        seeing: 0.9
        magnitude_zero_point: 30.0
        sky_brightness: 23.5
        num_exposures: 10

This indicates an input in magnitude/sq. arcsec. I believe the sky_brightness value is not manipulated and directly passed into lenstronomy.Util.data_util.bkg_noise() which takes input in counts/s/sq. arcsec, rather than magnitude/sq. arcsec. See function signature below (https://lenstronomy.readthedocs.io/en/latest/_modules/lenstronomy/Util/data_util.html):

def bkg_noise(
    readout_noise, exposure_time, sky_brightness, pixel_scale, num_exposures=1
):
    """Computes the expected Gaussian background noise of a pixel in units of
    counts/second.

    :param readout_noise: noise added per readout
    :param exposure_time: exposure time per exposure (in seconds)
    :param sky_brightness: counts per second per unit arcseconds square
    :param pixel_scale: size of pixel in units arcseonds
    :param num_exposures: number of exposures (with same exposure time) to be co-added
    :return: estimated Gaussian noise sqrt(variance)
    """
    exposure_time_tot = num_exposures * exposure_time
    readout_noise_tot = num_exposures * readout_noise**2  # square of readout noise
    sky_per_pixel = sky_brightness * pixel_scale**2
    sky_brightness_tot = exposure_time_tot * sky_per_pixel
    sigma_bkg = np.sqrt(readout_noise_tot + sky_brightness_tot) / exposure_time_tot
    return sigma_bkg

To test, I attach 3 images with only changing sky brightness in the .yaml file, and show that the sky brightness input is certainly acting like a counts/s input rather than a mag/sq. arcsec input.

Testing for a visible change in generated images with changing sky brightness required a sky_brightness: 10_000 or so before it could be seen, indicating inputs should be in flux. Entering a value of 0. or 1. sky brightness should produce a very bright image for magnitude/sq. arcsec but instead produces an image with no sky noise. Similarly, entering a negative value for the sky_brightness causes nan images, which should not occur for magnitude based inputs.

I also attach .yaml files which can be used to reproduce the same (uploaded as .txt at the end). The only difference in these files is the sky_brightness parameter.

Lastly, the default survey distributions provided in distributions.py (https://github.com/deepskies/deeplenstronomy/blob/9672080cd1a5792887341f4fb13c3b0bcf913600/deeplenstronomy/distributions.py) are also in mag/sq. arcsec. This means that any use of the default survey distributions such as des_sky_brightness causes() a too-low sky brightness input to replicate the survey.

Happy to provide more information or receive corrections if I have misunderstood this.

sky10000.txt
sky100.txt
sky1.txt

[bnord]

@Jasonpoh Did you say that you think this is a one-line fix? If so, do you have a guess about what that would be?

[Jasonpoh]

In this line of the code, the variable kwargs_single_band['sky_brightness'] should be converted from magnitude to counts per second per unit arcseconds squared first before being passed into data_util.bkg_noise.

This can be done by this lenstronomy function.

Basically, a hacky way to do this replacing the code snippet from line 273 with this:

#convert skybrightness from magnitude to counts per second per square arcsec
sky_brightness_cps = data_util.magnitude2cps(kwargs_single_band['sky_brightness'], kwargs_single_band['magnitude_zero_point'])

# generate image
image_sim = image_model.image(kwargs_lens_model_list, kwargs_source_list, kwargs_lens_light_list, kwargs_ps)
poisson = image_util.add_poisson(image_sim, exp_time=kwargs_single_band['exposure_time'])
sigma_bkg = data_util.bkg_noise(kwargs_single_band['read_noise'],
                     kwargs_single_band['exposure_time'],
                     sky_brightness_cps,
                    kwargs_single_band['pixel_scale'],
                    num_exposures=kwargs_single_band['num_exposures'])
bkg = image_util.add_background(image_sim, sigma_bkd=sigma_bkg)
image = image_sim + bkg + poisson

[bnord]

Thanks @Jasonpoh

@jsv1206 @ShrihanSolo @paxsonswierc This will affect our coming papers, so we should probably fix this soon. Do y'all have some time to work on this?

[paxsonswierc]

I implemented the fix as Jason wrote! It is now working as expected. I will push the branch with the fix today.

Here are lens from the same files Shrihan used:

Sky 1 --

Sky 100 --

Sky 10000 --

[bnord]

Thanks, Paxson. Please tag me in the pull request.