Skip to content

Commit

Permalink
Merge b4cde16 into fcb4089
Browse files Browse the repository at this point in the history
  • Loading branch information
@lpsinger
lpsinger committed Apr 19, 2018
2 parents fcb4089 + b4cde16 commit 9a270d2
Show file tree
Hide file tree
Showing 3 changed files with 85 additions and 61 deletions.
39 changes: 24 additions & 15 deletions astropy_healpix/core.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -267,6 +267,7 @@ def healpix_to_lonlat(healpix_index, nside, dx=None, dy=None, order='ring'):
raise ValueError('Either both or neither dx and dy must be specified')

healpix_index = np.asarray(healpix_index, dtype=np.int64)
nside = np.asarray(nside, dtype=np.int64)

if dx is None and dy is not None:
dx = 0.5
Expand All @@ -282,14 +283,16 @@ def healpix_to_lonlat(healpix_index, nside, dx=None, dy=None, order='ring'):
dx = dx.ravel()
dy = dy.ravel()

shape = healpix_index.shape
healpix_index = healpix_index.ravel()
nside = int(nside)

_validate_healpix_index('healpix_index', healpix_index, nside)
_validate_nside(nside)
order = _validate_order(order)

healpix_index, nside = np.broadcast_arrays(healpix_index, nside)

shape = healpix_index.shape
healpix_index = healpix_index.ravel()
nside = nside.ravel()

if dx is None:
lon, lat = core_cython.healpix_to_lonlat(healpix_index, nside, order)
else:
Expand Down Expand Up @@ -331,14 +334,14 @@ def lonlat_to_healpix(lon, lat, nside, return_offsets=False, order='ring'):
lon = lon.to(u.rad).value
lat = lat.to(u.rad).value

lon, lat = np.broadcast_arrays(lon, lat)
lon, lat, nside = np.broadcast_arrays(lon, lat, nside)

shape = np.shape(lon)

lon = lon.astype(float).ravel()
lat = lat.astype(float).ravel()
nside = nside.astype(np.int64).ravel()

nside = int(nside)
_validate_nside(nside)
order = _validate_order(order)

Expand Down Expand Up @@ -368,9 +371,11 @@ def nested_to_ring(nested_index, nside):
"""

nested_index = np.asarray(nested_index, dtype=np.int64)
nside = np.asarray(nside, dtype=np.int64)
nested_index, nside = np.broadcast_arrays(nested_index, nside)
shape = nested_index.shape
nested_index = nested_index.ravel()
nside = int(nside)
nside = nside.ravel()

_validate_healpix_index('nested_index', nested_index, nside)
_validate_nside(nside)
Expand All @@ -397,9 +402,11 @@ def ring_to_nested(ring_index, nside):
"""

ring_index = np.asarray(ring_index, dtype=np.int64)
nside = np.asarray(nside, dtype=np.int64)
ring_index, nside = np.broadcast_arrays(ring_index, nside)
shape = ring_index.shape
ring_index = ring_index.ravel()
nside = int(nside)
nside = nside.ravel()

_validate_healpix_index('ring_index', ring_index, nside)
_validate_nside(nside)
Expand Down Expand Up @@ -436,13 +443,13 @@ def bilinear_interpolation_weights(lon, lat, nside, order='ring'):
lon = lon.to(u.rad).value
lat = lat.to(u.rad).value

lon, lat = np.broadcast_arrays(lon, lat)
lon, lat, nside = np.broadcast_arrays(lon, lat, nside)

shape = np.shape(lon)

lon = lon.astype(float).ravel()
lat = lat.astype(float).ravel()
nside = int(nside)
nside = nside.astype(np.int64).ravel()

order = _validate_order(order)
_validate_nside(nside)
Expand Down Expand Up @@ -509,16 +516,18 @@ def neighbours(healpix_index, nside, order='ring'):
"""

healpix_index = np.asarray(healpix_index, dtype=np.int64)
nside = int(nside)

shape = np.shape(healpix_index)

healpix_index = healpix_index.ravel()
nside = np.asarray(nside, dtype=np.int64)

_validate_healpix_index('healpix_index', healpix_index, nside)
_validate_nside(nside)
order = _validate_order(order)

healpix_index, nside = np.broadcast_arrays(healpix_index, nside)
shape = np.shape(healpix_index)

healpix_index = healpix_index.ravel()
nside = nside.ravel()

neigh = core_cython.neighbours(healpix_index, nside, order)
return _restore_shape(neigh, shape=(8,) + shape)

Expand Down
91 changes: 49 additions & 42 deletions astropy_healpix/core_cython.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,8 @@ def _validate_order(str order):

@cython.boundscheck(False)
def healpix_to_lonlat(np.ndarray[int64_t, ndim=1, mode="c"] healpix_index,
int nside, str order):
np.ndarray[int64_t, ndim=1, mode="c"] nside,
str order):
"""
Convert HEALPix indices to longitudes/latitudes.
Expand All @@ -46,7 +47,7 @@ def healpix_to_lonlat(np.ndarray[int64_t, ndim=1, mode="c"] healpix_index,
----------
healpix_index : `~numpy.ndarray`
1-D array of HEALPix indices
nside : int
nside : `~numpy.ndarray`
Number of pixels along the side of each of the 12 top-level HEALPix tiles
order : { 'nested' | 'ring' }
Order of HEALPix pixels
Expand All @@ -71,12 +72,12 @@ def healpix_to_lonlat(np.ndarray[int64_t, ndim=1, mode="c"] healpix_index,

if order == 'nested':
for i in prange(n, nogil=True, schedule='static'):
xy_index = healpixl_nested_to_xy(healpix_index[i], nside)
healpixl_to_radec(xy_index, nside, dx, dy, &lon[i], &lat[i])
xy_index = healpixl_nested_to_xy(healpix_index[i], <int>nside[i])
healpixl_to_radec(xy_index, <int>nside[i], dx, dy, &lon[i], &lat[i])
elif order == 'ring':
for i in prange(n, nogil=True, schedule='static'):
xy_index = healpixl_ring_to_xy(healpix_index[i], nside)
healpixl_to_radec(xy_index, nside, dx, dy, &lon[i], &lat[i])
xy_index = healpixl_ring_to_xy(healpix_index[i], <int>nside[i])
healpixl_to_radec(xy_index, <int>nside[i], dx, dy, &lon[i], &lat[i])

return lon, lat

Expand All @@ -85,7 +86,8 @@ def healpix_to_lonlat(np.ndarray[int64_t, ndim=1, mode="c"] healpix_index,
def healpix_with_offset_to_lonlat(np.ndarray[int64_t, ndim=1, mode="c"] healpix_index,
np.ndarray[double_t, ndim=1, mode="c"] dx,
np.ndarray[double_t, ndim=1, mode="c"] dy,
int nside, str order):
np.ndarray[int64_t, ndim=1, mode="c"] nside,
str order):
"""
Convert HEALPix indices to longitudes/latitudes
Expand All @@ -100,7 +102,7 @@ def healpix_with_offset_to_lonlat(np.ndarray[int64_t, ndim=1, mode="c"] healpix_
dx, dy : `~numpy.ndarray`
1-D arrays of offsets inside the HEALPix pixel, which must be in the
range [0:1] (0.5 is the center of the HEALPix pixels)
nside : int
nside : `~numpy.ndarray`
Number of pixels along the side of each of the 12 top-level HEALPix tiles
order : { 'nested' | 'ring' }
Order of HEALPix pixels
Expand All @@ -122,20 +124,21 @@ def healpix_with_offset_to_lonlat(np.ndarray[int64_t, ndim=1, mode="c"] healpix_

if order == 'nested':
for i in prange(n, nogil=True, schedule='static'):
xy_index = healpixl_nested_to_xy(healpix_index[i], nside)
healpixl_to_radec(xy_index, nside, dx[i], dy[i], &lon[i], &lat[i])
xy_index = healpixl_nested_to_xy(healpix_index[i], <int>nside[i])
healpixl_to_radec(xy_index, <int>nside[i], dx[i], dy[i], &lon[i], &lat[i])
elif order == 'ring':
for i in prange(n, nogil=True, schedule='static'):
xy_index = healpixl_ring_to_xy(healpix_index[i], nside)
healpixl_to_radec(xy_index, nside, dx[i], dy[i], &lon[i], &lat[i])
xy_index = healpixl_ring_to_xy(healpix_index[i], <int>nside[i])
healpixl_to_radec(xy_index, <int>nside[i], dx[i], dy[i], &lon[i], &lat[i])

return lon, lat


@cython.boundscheck(False)
def lonlat_to_healpix(np.ndarray[double_t, ndim=1, mode="c"] lon,
np.ndarray[double_t, ndim=1, mode="c"] lat,
int nside, str order):
np.ndarray[int64_t, ndim=1, mode="c"] nside,
str order):
"""
Convert longitudes/latitudes to HEALPix indices
Expand All @@ -146,7 +149,7 @@ def lonlat_to_healpix(np.ndarray[double_t, ndim=1, mode="c"] lon,
----------
lon, lat : `~numpy.ndarray`
1-D arrays of longitude and latitude in radians
nside : int
nside : `~numpy.ndarray`
Number of pixels along the side of each of the 12 top-level HEALPix tiles
order : { 'nested' | 'ring' }
Order of HEALPix pixels
Expand All @@ -168,20 +171,21 @@ def lonlat_to_healpix(np.ndarray[double_t, ndim=1, mode="c"] lon,

if order == 'nested':
for i in prange(n, nogil=True, schedule='static'):
xy_index = radec_to_healpixlf(lon[i], lat[i], nside, &dx, &dy)
healpix_index[i] = healpixl_xy_to_nested(xy_index, nside)
xy_index = radec_to_healpixlf(lon[i], lat[i], <int>nside[i], &dx, &dy)
healpix_index[i] = healpixl_xy_to_nested(xy_index, <int>nside[i])
elif order == 'ring':
for i in prange(n, nogil=True, schedule='static'):
xy_index = radec_to_healpixlf(lon[i], lat[i], nside, &dx, &dy)
healpix_index[i] = healpixl_xy_to_ring(xy_index, nside)
xy_index = radec_to_healpixlf(lon[i], lat[i], <int>nside[i], &dx, &dy)
healpix_index[i] = healpixl_xy_to_ring(xy_index, <int>nside[i])

return healpix_index


@cython.boundscheck(False)
def lonlat_to_healpix_with_offset(np.ndarray[double_t, ndim=1, mode="c"] lon,
np.ndarray[double_t, ndim=1, mode="c"] lat,
int nside, str order):
np.ndarray[int64_t, ndim=1, mode="c"] nside,
str order):
"""
Convert longitudes/latitudes to healpix indices
Expand All @@ -192,7 +196,7 @@ def lonlat_to_healpix_with_offset(np.ndarray[double_t, ndim=1, mode="c"] lon,
----------
lon, lat : `~numpy.ndarray`
1-D arrays of longitude and latitude in radians
nside : int
nside : `~numpy.ndarray`
Number of pixels along the side of each of the 12 top-level HEALPix tiles
order : { 'nested' | 'ring' }
Order of HEALPix pixels
Expand All @@ -217,26 +221,27 @@ def lonlat_to_healpix_with_offset(np.ndarray[double_t, ndim=1, mode="c"] lon,

if order == 'nested':
for i in prange(n, nogil=True, schedule='static'):
xy_index = radec_to_healpixlf(lon[i], lat[i], nside, &dx[i], &dy[i])
healpix_index[i] = healpixl_xy_to_nested(xy_index, nside)
xy_index = radec_to_healpixlf(lon[i], lat[i], <int>nside[i], &dx[i], &dy[i])
healpix_index[i] = healpixl_xy_to_nested(xy_index, <int>nside[i])
elif order == 'ring':
for i in prange(n, nogil=True, schedule='static'):
xy_index = radec_to_healpixlf(lon[i], lat[i], nside, &dx[i], &dy[i])
healpix_index[i] = healpixl_xy_to_ring(xy_index, nside)
xy_index = radec_to_healpixlf(lon[i], lat[i], <int>nside[i], &dx[i], &dy[i])
healpix_index[i] = healpixl_xy_to_ring(xy_index, <int>nside[i])

return healpix_index, dx, dy


@cython.boundscheck(False)
def nested_to_ring(np.ndarray[int64_t, ndim=1, mode="c"] nested_index, int nside):
def nested_to_ring(np.ndarray[int64_t, ndim=1, mode="c"] nested_index,
np.ndarray[int64_t, ndim=1, mode="c"] nside):
"""
Convert a HEALPix 'nested' index to a HEALPix 'ring' index
Parameters
----------
nested_index : `~numpy.ndarray`
Healpix index using the 'nested' ordering
nside : int
nside : `~numpy.ndarray`
Number of pixels along the side of each of the 12 top-level HEALPix tiles
Returns
Expand All @@ -250,21 +255,22 @@ def nested_to_ring(np.ndarray[int64_t, ndim=1, mode="c"] nested_index, int nside
cdef np.ndarray[int64_t, ndim=1, mode="c"] ring_index = np.zeros(n, dtype=npy_int64)

for i in prange(n, nogil=True, schedule='static'):
ring_index[i] = healpixl_xy_to_ring(healpixl_nested_to_xy(nested_index[i], nside), nside)
ring_index[i] = healpixl_xy_to_ring(healpixl_nested_to_xy(nested_index[i], <int>nside[i]), <int>nside[i])

return ring_index


@cython.boundscheck(False)
def ring_to_nested(np.ndarray[int64_t, ndim=1, mode="c"] ring_index, int nside):
def ring_to_nested(np.ndarray[int64_t, ndim=1, mode="c"] ring_index,
np.ndarray[int64_t, ndim=1, mode="c"] nside):
"""
Convert a HEALPix 'ring' index to a HEALPix 'nested' index
Parameters
----------
ring_index : `~numpy.ndarray`
Healpix index using the 'ring' ordering
nside : int
nside : `~numpy.ndarray`
Number of pixels along the side of each of the 12 top-level HEALPix tiles
Returns
Expand All @@ -278,7 +284,7 @@ def ring_to_nested(np.ndarray[int64_t, ndim=1, mode="c"] ring_index, int nside):
cdef np.ndarray[int64_t, ndim=1, mode="c"] nested_index = np.zeros(n, dtype=npy_int64)

for i in prange(n, nogil=True, schedule='static'):
nested_index[i] = healpixl_xy_to_nested(healpixl_ring_to_xy(ring_index[i], nside), nside)
nested_index[i] = healpixl_xy_to_nested(healpixl_ring_to_xy(ring_index[i], <int>nside[i]), <int>nside[i])

return nested_index

Expand All @@ -287,7 +293,8 @@ def ring_to_nested(np.ndarray[int64_t, ndim=1, mode="c"] ring_index, int nside):
@cython.boundscheck(False)
def bilinear_interpolation_weights(np.ndarray[double_t, ndim=1, mode="c"] lon,
np.ndarray[double_t, ndim=1, mode="c"] lat,
int nside, str order):
np.ndarray[int64_t, ndim=1, mode="c"] nside,
str order):
"""
Get the four neighbours for each (lon, lat) position and the weight
associated with each one for bilinear interpolation.
Expand All @@ -296,7 +303,7 @@ def bilinear_interpolation_weights(np.ndarray[double_t, ndim=1, mode="c"] lon,
----------
lon, lat : `~numpy.ndarray`
1-D arrays of longitude and latitude in radians
nside : int
nside : `~numpy.ndarray`
Number of pixels along the side of each of the 12 top-level HEALPix tiles
order : { 'nested' | 'ring' }
Order of HEALPix pixels
Expand Down Expand Up @@ -345,10 +352,10 @@ def bilinear_interpolation_weights(np.ndarray[double_t, ndim=1, mode="c"] lon,
abort()

for i in range(n):
interpolate_weights(lon[i], lat[i], indices_indiv, weights_indiv, nside)
interpolate_weights(lon[i], lat[i], indices_indiv, weights_indiv, <int>nside[i])
for j in range(4):
if order_int == 0:
indices[j, i] = healpixl_xy_to_nested(healpixl_ring_to_xy(indices_indiv[j], nside), nside)
indices[j, i] = healpixl_xy_to_nested(healpixl_ring_to_xy(indices_indiv[j], <int>nside[i]), <int>nside[i])
else:
indices[j, i] = indices_indiv[j]
weights[j, i] = weights_indiv[j]
Expand All @@ -358,15 +365,15 @@ def bilinear_interpolation_weights(np.ndarray[double_t, ndim=1, mode="c"] lon,

@cython.boundscheck(False)
def neighbours(np.ndarray[int64_t, ndim=1, mode="c"] healpix_index,
int nside, str order):
np.ndarray[int64_t, ndim=1, mode="c"] nside, str order):
"""
Find all the HEALPix pixels that are the neighbours of a HEALPix pixel
Parameters
----------
healpix_index : `~numpy.ndarray`
1-D array of HEALPix indices
nside : int
nside : `~numpy.ndarray`
Number of pixels along the side of each of the 12 top-level HEALPix tiles
order : { 'nested' | 'ring' }
Order of HEALPix pixels
Expand Down Expand Up @@ -417,8 +424,8 @@ def neighbours(np.ndarray[int64_t, ndim=1, mode="c"] healpix_index,

for i in prange(n, schedule='static'):

xy_index = healpixl_nested_to_xy(healpix_index[i], nside)
healpixl_get_neighbours(xy_index, neighbours_indiv, nside)
xy_index = healpixl_nested_to_xy(healpix_index[i], <int>nside[i])
healpixl_get_neighbours(xy_index, neighbours_indiv, <int>nside[i])

for j in range(8):
k = 4 - j
Expand All @@ -427,15 +434,15 @@ def neighbours(np.ndarray[int64_t, ndim=1, mode="c"] healpix_index,
if neighbours_indiv[k] < 0:
neighbours[j, i] = -1
else:
neighbours[j, i] = healpixl_xy_to_nested(neighbours_indiv[k], nside)
neighbours[j, i] = healpixl_xy_to_nested(neighbours_indiv[k], <int>nside[i])

elif order_int == 1:

for i in prange(n, schedule='static'):

xy_index = healpixl_ring_to_xy(healpix_index[i], nside)
xy_index = healpixl_ring_to_xy(healpix_index[i], <int>nside[i])

healpixl_get_neighbours(xy_index, neighbours_indiv, nside)
healpixl_get_neighbours(xy_index, neighbours_indiv, <int>nside[i])

for j in range(8):
k = 4 - j
Expand All @@ -444,7 +451,7 @@ def neighbours(np.ndarray[int64_t, ndim=1, mode="c"] healpix_index,
if neighbours_indiv[k] < 0:
neighbours[j, i] = -1
else:
neighbours[j, i] = healpixl_xy_to_ring(neighbours_indiv[k], nside)
neighbours[j, i] = healpixl_xy_to_ring(neighbours_indiv[k], <int>nside[i])

free(neighbours_indiv)

Expand Down

0 comments on commit 9a270d2

Please sign in to comment.