Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with Healpy; if not, write to the Free Software
# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
#
# For more information about Healpy, see http://code.google.com/p/healpy
#
import healpy as H
from numpy import *
from pylab import *
nside = 32
npix = H.nside2npix(nside)
vec = array(H.pix2vec(nside,arange(npix)))
d = 3*array([0.5,0.3,0.2])
#d /= norm(d)
c = 2.
sig = dot(vec.T,d)+c
sig += randn(npix)
sig[randint(sig.size,size=10000)] = H.UNSEEN
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)
def match_hpx_pixel(nside, nest, nside_pix, ipix_ring):
"""
"""
ipix_in = np.arange(12 * nside * nside)
vecs = hp.pix2vec(nside, ipix_in, nest)
pix_match = hp.vec2pix(nside_pix, vecs[0], vecs[1], vecs[2]) == ipix_ring
return ipix_in[pix_match]
def principal_axes(prob, distmu, distsigma, nest=False):
npix = len(prob)
nside = hp.npix2nside(npix)
bad = ~(np.isfinite(prob) & np.isfinite(distmu) & np.isfinite(distsigma))
distmean, diststd, _ = parameters_to_moments(distmu, distsigma)
mass = prob * (np.square(diststd) + np.square(distmean))
mass[bad] = 0.0
xyz = np.asarray(hp.pix2vec(nside, np.arange(npix), nest=nest))
cov = np.dot(xyz * mass, xyz.T)
L, V = np.linalg.eigh(cov)
if np.linalg.det(V) < 0:
V = -V
return V
HEALPix array to a Mollweide projection and applying scipy.ndimage.label
Assumes non-nested HEALPix map.
Parameters:
pixels : Pixel values associated to (sparse) HEALPix array
values : (Sparse) HEALPix array of data values
nside : HEALPix dimensionality
threshold : Threshold value for object detection
xsize : Size of Mollweide projection
Returns:
labels, nlabels
"""
proj = healpy.projector.MollweideProj(xsize=xsize)
vec = healpy.pix2vec(nside,pixels)
xy = proj.vec2xy(vec)
ij = proj.xy2ij(xy)
xx,yy = proj.ij2xy()
# Convert to Mollweide
searchims = []
if values.ndim < 2: iterate = [values]
else: iterate = values.T
for i,value in enumerate(iterate):
logger.debug("Labeling slice %i...")
searchim = numpy.zeros(xx.shape,dtype=bool)
select = (value > threshold)
yidx = ij[0][select]; xidx = ij[1][select]
searchim[yidx,xidx] |= True
searchims.append( searchim )
searchims = numpy.array(searchims)
def match_hpx_pixel(nside, nest, nside_pix, ipix_ring):
"""
"""
import healpy as hp
ipix_in = np.arange(12 * nside * nside)
vecs = hp.pix2vec(nside, ipix_in, nest)
pix_match = hp.vec2pix(nside_pix, vecs[0], vecs[1], vecs[2]) == ipix_ring
return ipix_in[pix_match]
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.hpc.set_map(to_observe)
best_fit_shifts = self.hpc.minimize()
Draw vectors from a uniform distribution within a HEALpixel.
nside : healpix nside parameter
ipix : pixel number(s)
"""
if not(nest):
ipix = healpy.ring2nest(nside, ipix=ipix)
norder = nside2norder(nside)
nUp = 29 - norder
iUp = ipix * 4**nUp
if np.iterable(ipix):
iUp += np.random.randint(0, 4**nUp, size=np.size(ipix))
else:
iUp += np.random.randint(0, 4**nUp)
vec = healpy.pix2vec(nside=2**29, ipix=iUp, nest=True)
return vec