|
1
|
|
|
# Licensed under a 3-clause BSD style license - see LICENSE.rst |
|
2
|
|
|
from __future__ import absolute_import, division, print_function, unicode_literals |
|
3
|
|
|
import numpy as np |
|
4
|
|
|
import astropy.units as u |
|
5
|
|
|
from scipy.ndimage import map_coordinates |
|
|
|
|
|
|
6
|
|
|
|
|
7
|
|
|
|
|
8
|
|
|
__all__ = ["reproject_car_to_hpx", "reproject_car_to_wcs"] |
|
9
|
|
|
|
|
10
|
|
|
|
|
11
|
|
|
def _get_input_pix_celestial(wcs_in, wcs_out, shape_out): |
|
12
|
|
|
""" |
|
13
|
|
|
Get the pixel coordinates of the pixels in an array of shape ``shape_out`` |
|
14
|
|
|
in the input WCS. |
|
15
|
|
|
""" |
|
16
|
|
|
from reproject.wcs_utils import convert_world_coordinates |
|
|
|
|
|
|
17
|
|
|
|
|
18
|
|
|
# TODO: for now assuming that coordinates are spherical, not |
|
|
|
|
|
|
19
|
|
|
# necessarily the case. Also assuming something about the order of the |
|
20
|
|
|
# arguments. |
|
21
|
|
|
|
|
22
|
|
|
# Generate pixel coordinates of output image |
|
23
|
|
|
xp_out, yp_out = np.indices(shape_out, dtype=float)[::-1] |
|
24
|
|
|
|
|
25
|
|
|
# Convert output pixel coordinates to pixel coordinates in original image |
|
26
|
|
|
# (using pixel centers). |
|
27
|
|
|
xw_out, yw_out = wcs_out.wcs_pix2world(xp_out, yp_out, 0) |
|
28
|
|
|
|
|
29
|
|
|
xw_in, yw_in = convert_world_coordinates(xw_out, yw_out, wcs_out, wcs_in) |
|
30
|
|
|
|
|
31
|
|
|
xp_in, yp_in = wcs_in.wcs_world2pix(xw_in, yw_in, 0) |
|
32
|
|
|
|
|
33
|
|
|
return xp_in, yp_in |
|
34
|
|
|
|
|
35
|
|
|
|
|
36
|
|
|
def reproject_car_to_hpx(input_data, coord_system_out, nside, order=1, nested=False): |
|
|
|
|
|
|
37
|
|
|
import healpy as hp |
|
|
|
|
|
|
38
|
|
|
from reproject.wcs_utils import convert_world_coordinates |
|
|
|
|
|
|
39
|
|
|
from reproject.healpix.utils import parse_coord_system |
|
|
|
|
|
|
40
|
|
|
|
|
41
|
|
|
data, wcs_in = input_data |
|
42
|
|
|
|
|
43
|
|
|
npix = hp.nside2npix(nside) |
|
44
|
|
|
|
|
45
|
|
|
theta, phi = hp.pix2ang(nside, np.arange(npix), nested) |
|
46
|
|
|
lon_out = np.degrees(phi) |
|
47
|
|
|
lat_out = 90.0 - np.degrees(theta) |
|
48
|
|
|
|
|
49
|
|
|
# Convert between celestial coordinates |
|
50
|
|
|
coord_system_out = parse_coord_system(coord_system_out) |
|
51
|
|
|
with np.errstate(invalid="ignore"): |
|
52
|
|
|
lon_in, lat_in = convert_world_coordinates( |
|
53
|
|
|
lon_out, lat_out, (coord_system_out, u.deg, u.deg), wcs_in |
|
|
|
|
|
|
54
|
|
|
) |
|
55
|
|
|
|
|
56
|
|
|
# Look up pixels in input system |
|
57
|
|
|
yinds, xinds = wcs_in.wcs_world2pix(lon_in, lat_in, 0) |
|
58
|
|
|
|
|
59
|
|
|
# Interpolate |
|
60
|
|
|
data = np.pad(data, 3, mode="wrap") |
|
61
|
|
|
|
|
62
|
|
|
healpix_data = map_coordinates( |
|
63
|
|
|
data, [xinds + 3, yinds + 3], order=order, mode="wrap", cval=np.nan |
|
64
|
|
|
) |
|
65
|
|
|
|
|
66
|
|
|
return healpix_data, (~np.isnan(healpix_data)).astype(float) |
|
67
|
|
|
|
|
68
|
|
|
|
|
69
|
|
|
def reproject_car_to_wcs(input_data, wcs_out, shape_out, order=1): |
|
|
|
|
|
|
70
|
|
|
"""Reproject an all-sky CAR projection to another WCS projection. |
|
71
|
|
|
|
|
72
|
|
|
This method performs special handling of the projection edges to |
|
73
|
|
|
ensure that the interpolation of the CAR projection is correctly |
|
74
|
|
|
wrapped in longitude. |
|
75
|
|
|
""" |
|
76
|
|
|
slice_in, wcs_in = input_data |
|
77
|
|
|
|
|
78
|
|
|
array_new = np.zeros(shape_out) |
|
79
|
|
|
slice_out = array_new |
|
80
|
|
|
|
|
81
|
|
|
xp_in, yp_in = _get_input_pix_celestial( |
|
82
|
|
|
wcs_in.celestial, wcs_out.celestial, slice_out.shape |
|
83
|
|
|
) |
|
84
|
|
|
coordinates = np.array([yp_in.ravel(), xp_in.ravel()]) |
|
85
|
|
|
|
|
86
|
|
|
jmin, imin = np.floor(np.nanmin(coordinates, axis=1)).astype(int) - 1 |
|
87
|
|
|
jmax, imax = np.ceil(np.nanmax(coordinates, axis=1)).astype(int) + 1 |
|
88
|
|
|
|
|
89
|
|
|
ny, nx = slice_in.shape |
|
90
|
|
|
|
|
91
|
|
|
if imin >= nx or imax < 0 or jmin >= ny or jmax < 0: |
|
92
|
|
|
return array_new * np.nan, array_new.astype(float) |
|
93
|
|
|
|
|
94
|
|
|
# Pad by 3 pixels to ensure that cubic interpolation works |
|
95
|
|
|
slice_in = np.pad(slice_in, 3, mode="wrap") |
|
96
|
|
|
|
|
97
|
|
|
# Make sure image is floating point. We do this only now because |
|
98
|
|
|
# we want to avoid converting the whole input array if possible |
|
99
|
|
|
slice_in = np.asarray(slice_in, dtype=float) |
|
100
|
|
|
slice_out[:, :] = map_coordinates( |
|
101
|
|
|
slice_in, coordinates + 3, order=order, cval=np.nan, mode="constant" |
|
102
|
|
|
).reshape(slice_out.shape) |
|
103
|
|
|
|
|
104
|
|
|
return array_new, (~np.isnan(array_new)).astype(float) |
|
105
|
|
|
|
gammapy /
gammapy
Loading history...
