Skip to content

Commit

Permalink
Update healpy to upstream 1.10.3 release.
Browse files Browse the repository at this point in the history
  • Loading branch information
@tskisner
tskisner committed May 3, 2017
1 parent 340ae3d commit 815efae
Show file tree
Hide file tree
Showing 7 changed files with 191 additions and 41 deletions.
2 changes: 1 addition & 1 deletion README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ other work used to streamline that process.
Current basis of this source is:

Healpix_3.31_2016Aug26.tar.gz md5 = c0dc75e57f237b634fec97df55997918
healpy-1.10.1.tar.gz md5 = ee8750957a6fdfdafbef54e54788ec9b
healpy-1.10.3.tar.gz md5 = 6491777d1bcbd36d356551bf240d3a2f

Source changes from future releases will have to be merged by hand, so your
patience is appreciated.
Expand Down
2 changes: 1 addition & 1 deletion configure.ac
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ dnl
dnl +------------------------
dnl | Initialize package info
dnl +------------------------
AC_INIT([HEALPix], [3.31.1], [[email protected]], [healpix-autotools], [https://github.com/tskisner/healpix-autotools])
AC_INIT([HEALPix], [3.31.2], [[email protected]], [healpix-autotools], [https://github.com/tskisner/healpix-autotools])
AC_CONFIG_SRCDIR([Makefile.am])
AM_INIT_AUTOMAKE([foreign tar-ustar])
AC_CONFIG_HEADERS(config.h)
Expand Down
3 changes: 2 additions & 1 deletion src/healpy/healpy/fitsfunc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
"""Provides input and output functions for Healpix maps, alm, and cl.
"""
from __future__ import with_statement
from __future__ import division

import six
import gzip
Expand Down Expand Up @@ -208,7 +209,7 @@ def write_map(filename,m,nest=False,dtype=np.float32,fits_IDL=True,coord=None,pa
mm2 = np.asarray(mm)
cols.append(pf.Column(name=cn,
format='1024%s' % curr_fitsformat,
array=mm2.reshape(mm2.size/1024,1024),
array=mm2.reshape(mm2.size//1024,1024),
unit=cu))
else:
cols.append(pf.Column(name=cn,
Expand Down
112 changes: 99 additions & 13 deletions src/healpy/healpy/pixelfunc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,6 +127,40 @@ def check_theta_valid(theta):
if not((theta >= 0).all() and (theta <= np.pi + 1e-5).all()):
raise ValueError('THETA is out of range [0,pi]')

def lonlat2thetaphi(lon,lat):
""" Transform longitude and latitude (deg) into co-latitude and longitude (rad)
Parameters
----------
lon : int or array-like
Longitude in degrees
lat : int or array-like
Latitude in degrees
Returns
-------
theta, phi : float, scalar or array-like
The co-latitude and longitude in radians
"""
return np.pi/2. - np.radians(lat),np.radians(lon)

def thetaphi2lonlat(theta,phi):
""" Transform co-latitude and longitude (rad) into longitude and latitude (deg)
Parameters
----------
theta : int or array-like
Co-latitude in radians
phi : int or array-like
Longitude in radians
Returns
-------
lon, lat : float, scalar or array-like
The longitude and latitude in degrees
"""
return np.degrees(phi), 90. - np.degrees(theta)

def maptype(m):
"""Describe the type of the map (valid, single, sequence of maps).
Checks : the number of maps, that all maps have same length and that this
Expand Down Expand Up @@ -345,8 +379,8 @@ def ma(m, badval = UNSEEN, rtol = 1e-5, atol = 1e-8, copy = True):
else:
return tuple([ma(mm) for mm in m])

def ang2pix(nside,theta,phi,nest=False):
"""ang2pix : nside,theta[rad],phi[rad],nest=False -> ipix (default:RING)
def ang2pix(nside,theta,phi,nest=False,lonlat=False):
"""ang2pix : nside,theta[rad],phi[rad],nest=False,lonlat=False -> ipix (default:RING)
Parameters
----------
Expand All @@ -356,6 +390,9 @@ def ang2pix(nside,theta,phi,nest=False):
Angular coordinates of a point on the sphere
nest : bool, optional
if True, assume NESTED pixel ordering, otherwise, RING pixel ordering
lonlat : bool
If True, input angles are assumed to be longitude and latitude in degree,
otherwise, they are co-latitude and longitude in radians.
Returns
-------
Expand All @@ -381,16 +418,21 @@ def ang2pix(nside,theta,phi,nest=False):
>>> hp.ang2pix([1, 2, 4, 8, 16], np.pi/2, 0)
array([ 4, 12, 72, 336, 1440])
>>> hp.ang2pix([1, 2, 4, 8, 16], 0, 0, lonlat=True)
array([ 4, 12, 72, 336, 1440])
"""
if lonlat:
theta,phi = lonlat2thetaphi(theta,phi)
check_theta_valid(theta)
check_nside(nside)
if nest:
return pixlib._ang2pix_nest(nside,theta,phi)
else:
return pixlib._ang2pix_ring(nside,theta,phi)

def pix2ang(nside,ipix,nest=False):
"""pix2ang : nside,ipix,nest=False -> theta[rad],phi[rad] (default RING)
def pix2ang(nside,ipix,nest=False,lonlat=False):
"""pix2ang : nside,ipix,nest=False,lonlat=False -> theta[rad],phi[rad] (default RING)
Parameters
----------
Expand All @@ -400,6 +442,9 @@ def pix2ang(nside,ipix,nest=False):
Pixel indices
nest : bool, optional
if True, assume NESTED pixel ordering, otherwise, RING pixel ordering
lonlat : bool, optional
If True, return angles will be longitude and latitude in degree,
otherwise, angles will be longitude and co-latitude in radians (default)
Returns
-------
Expand All @@ -422,12 +467,20 @@ def pix2ang(nside,ipix,nest=False):
>>> hp.pix2ang([1, 2, 4, 8], 11)
(array([ 2.30052398, 0.84106867, 0.41113786, 0.2044802 ]), array([ 5.49778714, 5.89048623, 5.89048623, 5.89048623]))
>>> hp.pix2ang([1, 2, 4, 8], 11, lonlat=True)
(array([ 315. , 337.5, 337.5, 337.5]), array([-41.8103149 , 41.8103149 , 66.44353569, 78.28414761]))
"""
check_nside(nside)
if nest:
return pixlib._pix2ang_nest(nside, ipix)
theta,phi = pixlib._pix2ang_nest(nside, ipix)
else:
theta,phi = pixlib._pix2ang_ring(nside,ipix)

if lonlat:
return thetaphi2lonlat(theta,phi)
else:
return pixlib._pix2ang_ring(nside,ipix)
return theta, phi

def xyf2pix(nside,x,y,face,nest=False):
"""xyf2pix : nside,x,y,face,nest=False -> ipix (default:RING)
Expand Down Expand Up @@ -588,7 +641,7 @@ def pix2vec(nside,ipix,nest=False):
else:
return pixlib._pix2vec_ring(nside,ipix)

def ang2vec(theta, phi):
def ang2vec(theta, phi, lonlat=False):
"""ang2vec : convert angles to 3D position vector
Parameters
Expand All @@ -597,6 +650,9 @@ def ang2vec(theta, phi):
colatitude in radians measured southward from north pole (in [0,pi]).
phi : float, scalar or array-like
longitude in radians measured eastward (in [0, 2*pi]).
lonlat : bool
If True, input angles are assumed to be longitude and latitude in degree,
otherwise, they are co-latitude and longitude in radians.
Returns
-------
Expand All @@ -608,19 +664,24 @@ def ang2vec(theta, phi):
--------
vec2ang, rotator.dir2vec, rotator.vec2dir
"""
if lonlat:
theta,phi = lonlat2thetaphi(theta,phi)
check_theta_valid(theta)
sintheta = np.sin(theta)
return np.array([sintheta*np.cos(phi),
sintheta*np.sin(phi),
np.cos(theta)]).T

def vec2ang(vectors):
def vec2ang(vectors, lonlat=False):
"""vec2ang: vectors [x, y, z] -> theta[rad], phi[rad]
Parameters
----------
vectors : float, array-like
the vector(s) to convert, shape is (3,) or (N, 3)
lonlat : bool, optional
If True, return angles will be longitude and latitude in degree,
otherwise, angles will be longitude and co-latitude in radians (default)
Returns
-------
Expand All @@ -636,7 +697,10 @@ def vec2ang(vectors):
theta = np.arccos(vectors[:,2]/dnorm)
phi = np.arctan2(vectors[:,1],vectors[:,0])
phi[phi < 0] += 2*np.pi
return theta, phi
if lonlat:
return thetaphi2lonlat(theta,phi)
else:
return theta, phi

def ring2nest(nside, ipix):
"""Convert pixel number from RING ordering to NESTED ordering.
Expand Down Expand Up @@ -1122,7 +1186,7 @@ def isnpixok(npix):
nside = np.sqrt(npix/12.)
return isnsideok(nside)

def get_interp_val(m,theta,phi,nest=False):
def get_interp_val(m,theta,phi,nest=False,lonlat=False):
"""Return the bi-linear interpolation value of a map using 4 nearest neighbours.
Parameters
Expand All @@ -1133,6 +1197,9 @@ def get_interp_val(m,theta,phi,nest=False):
angular coordinates of point at which to interpolate the map
nest : bool
if True, the is assumed to be in NESTED ordering.
lonlat : bool
If True, input angles are assumed to be longitude and latitude in degree,
otherwise, they are co-latitude and longitude in radians.
Returns
-------
Expand All @@ -1153,11 +1220,16 @@ def get_interp_val(m,theta,phi,nest=False):
>>> hp.get_interp_val(np.arange(12.), np.pi/2, np.pi/2 + 2*np.pi)
5.0
>>> hp.get_interp_val(np.arange(12.), np.linspace(0, np.pi, 10), 0)
array([ 1.5 , 1.5 , 1.5 , 2.20618428, 3.40206143,
5.31546486, 7.94639458, 9.5 , 9.5 , 9.5 ])
>>> hp.get_interp_val(np.arange(12.), 0, np.linspace(90, -90, 10), lonlat=True)
array([ 1.5 , 1.5 , 1.5 , 2.20618428, 3.40206143,
5.31546486, 7.94639458, 9.5 , 9.5 , 9.5 ])
"""
m2=m.ravel()
nside=npix2nside(m2.size)
if lonlat:
theta,phi = lonlat2thetaphi(theta,phi)
if nest:
r=pixlib._get_interpol_nest(nside,theta,phi)
else:
Expand All @@ -1170,7 +1242,7 @@ def get_interp_val(m,theta,phi,nest=False):
def get_neighbours(nside, theta, phi=None, nest=False):
raise NameError("get_neighbours has been renamed to get_interp_weights")

def get_interp_weights(nside,theta,phi=None,nest=False):
def get_interp_weights(nside,theta,phi=None,nest=False,lonlat=False):
"""Return the 4 closest pixels on the two rings above and below the
location and corresponding weights.
Weights are provided for bilinear interpolation along latitude and longitude
Expand All @@ -1184,6 +1256,9 @@ def get_interp_weights(nside,theta,phi=None,nest=False):
otherwise theta[rad],phi[rad] are angular coordinates
nest : bool
if ``True``, NESTED ordering, otherwise RING ordering.
lonlat : bool
If True, input angles are assumed to be longitude and latitude in degree,
otherwise, they are co-latitude and longitude in radians.
Returns
-------
Expand All @@ -1204,6 +1279,9 @@ def get_interp_weights(nside,theta,phi=None,nest=False):
>>> hp.get_interp_weights(1, 0, 0)
(array([1, 2, 3, 0]), array([ 0.25, 0.25, 0.25, 0.25]))
>>> hp.get_interp_weights(1, 0, 90, lonlat=True)
(array([1, 2, 3, 0]), array([ 0.25, 0.25, 0.25, 0.25]))
>>> hp.get_interp_weights(1, [0, np.pi/2], 0)
(array([[ 1, 4],
[ 2, 5],
Expand All @@ -1216,6 +1294,8 @@ def get_interp_weights(nside,theta,phi=None,nest=False):
check_nside(nside)
if phi == None:
theta,phi = pix2ang(nside,theta,nest=nest)
elif lonlat:
theta,phi = lonlat2thetaphi(theta,phi)
if nest:
r=pixlib._get_interpol_nest(nside,theta,phi)
else:
Expand All @@ -1224,7 +1304,7 @@ def get_interp_weights(nside,theta,phi=None,nest=False):
w=np.array(r[4:8])
return (p,w)

def get_all_neighbours(nside, theta, phi=None, nest=False):
def get_all_neighbours(nside, theta, phi=None, nest=False, lonlat=False):
"""Return the 8 nearest pixels.
Parameters
Expand All @@ -1236,6 +1316,9 @@ def get_all_neighbours(nside, theta, phi=None, nest=False):
otherwise, theta[rad],phi[rad] are angular coordinates
nest : bool
if ``True``, pixel number will be NESTED ordering, otherwise RING ordering.
lonlat : bool
If True, input angles are assumed to be longitude and latitude in degree,
otherwise, they are co-latitude and longitude in radians.
Returns
-------
Expand All @@ -1257,10 +1340,13 @@ def get_all_neighbours(nside, theta, phi=None, nest=False):
>>> hp.get_all_neighbours(1, np.pi/2, np.pi/2)
array([ 8, 4, 0, -1, 1, 6, 9, -1])
>>> hp.get_all_neighbours(1, 90, 0, lonlat=True)
array([ 8, 4, 0, -1, 1, 6, 9, -1])
"""
check_nside(nside)
if phi is not None:
theta = ang2pix(nside,theta, phi,nest=nest)
theta = ang2pix(nside, theta, phi, nest=nest, lonlat=lonlat)
if nest:
r=pixlib._get_neighbors_nest(nside,theta)
else:
Expand Down
Loading

0 comments on commit 815efae

Please sign in to comment.