Adding erosion step to source density map creation

[cmurray-astro]

Addressing issue #693, I added an "erode_boundary" parameter to the make_source_dens_map function in create_background_density_map. This (I think) will mean that the source density is computed for the same region that ASTs are placed in... in other words, the RA/Dec grid used for computing source density across each field will be the same (i.e., generated using the full catalog), but the value of the source density in each pixel is based on the sources that overlap with the eroded region.

Very much looking forward to thoughts from @galaxyumi and @karllark!

[cmurray-astro]

this was admittedly a premature PR. working on bugs!

[cmurray-astro]

I changed the strategy so that you can pass an optional parameter to create_background_density_map.main_make_map() called "erode_boundary" (which can be set to the same value as the beast_settings.txt parameter "ast_erode_selection_region", e.g.) which will be used to determine the RA/Dec grid on which the source density map is created (previously, the grid used the min/max values in the catalog).

For Scylla SMC_6 this results in a change in SD maps from the first to the second plots:

[lgtm-com[bot]]

This pull request introduces 1 alert when merging 32f25587d5576862a6924f2efb139f441fe6a74c into 2f67e824ea1a4d2dbb7b32e209db7a17b3a8e624 - view on LGTM.com

new alerts:

• 1 for Unused import

[galaxyumi]

Thanks @cmurray-astro for stepping out for this!

It is not super clear to me why the higher SD bins still appear on the boundary. Maybe they are very real this time?

[cmurray-astro]

Yes -- the troubling thing I've found is the dynamic range in source density depends strongly on the selected parameters for the map (pixel size (e.g., 5''), number of pixels, and now erosion boundary)

[karllark]

Is there something at the boundary that causes high values - maybe in the computation of the area of the "spatial pixels" at the edges? Maybe we just should not do any regions that are not fully covered inside the erosion boundary? Or maybe I don't really understand this issue - if so ignore this comment entirely.

[cmurray-astro]

these are good questions! @christinawlindberg has graciously agreed to test this PR and check how the pixel coverages are affected :)

[galaxyumi]

@cmurray-astro Yes, the dynamic range in source density also strongly depends on the magnitude range over which the source density is calculated. As long as the relative source density profile shape (# of pixels vs. source density) looks more or less the same, the varying absolute dynamic range, which depends on many factors, shouldn't be a problem. In other words, we should not rely on the absolute source density values. From conversation with @karllark, he pointed out that we don't know the underlying true source densities anyway without a perfect completeness correction. The key might be not reducing the number of source density bins by too much just based on the dynamic range in SD. For example, for one of the LUVIT galaxies, I had two (actually tried many more) versions of SD maps with different combinations of SD map parameters. I noticed that the dynamic range of one SD map goes from 0 to 15 while that of the other map goes from 0 to 7. I wondered if I want to reduce the number of SD bins for the second case because of a narrow dynamic range, and my answer was no. This was mainly because the way of distributing input ASTs across the image. To better explore the observed stellar density profile, to me it was more reasonable to keep the same number SD bins for those two cases. Does it make sense to you as well?

For the moment, my conclusion about creating a SD map has been reached to the point where I might want to calculate SD twice. First, create an initial SD map that well reflects the observed 2d stellar density profile just to place input ASTs. Once we compute the completeness from real AST results, we can create the final SD map, which will be used to split catalogs, compute noise models, and fit SEDs. In short, I won't worry too much about getting a precise SD map in the step of generating input ASTs. But, of course, we want to make sure the eroding is happening in the expected way!

[christinawlindberg]

I looked at where the SD map tiles/pixels were located for high SD ( > 2.32) and it looks like all the tiles in that SD bin are edge cases. Below, I've created plots for each of the tiles, and it looks like these elevated SD values are just a side effect of the tiny areas.

If I knew the indices for all of the edge tiles, I would suspect that there are also a lot of tiles with lower than average source density.

To avoid these cases, we thought it might be best to only compute the source density in tiles with a fractional overlap area > = 10%. @cmurray-astro Has an updated source density plot with this implemented and the SD range is a lot smaller!

[cmurray-astro]

yep! excluding pixels in the source density map grid which overlap with the catalog by less than 10% results in this:

[cmurray-astro]

@galaxyumi I agree with your conclusions about considering SD twice -- first for placing input ASTs, and then for splitting/fitting. Based on Christina's analysis, I also think it makes sense to exclude pixels which have <10% coverage by sources, because our measurement of source density for these is based on very small number statistics. I'm wondering what you both think (@karllark too) -- hard code a limit or add another settable parameter?

[karllark]

@christinawlindberg : if you are happy with the revised PR, please formally approve.