|
1
|
|
|
import healpy as hp |
|
|
|
|
|
|
2
|
|
|
import numpy as np |
|
|
|
|
|
|
3
|
|
|
import six |
|
4
|
|
|
|
|
5
|
|
|
from scipy.ndimage import map_coordinates |
|
|
|
|
|
|
6
|
|
|
|
|
7
|
|
|
from astropy.coordinates import Galactic, ICRS |
|
|
|
|
|
|
8
|
|
|
from astropy import units as u |
|
|
|
|
|
|
9
|
|
|
from astropy.coordinates import UnitSphericalRepresentation |
|
|
|
|
|
|
10
|
|
|
from astropy.wcs.utils import wcs_to_celestial_frame |
|
|
|
|
|
|
11
|
|
|
|
|
12
|
|
|
from ..special_values import UNSEEN |
|
13
|
|
|
from ..interpolation import FastBilinearInterpolation |
|
14
|
|
|
|
|
15
|
|
|
|
|
16
|
|
|
ORDER = {} |
|
17
|
|
|
ORDER['nearest-neighbor'] = 0 |
|
18
|
|
|
ORDER['bilinear'] = 1 |
|
19
|
|
|
ORDER['biquadratic'] = 2 |
|
20
|
|
|
ORDER['bicubic'] = 3 |
|
21
|
|
|
|
|
22
|
|
|
|
|
23
|
|
|
COORDSYS = { |
|
24
|
|
|
'g': Galactic(), |
|
25
|
|
|
'c': ICRS(), |
|
26
|
|
|
'icrs': ICRS(), |
|
27
|
|
|
} |
|
28
|
|
|
|
|
29
|
|
|
|
|
30
|
|
|
def _parse_coord_system(system): |
|
31
|
|
|
|
|
32
|
|
|
try: |
|
33
|
|
|
|
|
34
|
|
|
return COORDSYS[system.lower()] |
|
35
|
|
|
|
|
36
|
|
|
except KeyError: # pragma: no cover |
|
37
|
|
|
|
|
38
|
|
|
raise ValueError("Coordinate system %s is not known" % system) |
|
39
|
|
|
|
|
40
|
|
|
|
|
41
|
|
|
def _convert_world_coordinates(lon_in, lat_in, wcs_in, wcs_out): |
|
42
|
|
|
|
|
43
|
|
|
frame_in, lon_in_unit, lat_in_unit = wcs_in |
|
44
|
|
|
|
|
45
|
|
|
wcs_out = wcs_out.celestial |
|
46
|
|
|
frame_out = wcs_to_celestial_frame(wcs_out) |
|
47
|
|
|
lon_out_unit = u.Unit(wcs_out.wcs.cunit[0]) |
|
48
|
|
|
lat_out_unit = u.Unit(wcs_out.wcs.cunit[1]) |
|
49
|
|
|
|
|
50
|
|
|
data = UnitSphericalRepresentation(lon_in * lon_in_unit, |
|
51
|
|
|
lat_in * lat_in_unit) |
|
52
|
|
|
|
|
53
|
|
|
coords_in = frame_in.realize_frame(data) |
|
54
|
|
|
coords_out = coords_in.transform_to(frame_out) |
|
55
|
|
|
|
|
56
|
|
|
lon_out = coords_out.represent_as('unitspherical').lon.to(lon_out_unit).value |
|
57
|
|
|
lat_out = coords_out.represent_as('unitspherical').lat.to(lat_out_unit).value |
|
58
|
|
|
|
|
59
|
|
|
return lon_out, lat_out |
|
60
|
|
|
|
|
61
|
|
|
|
|
62
|
|
|
class FlatSkyToHealpixTransform(object): |
|
63
|
|
|
""" |
|
64
|
|
|
A class to perform transformation from a flat sky projection to Healpix optimized to be used for the same |
|
|
|
|
|
|
65
|
|
|
transformation over and over again. |
|
66
|
|
|
|
|
67
|
|
|
The constructor will pre-compute all needed quantities for the transformation, and the __call__ method just applies |
|
|
|
|
|
|
68
|
|
|
the transformation. This avoids to re-compute the same quantities over and over again. |
|
69
|
|
|
""" |
|
70
|
|
|
|
|
71
|
|
|
def __init__(self, wcs_in, coord_system_out, nside, pixels_id, input_shape, order='bilinear', nested=False): |
|
|
|
|
|
|
72
|
|
|
|
|
73
|
|
|
# Look up lon, lat of pixels in output system and convert colatitude theta |
|
74
|
|
|
# and longitude phi to longitude and latitude. |
|
75
|
|
|
theta, phi = hp.pix2ang(nside, pixels_id, nested) |
|
76
|
|
|
|
|
77
|
|
|
lon_out = np.degrees(phi) |
|
78
|
|
|
lat_out = 90. - np.degrees(theta) |
|
79
|
|
|
|
|
80
|
|
|
# Convert between celestial coordinates |
|
81
|
|
|
coord_system_out = _parse_coord_system(coord_system_out) |
|
82
|
|
|
|
|
83
|
|
|
with np.errstate(invalid='ignore'): |
|
84
|
|
|
lon_in, lat_in = _convert_world_coordinates(lon_out, lat_out, (coord_system_out, u.deg, u.deg), wcs_in) |
|
|
|
|
|
|
85
|
|
|
|
|
86
|
|
|
# Look up pixels in input system |
|
87
|
|
|
yinds, xinds = wcs_in.wcs_world2pix(lon_in, lat_in, 0) |
|
88
|
|
|
|
|
89
|
|
|
self._coords = [xinds, yinds] |
|
90
|
|
|
|
|
91
|
|
|
# Interpolate |
|
92
|
|
|
|
|
93
|
|
|
if isinstance(order, six.string_types): |
|
94
|
|
|
order = ORDER[order] |
|
95
|
|
|
|
|
96
|
|
|
self._order = order |
|
97
|
|
|
|
|
98
|
|
|
self._interpolator = FastBilinearInterpolation(input_shape, self._coords) |
|
99
|
|
|
|
|
100
|
|
|
def __call__(self, data, fill_value=UNSEEN): |
|
101
|
|
|
|
|
102
|
|
|
# healpix_data = map_coordinates(data, self._coords, |
|
103
|
|
|
# order=self._order, |
|
104
|
|
|
# mode='constant', cval=fill_value) |
|
105
|
|
|
|
|
106
|
|
|
healpix_data = self._interpolator(data) |
|
107
|
|
|
|
|
108
|
|
|
return healpix_data |
|
109
|
|
|
|
giacomov /
hawc_hal
Loading history...
The coding style of this project requires that you add a docstring to this code element. Below, you find an example for methods:
If you would like to know more about docstrings, we recommend to read PEP-257: Docstring Conventions.