Add new skyhealpixs.py file

[araichoor]

This PR adds a new script skyhealpixs.py to use the lookup-table for stuck-sky positions, based on healpix-based files.
A companion PR should be done in fiberassign, to make use of that script.
See this commit: desihub/fiberassign@fa681fa (mostly adapting fiberassign/stucksky.py written by @dstndstn).

Development done with SKYHEALPIXS_DIR=/global/cfs/cdirs/desi/users/koposov/gaia_sky_mask (provided by @segasai).

It is twinning the skybricks.py script, with replacing the per-brick items by per-healpix-pixel items.
It assumes a SKYHEALPIXS_DIR environmental variable, similar to SKYBRICKS_DIR.
I would appreciate feedback on those two choices I made, i.e. is it ok? or would it be better to modify skybricks.py and keep using SKYBRICKS_DIR?
I chose to duplicate the information to minimize confusion / code complexity.

I assumed that the (ra, dec) -> (x, y) conversion works similarly, and kept the rounding of (x, y) to integers:

desitarget/py/desitarget/skyhealpixs.py

Lines 81 to 91 in e0307a1

	H, W = skymap.shape
	# create WCS object
	w = WCS(naxis=2)
	w.wcs.ctype = [hdr['CTYPE1'], hdr['CTYPE2']]
	w.wcs.crpix = [hdr['CRPIX1'], hdr['CRPIX2']]
	w.wcs.crval = [hdr['CRVAL1'], hdr['CRVAL2']]
	w.wcs.cd = [[hdr['CD1_1'], hdr['CD1_2']],
	[hdr['CD2_1'], hdr['CD2_2']]]
	x, y = w.wcs_world2pix(ras.flat[I], decs.flat[I], 0)
	x = np.round(x).astype(int)
	y = np.round(y).astype(int)

For sanity-check, I have run it on the 770k stuck-sky fibers of the 1.5k observed main tiles so far; I identify stuck-sky fibers with (OBJTYPE=="SKY") & (TARGETID<0).
Only 85 positions are returned as not-good, and those all are close to some bright star, so this looks ok to me; e.g.:

(https://www.legacysurvey.org/viewer-desi/?ra=213.3179&dec=-0.5157&catalog=tmpia0ogmvv#)

One caveat:
Reading those files is a bit slow, takes ~5s, whereas the per-brick with skybricks.py is much faster, I feel.
@schlafly : I don t know if that is an issue from the operations point-of-view.
If so, I can make a more precise diagnosis.
But already one option: the current healpix-split use nside=64, i.e. a pixel size of ~0.83 deg2; one could imagine that using nside=256, i.e. a pixel size of ~0.05 deg2, close to the size of a brick, could be a reasonable choice?

[coveralls]

Coverage decreased (-0.2%) to 57.116% when pulling bb0c13c on skyhealpixs_gaia into f7aef60 on master.

[segasai]

Thanks Anand, it looks good, but with a couple of comments.

[geordie666]

Don't forget to update the doc/changes.rst file!

[araichoor]

Thanks for reminding that, I have just updated the doc/changes.rst.

[sbailey]

I'm catching up on the context here: is this intended to replace skybricks, or be in addition to it? I see that /global/cfs/cdirs/desi/users/koposov/gaia_sky_mask is 36 GB while $SKYBRICKS_DIR=/global/cfs/cdirs/desi/target/skybricks/v3/ is 52 GB, which implies that this is in addition to skybricks, but maybe you've found a more clever packing. I'd really like to have one sky lookup that automatically works both on- and off- footprint rather than having to have fiberassign and desispec code to on/off footprint logic to know which one to call.

[segasai]

I would assume that the overlaps of healpixels are smaller than for the bricks. Also my bit-packing is I convert 4x4 0.25" single bit pixels into 16 bits rather than just counting how many of 4x4 pixels are occupied (which I would assume make it less compressible). Also I think gaia maps must be more empty, so therefore more easily compressed, so I don't know if the size difference is really a reflection of better packing or not.

One big advantage of having one consistent skymap pixelization is that one could in principle just combine those from multiple sources I.e. we have gaia sky map and decasls and wise sky map and we could just combine then pixel by pixel.

[sbailey]

@segasai confirming:

• do your skyhealpix already cover all sky, or not yet?
• if they cover all sky, are they only Gaia-based or do they use legacy survey data on-footprint when available?

I agree with the benefits of "one consistent skymap pixelization", but I don't want to switch to that until it also incorporates the on-footprint information that we already have, and I'd also like to avoid an intermediate state where we have one pixelation+library on-footprint and a different pixelation+library off-footprint.

[segasai]

• It covers all sky with dec>=-30 (to the best of my knowledge)
• It is purely Gaia based.
  generated by this https://github.com/segasai/gaia_sky_mask/

I guess I could postprocess my files by merging them with the current v3 skys. (that's not the ideal approach, but at least easy)

I was hoping that we would use pure gaia maps for all the backup tiles and the existing skies for non-backup and merge them later ... AFAIU the interface is basically the same. But I agree it'd be nice to avoid the state with two sets of skies.

[djschlegel]

Sergey's files cover the full sky at dec > -33 deg, plus some sky further south at RA~300 that may be unintentional. The biggest difference in file size is from Gaia (right image) masking far fewer sources than the DR9 masks (left image). But I'll note that gzip is ~30% more efficient at compressing these images than fpack. It would be convenient if these were in the same brick format such that they could be compared more easily, and then we could even combine them within the DR9 footprint, as suggested.

[segasai]

@dstndstn Do I understand correctly that your skies basically rely on blobs>-1 and nexp>0 from legacysurvey/metrics ?
I think I can probably quickly extract that to my pixels

[dstndstn]

My $0.02 is I think these healpix tiles are not as nice as the desiutil Brick tiles :)

Yes, my skybricks are defined by blobs > -1, or sum(nexp) == 0:
https://github.com/dstndstn/desi-misc/blob/main/desi-sky-locations.py#L68
https://github.com/dstndstn/desi-misc/blob/main/desi-sky-locations.py#L76

Regarding the compression: I'm using FPACK tiled compression: it's compressed in 256x256 blocks to allow for quicker reading of subimages, though I guess that the way it's used by fiberassign, that's not used. The tiling results in slightly worse compression also.

[segasai]

I've started to produce the healpix-shaped legacy-survey based masks here /global/homes/k/koposov/desi_koposov/make_sky/gaia_sky_mask/
they are basically the same as gaia-ones, but rely on the legacy data as Dustin did. The only thing that I didn't do here is the masking out of regions without the DECALS coverage (i.e. nexp==0) as it is not 100% obvious we want that if we want to AND combine the masks.

The masks can be directly compared (top left is gaia, top right is decals, bottom is a dss image)

[segasai]

I've finished creating the legacy survey sky healpixels here
/global/cfs/cdirs/desi/users/koposov/sky_mask/ls_data/
moved the gaia pixels here
/global/cfs/cdirs/desi/users/koposov/sky_mask/gaia_data/
They are now all generate in the same format/tiling.
The size of the legacy survey skys is 39gb
(and again I didn't mask out completely the out of footprint areas for DECals maps)

[sbailey]

Documenting discussion from the Monday survey ops call into this ticket:

Eddie pointed out that for backwards compatibility reasons, we will need to continue using existing skybricks code for main survey tiles. As David S noted, even if this new code exactly reproduced sky decisions on-footprint, tiles that hang off the footprint would still have different assignements because new skies would become available at the edge, thus changing the sky assignment bookkeeping and bumping.

i.e. ugh, it looks like fiberassign will need some if/then logic to chose which stucksky routine to use based upon flavor of tile.

Also, see desihub/desispec#1483 for the impact on offline pipeline operations for backup tiles.

[araichoor]

for what is worth, I ve already coded up in some fiberassign branch the if/then logic to pick gaia-based or ls-based look up table: desihub/fiberassign@fa681fa
I chose to add a new --lookup_sky_source argument to the (~core) fba_run function:

    parser.add_argument("--lookup_sky_source", required=False, default="ls",
                        choices={"ls", "gaia"},
                        help="Source for the look-up table for sky positions for stuck fibers:"
                        " 'ls': uses $SKYBRICKS_DIR; 'gaia': uses $SKYHEALPIXS_DIR (default=ls)")

this is set in fba_launch_io.launch_onetile_fa(), based on the program.

[schlafly]

for what is worth, I ve already coded up in some fiberassign branch the if/then logic to pick gaia-based or ls-based look up table: ***@***.*** <desihub/fiberassign@fa681fa> I chose to add a new --lookup_sky_source argument to the (~core) fba_run function: parser.add_argument("--lookup_sky_source", required=False, default="ls", choices={"ls", "gaia"}, help="Source for the look-up table for sky positions for stuck fibers:" " 'ls': uses $SKYBRICKS_DIR; 'gaia': uses $SKYHEALPIXS_DIR (default=ls)") this is set in fba_launch_io.launch_onetile_fa(), based on the program.

You probably have already done this, but I didn't immediately spot it: let's store this in the header so that the offline pipeline can duplicate. i.e., when desispec goes to verify whether positioners landed on sky, it can do something like stuck_on_sky(..., lookup_sky_source=hdr['SKY_SOURCE'])

…

[araichoor]

that is a very good suggestion, I ll do that (I haven t done it yet), thanks!

[araichoor]

For the record, the fiberassign companion PR mentioned in the first message is: desihub/fiberassign#416.