/usr/lib/python3/dist-packages/reproject/array_utils.py is in python3-reproject 0.3.1-4.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 | import numpy as np
from functools import reduce
def iterate_over_celestial_slices(array_in, array_out, wcs):
"""
Given two arrays with the same number of dimensions, iterate over the
celestial slices in each.
The index of the celestial axes should match between the two arrays and
will be taken from ``wcs``. The iterator returns views, so these can be
used to modify the original arrays.
Parameters
----------
array_in : `~numpy.ndarray`
The input array for the reprojection
array_out : `~numpy.ndarray`
The output array for the reprojection
wcs : `~astropy.wcs.WCS`
The WCS for the input array (only used to identify the celestial axes).
Returns
-------
slice_in : `~numpy.ndarray`
A view to a celestial slice in ``array_in``
slice_out : `~numpy.ndarray`
A view to the corresponding slice in ``array_out``
"""
# First put lng/lat as first two dimensions in WCS/last two in Numpy
if wcs.wcs.lng == 0 and wcs.wcs.lat == 1:
array_in_view = array_in
array_out_view = array_out
elif wcs.wcs.lng == 1 and wcs.wcs.lat == 0:
array_in_view = array_in.swapaxes(-1, -2)
array_out_view = array_out.swapaxes(-1, -2)
else:
array_in_view = array_in.swapaxes(-2, -1 - wcs.wcs.lat).swapaxes(-1, -1 - wcs.wcs.lng)
array_out_view = array_out.swapaxes(-2, -1 - wcs.wcs.lat).swapaxes(-1, -1 - wcs.wcs.lng)
# Flatten remaining dimensions to make it easier to loop over
from operator import mul
nx_in = array_in_view.shape[-1]
ny_in = array_in_view.shape[-2]
n_remaining_in = reduce(mul, array_in_view.shape, 1) // nx_in // ny_in
nx_out = array_out_view.shape[-1]
ny_out = array_out_view.shape[-2]
n_remaining_out = reduce(mul, array_out_view.shape, 1) // nx_out // ny_out
if n_remaining_in != n_remaining_out:
raise ValueError("Number of non-celestial elements should match")
array_in_view = array_in_view.reshape(n_remaining_in, ny_in, nx_in)
array_out_view = array_out_view.reshape(n_remaining_out, ny_out, nx_out)
for slice_index in range(n_remaining_in):
yield array_in_view[slice_index], array_out_view[slice_index]
def pad_edge_1(array):
try:
return np.pad(array, 1, mode='edge')
except: # numpy < 1.7 workaround pragma: no cover
new_array = np.zeros((array.shape[0] + 2,
array.shape[1] + 2),
dtype=array.dtype)
new_array[1:-1, 1:-1] = array
new_array[0, 1:-1] = new_array[1, 1:-1]
new_array[-1, 1:-1] = new_array[-2, 1:-1]
new_array[:, 0] = new_array[:, 1]
new_array[:, -1] = new_array[:, -2]
return new_array
def map_coordinates(image, coords, **kwargs):
# In the built-in scipy map_coordinates, the values are defined at the
# center of the pixels. This means that map_coordinates does not
# correctly treat pixels that are in the outer half of the outer pixels.
# We solve this by extending the array, updating the pixel coordinates,
# then getting rid of values that were sampled in the range -1 to -0.5
# and n to n - 0.5.
from scipy.ndimage import map_coordinates as scipy_map_coordinates
image = pad_edge_1(image)
values = scipy_map_coordinates(image, coords + 1, **kwargs)
reset = np.zeros(coords.shape[1], dtype=bool)
for i in range(coords.shape[0]):
reset |= (coords[i] < -0.5)
reset |= (coords[i] > image.shape[i] - 0.5)
values[reset] = kwargs.get('cval', 0.)
return values
|