Source code for desitarget.mock.sky

import numpy as np
import healpy as hp

def _random_theta_phi(nside, pix):
    theta, phi = hp.pix2ang(nside, pix, nest=True)
    dpix = np.sqrt(hp.nside2pixarea(nside))
    theta += np.random.uniform(-dpix/2, dpix/2, size=len(theta))
    phi += np.random.uniform(-dpix/2, dpix/2, size=len(phi)) * np.cos(np.pi/2 - theta)
    return theta % np.pi, phi % (2*np.pi)

[docs]def random_sky(nside=2048, allsky=True, tiles=None, maxiter=20, outfile=None): '''Returns sky locations within healpixels covering tiles Options: nside (int): healpixel nside; coverage is uniform at this scale allsky (bool): generate sky positions over the full sky tiles: DESI tiles to cover, for desimodel.footprint.tiles2pix(); only used if allsky=False. maxiter (int): maximum number of iterations to ensure coverage Generates sky locations that are more uniform than true random, such that every healpixel has a point within it. Note that this should *not* be used for mock randoms. nside=2048 corresponds to about half of a DESI positioner patrol area and results in ~18M sky locations over the full footprint. ''' if allsky: pix = np.arange( hp.nside2npix(nside) ) else: import desimodel.io import desimodel.footprint if tiles is None: tiles = desimodel.io.load_tiles() pix = desimodel.footprint.tiles2pix(nside, tiles) theta, phi = _random_theta_phi(nside, pix) #- there is probably a more clever way to do this, but iteratively #- replace points that scatter outside their original healpixel until #- all healpixels are covered for i in range(maxiter): skypix = hp.ang2pix(nside, theta, phi, nest=True) missing = np.in1d(pix, skypix, invert=True) ii = np.where(missing)[0] if len(ii) == 0: break tx, px = _random_theta_phi(nside, pix[ii]) theta[ii] = tx phi[ii] = px ra = np.degrees(phi) dec = 90 - np.degrees(theta) if outfile is not None: from astropy.table import Table keep = dec >= -30 out = Table() out['HPXPIXEL'] = pix[keep] out['RA'] = ra[keep] out['DEC'] = dec[keep] print('Writing {}'.format(outfile)) out.write(outfile, overwrite=True) return ra, dec, pix