Reweight area estimation in QA plots to better reflect the imaging survey footprint

[geordie666]

This PR creates the tools to better handle imaging survey areas for QA. It includes:

• A new module desitarget.imagefootprint to build HEALPix weight maps of the imaging footprint.
• An executable (bin) interface for that module to make weight maps from the command line.
• A loader in desitarget.io to resample maps to any HEALPix nside.
• Updates to desitarget.QA to handle the new imaging area weight maps.

The code was developed from an earlier script by Anand Raichoor. The main algorithm works by Monte Carlo sampling the Legacy Surveys coadd files for a given Data Release at the pixel-level. A large density of RA/Dec points are generated in order to determine the area-weighted fraction of locations that are in pixels that have one or more observations.

Generating 10million random points per sq. deg., the weight map can be produced at a rate of about 2500 bricks per minute across 32 processors on a single node on cori.

Here is a comparison of the DR4 QA plots before

and after

the adoption of the new imaging weights.

[moustakas]

This is an outstanding development, @geordie666, thank you! But can you please provide a couple examples of how you're calling the code differently? For example, it looks like desitarget.QA is still just reading the old desi-healpix-weights.fits file from desimodel.io.load_pixweight.

I'm asking because I use your pixel weights in the spectral mocks code to determine the target density and the downsampling factor from the mocks...

[geordie666]

@moustakas:

The mocks part of the code will still call the desimodel.footprint weights map. In desitarget.QA.mock_make_qa_plots you still have the lines of code:

    from desimodel import io, footprint
    pix = footprint.radec2pix(nside, targs["RA"], targs["DEC"])
    log.info('Calculated HEALPixel for targets and truths...t = {:.1f}s'
             .format(time()-start))

    #ADM set up the weight of each HEALPixel, if requested.
    weights = np.ones(len(targs))
    if weight:
        #ADM retrieve the map of what HEALPixels are actually in the DESI footprint
        pixweight = io.load_pixweight(nside)
        #ADM determine what HEALPixels each target is in, to set the weights
        fracarea = pixweight[pix]

but the real-imaging-data part of the code will call the imaging-based weights map, e.g., in desitarget.QA.make_qa_plots you now have:

    from desimodel import footprint
    pix = footprint.radec2pix(nside, targs["RA"], targs["DEC"])
    log.info('Calculated HEALPixel for each target...t = {:.1f}s'
             .format(time()-start))
    #ADM set up the weight of each HEALPixel, if requested.
    weights = np.ones(len(targs))
    if weight:
        #ADM load the imaging weights file
        if imaging_map_file is not None:
            from desitarget import io as dtio
            pixweight = dtio.load_pixweight(imaging_map_file,nside)

To make the weights map, follow the help message in the bin/make_imaging_weight_map script, for example:

make_imaging_weight_map /global/project/projectdirs/cosmo/data/legacysurvey/dr4/ $CSCRATCH/pixmaptemp.fits

To then use that weight map in the (real-imaging-data) QA pass it to the existing bin/run_target_qa script, e.g.:

run_target_qa /project/projectdirs/desi/target/catalogs/targets-dr4-0.20.0.fits /global/project/projectdirs/desi/www/users/adamyers/garbage -w $CSCRATCH/pixmaptemp.fits

That's the top-line information, at least, although it's nothing more than is in the help messages for those scripts. Let me know if you need more information.