<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Basic demonstration of creating and using masks for bright sources
Author
Step1:  You may have to set up your $CSCRATCH environment variable so that Python can find it, e.g.
Step2:  These are some circular regions that could be masked, but let's check the capability also for elliptical regions
Step3:  Note that the ellipticity components, here, which are defined at the bottom of the page at, e.g., http
Step4:  Creating a file of targets and masking at the command line
Now let's compile a file of targets from the same set of sweeps. Again, I'll dump the targets file in your $CSCRATCH 
directory. This took just under 4 minutes, for me
Step5:  Masking targets
Our targets will already have DESI_TARGET bits set (indicating whether something is an LRG, ELG, QSO etc.) but we can update those bits to indicate which targets are in or are not in a mask. The IN_BRIGHT_OBJECT bit in DESI_TARGET indicates whether something is in a mask.
Step6:  Let's plot which objects are in masks and which are not, against the backdrop of the mask (in a small region of the sky)
Step7:  Note that the BADSKY locations are just outside the perimeter of the masks, and are quite obvious in the plot.
Masking points from a random catalog
The brightmask.set_target_bits() function wraps the main masking code in brightmask, called is_in_bright_mask(). The is_in_bright_mask() function can be used to compare coordinate locations to a mask and return which objects are in the mask (within the IN_RADIUS of the mask) and which objects are close to the mask (within the NEAR_RADIUS of the mask).
Let's create a random catalog over the small area of sky that we've been considering (154.8<sup>o</sup> < RA < 155.1<sup>o</sup>) and (19.7<sup>o</sup> < Dec < 20.0<sup>o</sup>)
Step8:  Now let's mask that random catalog
Step9:  and plot the random points that are and are not in the mask, both for the IN_RADIUS and the NEAR_RADIUS
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import matplotlib.pyplot as plt
import os
import numpy as np
import fitsio
from desitarget import desi_mask, brightmask
Explanation: Basic demonstration of creating and using masks for bright sources
Author: Adam D. Myers, University of Wyoming
Getting Started
Everything should work fine using a DESI kernel loaded from NERSC's jupyter-dev server. See the instuctions under Getting Started at:
https://github.com/desihub/tutorials/blob/master/Intro_to_DESI_spectra.ipynb
If you're running locally on your own machine, you may also need to set up your own development kernel as described here:
https://desi.lbl.gov/trac/wiki/Computing/JupyterAtNERSC
under defining your own kernel.
Preliminaries
Now, retrieve a sample of the sweeps files, say, everything in DR5 that starts sweep-150. I'll assume you put this in your $CSCRATCH directory at NERSC, as follows:
mkdir $CSCRATCH/sweep150
cp /global/project/projectdirs/cosmo/data/legacysurvey/dr5/sweep/5.0/sweep-150* $CSCRATCH/sweep150
Creating a mask based on bright sources from a sweep-like file or files at the command line
To create a bright source mask by scraping bright sources (of all types) from the sweeps:
make_bright_mask $CSCRATCH/sweep150 $CSCRATCH/sourcemask150.fits
this mask will use default settings and create masks for everything that is brighter than 10th magnitude in each of the g, r and z bands. This took just under 1.5 minutes, for me, on edison, and returned the following informative message:
INFO:make_bright_mask:48:&lt;module&gt;: wrote a file of 1292 masks to $CSCRATCH/sourcemask150.fits
You can also customize the limiting magnitude of source to be masked, e.g.:
make_bright_mask $CSCRATCH/sweep150 --bands GZ --maglim 9,10 $CSCRATCH/blat.fits
will mask any source that is brighter than g < 9 or z < 10.
This took about 1.5 minutes for me on edison and returned the following message:
INFO:make_bright_mask:48:&lt;module&gt;: wrote a file of 1164 masks to $CSCRATCH/blat.fits
Plotting masks
Let's examine the mask and plot it. Note that NEAR_RADIUS is a warning that you're near a mask (which could be useful for large-scale structure analyses) whereas IN_RADIUS is the radius to which DESI will certainly not place any fibers (for fear of saturating adjacent sources).
End of explanation
os.environ["CSCRATCH"] = '/global/cscratch1/sd/adamyers'
sourcemask = fitsio.read("$CSCRATCH/sourcemask150.fits")
brightmask.plot_mask(sourcemask,limits=[151,150,1,2])
brightmask.plot_mask(sourcemask,limits=[151,150,1,2],radius="NEAR_RADIUS")
Explanation: You may have to set up your $CSCRATCH environment variable so that Python can find it, e.g.:
End of explanation
from desitarget.brightmask import _rexlike
from desitarget.cuts import _psflike
rex_or_psf = _rexlike(sourcemask["TYPE"]) |  _psflike(sourcemask["TYPE"])
wcircle = np.where(rex_or_psf)
wellipse = np.where(~rex_or_psf)
sourcemask[wcircle][20:25]
sourcemask[wellipse][20:25]
Explanation: These are some circular regions that could be masked, but let's check the capability also for elliptical regions:
End of explanation
brightmask.plot_mask(sourcemask,limits=[155.1,154.8,19.7,20.0],radius="NEAR_RADIUS")
Explanation: Note that the ellipticity components, here, which are defined at the bottom of the page at, e.g., http://legacysurvey.org/dr5/catalogs/ are 0.0 for circular masks. 
Let's plot a region that has a mix of circular and elliptical masks:
End of explanation
targs = fitsio.read("$CSCRATCH/targs150.fits")
print(len(targs))
print(len(np.where( (targs["DESI_TARGET"] & desi_mask.BADSKY) != 0 )[0]))
targs = brightmask.append_safe_targets(targs,sourcemask)
print(len(targs))
print(len(np.where( (targs["DESI_TARGET"] & desi_mask.BADSKY) != 0 )[0]))
w = np.where( (targs["DESI_TARGET"] & desi_mask.BADSKY) != 0 )
badskies= targs[w]
brightmask.plot_mask(sourcemask,show=False)
plt.axis([155.1,154.8,19.7,20.0])
plt.plot(badskies["RA"],badskies["DEC"],'k,')
plt.xlabel('RA (o)')
plt.ylabel('Dec (o)')
plt.show()
Explanation: Creating a file of targets and masking at the command line
Now let's compile a file of targets from the same set of sweeps. Again, I'll dump the targets file in your $CSCRATCH 
directory. This took just under 4 minutes, for me:
select_targets $CSCRATCH/sweep150 $CSCRATCH/targs150.fits
and informed me that:
INFO:select_targets:73:&lt;module&gt;: 1741104 targets written to /global/cscratch1/sd/adamyers/targs150.fits
Note that at NERSC, the whole masking procedure can be applied to the targets in one line of code via, e.g.:
select_targets $CSCRATCH/sweep150 $CSCRATCH/targsblat150.fits --mask $CSCRATCH/sourcemask150.fits
Generating BADSKY ("SAFE") locations
First, let's generate BADSKY ("SAFE") locations around the periphery of each mask and plot them against the backdrop of the mask. Note that drstring, if set, is used to assign the Data Release (DR) bit to BADSKY targets and to assign BADSKY targets an OBJID that is higher than any target in a given brick of that DR. (That's the part that requires a DR survey-brick file, though).
End of explanation
desi_mask
dt = brightmask.set_target_bits(targs,sourcemask)
inmask = np.where( (dt & desi_mask.IN_BRIGHT_OBJECT) != 0) 
masked = targs[inmask]
notinmask = np.where( (dt & desi_mask.IN_BRIGHT_OBJECT) == 0) 
unmasked = targs[notinmask]
Explanation: Masking targets
Our targets will already have DESI_TARGET bits set (indicating whether something is an LRG, ELG, QSO etc.) but we can update those bits to indicate which targets are in or are not in a mask. The IN_BRIGHT_OBJECT bit in DESI_TARGET indicates whether something is in a mask.
End of explanation
brightmask.plot_mask(sourcemask,show=False)
plt.axis([155.1,154.8,19.7,20.0])
plt.xlabel('RA (o)')
plt.ylabel('Dec (o)')
plt.plot(masked["RA"],masked["DEC"],'kx')
plt.plot(unmasked["RA"],unmasked["DEC"],'r.')
plt.show()
Explanation: Let's plot which objects are in masks and which are not, against the backdrop of the mask (in a small region of the sky):
End of explanation
from numpy.random import random
Nran = 100000
rancat = np.zeros(Nran, dtype=[('RA', '>f8'), ('DEC', '>f8')])
rancat["RA"] = 154.8+0.3*(random(Nran))
rancat["DEC"] = np.degrees(np.arcsin(np.sin(np.radians(20))-random(Nran)*0.05))
Explanation: Note that the BADSKY locations are just outside the perimeter of the masks, and are quite obvious in the plot.
Masking points from a random catalog
The brightmask.set_target_bits() function wraps the main masking code in brightmask, called is_in_bright_mask(). The is_in_bright_mask() function can be used to compare coordinate locations to a mask and return which objects are in the mask (within the IN_RADIUS of the mask) and which objects are close to the mask (within the NEAR_RADIUS of the mask).
Let's create a random catalog over the small area of sky that we've been considering (154.8<sup>o</sup> < RA < 155.1<sup>o</sup>) and (19.7<sup>o</sup> < Dec < 20.0<sup>o</sup>):
End of explanation
inmask, nearmask = brightmask.is_in_bright_mask(rancat,sourcemask)
masked = rancat[np.where(inmask)]
notmasked = rancat[np.where(~inmask)]
near = rancat[np.where(nearmask)]
notnear = rancat[np.where(~nearmask)]
Explanation: Now let's mask that random catalog:
End of explanation
brightmask.plot_mask(sourcemask,show=False)
plt.axis([155.1,154.8,19.7,20.0])
plt.xlabel('RA (o)')
plt.ylabel('Dec (o)')
plt.plot(masked["RA"],masked["DEC"],'r.')
plt.plot(notmasked["RA"],notmasked["DEC"],'g,')
plt.show()
brightmask.plot_mask(sourcemask,show=False,radius="NEAR_RADIUS")
plt.axis([155.1,154.8,19.7,20.0])
plt.xlabel('RA (o)')
plt.ylabel('Dec (o)')
plt.plot(near["RA"],near["DEC"],'r.')
plt.plot(notnear["RA"],notnear["DEC"],'g,')
plt.show()
Explanation: and plot the random points that are and are not in the mask, both for the IN_RADIUS and the NEAR_RADIUS:
End of explanation
<END_TASK>