How to use the healpy.UNSEEN function in healpy

thb,phb = H.pix2ang(nside,arange(npix))
glat=90.-thb*180/pi

w=abs(glat)<10
sig[w] += 3.*randn(w.sum())+10

H.mollview(sig)

print 'ini: mono=%s, dipole=%s'%(c,d)
mono,dipole = H.pixelfunc.fit_dipole(sig)
print 'fit: mono=%s, dipole=%s'%(mono,dipole)

resfit = dot(vec.T,dipole)+mono

diff = sig.copy()
diff[sig != H.UNSEEN] -= resfit[sig != H.UNSEEN]

H.mollview(diff)

# use remove_dipole
m2 = H.pixelfunc.remove_dipole(sig)
m3 = H.pixelfunc.remove_monopole(sig)

H.mollview(m2)
H.mollview(m3)

def mergeSparseHealpixMaps(infiles, outfile=None,
                           pix_data_extension='PIX_DATA',
                           pix_field='PIX',
                           distance_modulus_extension='DISTANCE_MODULUS',
                           distance_modulus_field='DISTANCE_MODULUS',
                           default_value=healpy.UNSEEN):
    """
    Use the first infile to determine the basic contents to expect for the other files.
    """
    # Setup
    if isinstance(infiles,basestring): infiles = [infiles]
    
    distance_modulus_array = None
    pix_array = []
    data_dict = {}

    reader = pyfits.open(infiles[0])
    nside = reader[pix_data_extension].header['NSIDE']

    for ii in range(0, len(reader)):
        if reader[ii].name == distance_modulus_extension:
            distance_modulus_array = reader[distance_modulus_extension].data.field(distance_modulus_field)

Parameters
    ----------
    in_map : np.array
        A valie healpix map
    nside_out : int
        The desired resolution to convert in_map to
    UNSEEN2nan : bool (True)
        If True, convert any hp.UNSEEN values to np.nan
    """
    current_nside = hp.npix2nside(np.size(in_map))
    if current_nside != nside_out:
        out_map = hp.ud_grade(in_map, nside_out=nside_out)
    else:
        out_map = in_map
    if UNSEEN2nan:
        out_map[np.where(out_map == hp.UNSEEN)] = np.nan
    return out_map

import lsst.sims.featureScheduler as fs
from lsst.sims.speedObservatory import Speed_observatory
import healpy as hp


if __name__ == "__main__":

    survey_length = 365.25  # days
    # Define what we want the final visit ratio map to look like
    target_map = fs.standard_goals()['r']
    filtername = 'r'

    bfs = []
    bfs.append(fs.M5_diff_basis_function(filtername=filtername, teff=False))
    bfs.append(fs.Target_map_basis_function(target_map=target_map, filtername=filtername,
                                            out_of_bounds_val=hp.UNSEEN))
    bfs.append(fs.Quadrant_basis_function(quadrants='S', azWidth=15., maxAlt=75.))
    bfs.append(fs.Slewtime_basis_function(filtername=filtername))

    weights = np.array([1., 0.2, 1., 4.])
    surveys = []
    surveys.append(fs.Greedy_survey_fields(bfs, weights, block_size=1, filtername=filtername))
    surveys.append(fs.Pairs_survey_scripted([], []))
    scheduler = fs.Core_scheduler(surveys)

    observatory = Speed_observatory()
    observatory, scheduler, observations = fs.sim_runner(observatory, scheduler,
                                                         survey_length=survey_length,
                                                         filename='one_filter_south_pairs.db',
                                                         delete_past=True)

Return pointings for a block of sky
        """
        if not self.reward_checked:
            reward_smooth = self.calc_reward_function()
        else:
            reward_smooth = self.reward_smooth

        # Pick the top healpixels to observe
        reward_max = np.where(reward_smooth == np.max(reward_smooth))[0].min()
        unmasked = np.where(self.reward_smooth != hp.UNSEEN)[0]
        selected = np.where(reward_smooth[unmasked] >= np.percentile(reward_smooth[unmasked],
                                                                     self.percentile_clip))
        selected = unmasked[selected]

        to_observe = np.empty(reward_smooth.size, dtype=float)
        to_observe.fill(hp.UNSEEN)
        # Only those within max_region_size of the maximum
        max_vec = hp.pix2vec(self.nside, reward_max)
        pix_in_disk = hp.query_disc(self.nside, max_vec, self.max_region_size)

        # Select healpixels that have high reward, and are within
        # radius of the maximum pixel
        # Selected pixels are above the percentile threshold and within the radius
        selected = np.intersect1d(selected, pix_in_disk)
        if np.size(selected) > self.block_size:
            order = np.argsort(reward_smooth[selected])
            selected = selected[order[-self.block_size:]]

        to_observe[selected] = self.extra_features[0].feature[selected]

        # Find the pointings that observe the given pixels, and minimize the cross-correlation
        # between pointing overlaps regions and co-added depth

self.reward = 0
            indx = np.arange(hp.nside2npix(self.nside))
            # keep track of masked pixels
            mask = np.zeros(indx.size, dtype=bool)
            for bf, weight in zip(self.basis_functions, self.basis_weights):
                basis_value = bf(indx=indx)
                mask[np.where(basis_value == hp.UNSEEN)] = True
                if hasattr(basis_value, 'mask'):
                    mask[np.where(basis_value.mask == True)] = True
                self.reward += basis_value*weight
                # might be faster to pull this out into the feasabiliity check?
                if hasattr(self.reward, 'mask'):
                    indx = np.where(self.reward.mask == False)[0]
            self.reward[mask] = hp.UNSEEN
            self.reward.mask = mask
            self.reward.fill_value = hp.UNSEEN
            # inf reward means it trumps everything.
            if np.any(np.isinf(self.reward)):
                self.reward = np.inf

            if self.smoothing_kernel is not None:
                self.smooth_reward()

            # With the full reward map, now apply the masks and see which one is the winner
            potential_rewards = []
            potential_pointings = []
            potential_reward_maps = []
            for alt_az_mask in self.alt_az_masks:
                # blank reward map
                reward_map = self.reward*0 + hp.UNSEEN
                # convert the alt,az mask to ra,dec like the reward
                alt = self.extra_features['altaz'].feature['alt']

self.reward.mask = mask
            self.reward.fill_value = hp.UNSEEN
            # inf reward means it trumps everything.
            if np.any(np.isinf(self.reward)):
                self.reward = np.inf

            if self.smoothing_kernel is not None:
                self.smooth_reward()

            # With the full reward map, now apply the masks and see which one is the winner
            potential_rewards = []
            potential_pointings = []
            potential_reward_maps = []
            for alt_az_mask in self.alt_az_masks:
                # blank reward map
                reward_map = self.reward*0 + hp.UNSEEN
                # convert the alt,az mask to ra,dec like the reward
                alt = self.extra_features['altaz'].feature['alt']
                az = self.extra_features['altaz'].feature['az']
                mask_to_radec = hp.get_interp_val(alt_az_mask, np.pi/2 - alt, az)
                # Potential healpixels to observe, unmasked in alt,az and reward function
                good = np.where((self.reward != hp.UNSEEN) & (mask_to_radec > 0.))
                reward_map[good] = self.reward[good]
                potential_reward_maps.append(reward_map)
                # All the potential pointings in the alt,az block
                fields = self.hp2fields[self.hpids[good]]
                # compute a reward for each of the potential fields
                ufields = np.unique(fields)
                # Gather all the healpixels that are available
                observed_hp = np.isin(self.hp2fields, ufields)
                # now to bin up the reward for each field pointing
                # Use the max reward, I think summing ran into issues where

Return pointings for a block of sky
        """
        if not self.reward_checked:
            reward_smooth = self.calc_reward_function()
        else:
            reward_smooth = self.reward_smooth

        # Pick the top healpixels to observe
        reward_max = np.where(reward_smooth == np.max(reward_smooth))[0].min()
        unmasked = np.where(self.reward_smooth != hp.UNSEEN)[0]
        selected = np.where(reward_smooth[unmasked] >= np.percentile(reward_smooth[unmasked],
                                                                     self.percentile_clip))
        selected = unmasked[selected]

        to_observe = np.empty(reward_smooth.size, dtype=float)
        to_observe.fill(hp.UNSEEN)
        # Only those within max_region_size of the maximum
        max_vec = hp.pix2vec(self.nside, reward_max)
        pix_in_disk = hp.query_disc(self.nside, max_vec, self.max_region_size)

        # Select healpixels that have high reward, and are within
        # radius of the maximum pixel
        # Selected pixels are above the percentile threshold and within the radius
        selected = np.intersect1d(selected, pix_in_disk)
        if np.size(selected) > self.block_size:
            order = np.argsort(reward_smooth[selected])
            selected = selected[order[-self.block_size:]]

        to_observe[selected] = self.extra_features[0].feature[selected]

        # Find the pointings that observe the given pixels, and minimize the cross-correlation
        # between pointing overlaps regions and co-added depth

cls_s=hp.read_map("cont_lss_star_ns%d.fits"%nside_lss,field=0,verbose=False)
    cwp_q_s,cwp_u_s=hp.read_map("cont_wl_psf_ns%d.fits"%nside_lss,field=[0,1],verbose=False)
    cws_q_s,cws_u_s=hp.read_map("cont_wl_ss_ns%d.fits"%nside_lss,field=[0,1],verbose=False)
    tonull=np.where(msk_s==0)
    plt.figure(figsize=(12,13))
    fs_or=matplotlib.rcParams['font.size']
    ax=plt.subplot(421); mpp=msk_s.copy();# mpp[tonull]=hp.UNSEEN
    hp.mollview(mpp,hold=True,notext=True,coord=['G','C'],title='Mask',cbar=False,min=0,max=1.09)
    ax=plt.subplot(422); mpp=msk_s*mpt_s; mpp[tonull]=hp.UNSEEN
    matplotlib.rcParams.update({'font.size':14})
    hp.mollview(mpp,hold=True,notext=True,coord=['G','C'],title='$\\delta$',cbar=False)
    ax=plt.subplot(423); mpp=msk_s*mpq_s; mpp[tonull]=hp.UNSEEN
    hp.mollview(mpp,hold=True,notext=True,coord=['G','C'],title='$\\gamma_1$',cbar=False)
    ax=plt.subplot(424); mpp=msk_s*mpu_s; mpp[tonull]=hp.UNSEEN
    hp.mollview(mpp,hold=True,notext=True,coord=['G','C'],title='$\\gamma_2$',cbar=False)
    ax=plt.subplot(425); mpp=msk_s*cld_s; mpp[tonull]=hp.UNSEEN
    matplotlib.rcParams.update({'font.size':fs_or})
    hp.mollview(mpp,hold=True,notext=True,coord=['G','C'],title='Dust',cbar=False)
    ax=plt.subplot(426); mpp=msk_s*cls_s; mpp[tonull]=hp.UNSEEN
    hp.mollview(mpp,hold=True,notext=True,coord=['G','C'],title='Stars',cbar=False)
    ax=plt.subplot(427); mpp=msk_s*cwp_q_s; mpp[tonull]=hp.UNSEEN
    hp.mollview(mpp,hold=True,notext=True,coord=['G','C'],title='PSF',cbar=False)
    ax=plt.subplot(428); mpp=msk_s*cws_u_s; mpp[tonull]=hp.UNSEEN
    hp.mollview(mpp,hold=True,notext=True,coord=['G','C'],title='Small-scale',cbar=False)
    plt.savefig("../plots_paper/maps_lss_sph.pdf",bbox_inches='tight')
    plt.show()

if (o.whichfig==1) or (o.whichfig==-1) :
    #Plotting LSS power spectra
    l,cltt,clee,clbb,clte,nltt,nlee,nlbb,nlte=np.loadtxt("cls_lss.txt",unpack=True)

    plt.figure()

cross_corr : float
            If return_pointings_map is False, return the sum of the pointing map multipled
            with the
        """
        # XXX-check the nside

        # Unpack the x variable
        ra_rot = x[0]
        dec_rot = x[1]
        im_rot = x[2]

        # Rotate pointings to desired position
        final_ra, final_dec = rotate_ra_dec(self.ra, self.dec, ra_rot, dec_rot, init_rotate=im_rot)
        # Find the number of observations at each healpixel
        obs_map = self.p2hp_search(final_ra, final_dec)
        good = np.where(self.inmap != hp.UNSEEN)[0]

        if return_pointings_map:
            obs_indx = self.p2hp_search(final_ra, final_dec, stack=False)
            good_pointings = np.array([True if np.intersect1d(indxes, good).size > 0
                                      else False for indxes in obs_indx])
            if True not in good_pointings:
                raise ValueError('No pointings overlap requested pixels')
            obs_map = self.p2hp(final_ra[good_pointings], final_dec[good_pointings])
            return final_ra[good_pointings], final_dec[good_pointings], obs_map
        else:
            # If some requested pixels are not observed
            if np.min(obs_map[good]) == 0:
                return np.inf
            else:
                result = np.sum(self.inmap[good] *
                                obs_map[good])/float(np.sum(self.inmap[good] + obs_map[good]))