lsst_desi_overlap.py

import fitsio
import healpy
import os
import numpy as np
import pylab

nside = 512
npix = healpy.nside2npix(nside)

# download DESI
if not os.path.exists('desi-healpix-weights.fits'):
    os.system("wget https://desi.lbl.gov/svn/code/desimodel/trunk/data/footprint/desi-healpix-weights.fits")

if not os.path.exists('opsim_nvisits_g.fits'):
    os.system("wget http://lambda.gsfc.nasa.gov/data/footprint-maps/opsim_nvisits_g.fits.gz")
    os.system("gunzip opsim_nvisits_g.fits.gz")

f = fitsio.FITS("desi-healpix-weights.fits")
desi = f[0][:]
f.close()
desi = healpy.ud_grade(desi, nside_out=nside, order_in='nest', order_out='ring')

# read lsst
lsst = healpy.read_map("opsim_nvisits_g.fits")
lsst = healpy.ud_grade(lsst, nside_out=nside, order_in='ring', order_out='ring')

# get region in both
overlap = (lsst>0) & (desi > 0)

# print region
frac = overlap.sum() / overlap.size
sq_deg = frac * 41253.
print(f'Sky fraction in non-galactic overlap = {frac:.2f}')
print(f'Area in non-galactic overlap = {sq_deg:.1f}')

# plot overlaps
healpy.mollview(overlap)
healpy.graticule()
pylab.savefig('overlap.png')