Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def _hp2fieldsetup(self, ra, dec, leafsize=100):
"""Map each healpixel to nearest field. This will only work if healpix
resolution is higher than field resolution.
"""
if self.camera == 'LSST':
pointing2hpindx = hp_in_lsst_fov(nside=self.nside)
elif self.camera == 'comcam':
pointing2hpindx = hp_in_comcam_fov(nside=self.nside)
self.hp2fields = np.zeros(hp.nside2npix(self.nside), dtype=np.int)
for i in range(len(ra)):
hpindx = pointing2hpindx(ra[i], dec[i], rotSkyPos=0.)
self.hp2fields[hpindx] = i
def mapFromClm(clm, nside):
"""
Given an array of C_{lm} values, produce a pixel-power-map (non-Nested) for
healpix pixelation with nside
@param clm: Array of C_{lm} values (inc. 0,0 element)
@param nside: Nside of the healpix pixelation
return: Healpix pixels
Use real_sph_harm for the map
"""
npixels = hp.nside2npix(nside)
pixels = hp.pix2ang(nside, np.arange(npixels), nest=False)
h = np.zeros(npixels)
ind = 0
maxl = int(np.sqrt(len(clm))) - 1
for ll in range(maxl + 1):
for mm in range(-ll, ll + 1):
h += clm[ind] * real_sph_harm(mm, ll, pixels[1], pixels[0])
ind += 1
return h
if __name__ == "__main__":
from optparse import OptionParser
usage = "Usage: %prog [options] input"
description = "python script"
parser = OptionParser(usage=usage,description=description)
parser.add_option('-o','--outfile',default='allsky_results.png')
parser.add_option('-t','--targets',default=None)
parser.add_option('-c','--coord',default='GAL')
parser.add_option('-p','--proj',default='MOL',choices=['MOL','CAR'])
parser.add_option('-f','--field',default='LOG_LIKELIHOOD')
(opts, args) = parser.parse_args()
nside = pyfits.open(args[0])[1].header['NSIDE']
map = healpy.UNSEEN * numpy.ones( healpy.nside2npix(nside) )
pix,vals = ugali.utils.skymap.readSparseHealpixMap(args[0],opts.field,construct_map=False)
map[pix] = vals[0]
if opts.coord.upper() == "GAL":
coord = 'G'
elif opts.coord.upper() == "CEL":
coord = 'GC'
if opts.proj.upper() == "MOL":
healpy.mollview(map,coord=coord,xsize=1000)
elif opts.proj.upper() == "CAR":
healpy.cartview(map,coord=coord,xsize=1000)
else:
raise Exception("...")
healpy.graticule()
if opts.targets:
def calc_reward_function(self):
self.reward_count += 1
self.reward_checked = True
if self._check_feasability():
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)):
for ii in range(0, len(self.distance_modulus_array)):
if distance_modulus_index is not None and ii != distance_modulus_index:
continue
distance = ugali.utils.projector.distanceModulusToDistance(self.distance_modulus_array[ii])
if mode == 'ts':
title = 'Test Statistic (%.1f kpc)'%(distance)
map = healpy.UNSEEN * numpy.ones(healpy.nside2npix(self.nside))
map[self.pixels] = self.log_likelihood_sparse[ii]
index = map != healpy.UNSEEN
map[index] = 2. * map[index]
if mode == 'lim':
title = 'Upper Limit (%.1f kpc)'%(distance)
map = healpy.UNSEEN * numpy.ones(healpy.nside2npix(self.nside))
map[self.pixels] = self.richness_lim_sparse[ii]
ugali.utils.plotting.zoomedHealpixMap(title, map, lon_median, lat_median, radius, **kwargs)
#except:
# logger.warn("Couldn't open %s; will try again."%footfile)
# footprint = footfile
mag_column = self.config['catalog']['mag_%i_field'%field]
magerr_column = self.config['catalog']['mag_err_%i_field'%field]
# For simple maglims
release = self.config['data']['release'].lower()
band = self.config['catalog']['mag_%i_band'%field]
pixel_pix_name = 'PIX%i'%self.nside_pixel
data = fitsio.read(infile,columns=[pixel_pix_name])
#mask_pixels = numpy.arange( healpy.nside2npix(self.nside_mask), dtype='int')
mask_maglims = numpy.zeros( healpy.nside2npix(self.nside_mask) )
out_pixels = numpy.zeros(0,dtype='int')
out_maglims = numpy.zeros(0)
# Find the objects in each pixel
pixel_pix = data[pixel_pix_name]
mask_pix = ugali.utils.skymap.superpixel(pixel_pix,self.nside_pixel,self.nside_mask)
count = Counter(mask_pix)
pixels = sorted(count.keys())
pix_digi = numpy.digitize(mask_pix,pixels).argsort()
idx = 0
min_num = 500
signal_to_noise = 10.
magerr_lim = 1/signal_to_noise
for pix in pixels:
# Calculate the magnitude limit in each pixel
def _hp2fieldsetup(self, ra, dec, leafsize=100):
"""Map each healpixel to nearest field. This will only work if healpix
resolution is higher than field resolution.
"""
pointing2hpindx = hp_in_lsst_fov(nside=self.nside)
self.hp2fields = np.zeros(hp.nside2npix(self.nside), dtype=np.int)
for i in range(len(ra)):
hpindx = pointing2hpindx(ra[i], dec[i])
self.hp2fields[hpindx] = i
def calc_reward_function(self, conditions):
self.reward_checked = True
if self._check_feasibility(conditions):
self.reward = 0
indx = np.arange(hp.nside2npix(self.nside))
for bf, weight in zip(self.basis_functions, self.basis_weights):
basis_value = bf(conditions, indx=indx)
self.reward += basis_value*weight
if np.any(np.isinf(self.reward)):
self.reward = np.inf
else:
# If not feasable, negative infinity reward
self.reward = -np.inf
if self.smoothing_kernel is not None:
self.smooth_reward()
return self.reward_smooth
else:
return self.reward
def __init__(self, nside=None, out_of_bounds_val=np.nan, az_min=0., az_max=180.):
super(Mask_azimuth_basis_function, self).__init__(nside=nside)
self.az_min = int_rounded(np.radians(az_min))
self.az_max = int_rounded(np.radians(az_max))
self.out_of_bounds_val = out_of_bounds_val
self.result = np.ones(hp.nside2npix(self.nside))
from astroML.datasets import fetch_wmap_temperatures
#----------------------------------------------------------------------
# This function adjusts matplotlib settings for a uniform feel in the textbook.
# Note that with usetex=True, fonts are rendered with LaTeX. This may
# result in an error if LaTeX is not installed on your system. In that case,
# you can set usetex to False.
from astroML.plotting import setup_text_plots
setup_text_plots(fontsize=8, usetex=True)
#------------------------------------------------------------
# First plot an example pixellization
# Prepare the healpix pixels
NSIDE = 4
m = np.arange(hp.nside2npix(NSIDE))
print("number of pixels:", len(m))
# Plot the pixelization
fig = plt.figure(1, figsize=(5, 3.75))
hp.mollview(m, nest=True, title="HEALPix Pixels (Mollweide)", fig=1)
# remove colorbar: we don't need it for this plot
fig.delaxes(fig.axes[1])
#------------------------------------------------------------
# Next plot the wmap pixellization
wmap_unmasked = fetch_wmap_temperatures(masked=False)
# plot the unmasked map
fig = plt.figure(2, figsize=(5, 3.75))
hp.mollview(wmap_unmasked, min=-1, max=1, title='Raw WMAP data',