1 | # Licensed under a 3-clause BSD style license - see LICENSE.rst |
||
0 ignored issues
–
show
coding-style
introduced
by
Loading history...
|
|||
2 | """Utilities for dealing with HEALPix projections and mappings.""" |
||
3 | from __future__ import absolute_import, division, print_function, unicode_literals |
||
4 | from collections import OrderedDict |
||
5 | import re |
||
6 | import copy |
||
7 | import numpy as np |
||
8 | from ..extern import six |
||
9 | from astropy.io import fits |
||
10 | from astropy.coordinates import SkyCoord |
||
11 | from astropy.units import Quantity |
||
12 | from ..utils.scripts import make_path |
||
13 | from .wcs import WcsGeom |
||
14 | from .geom import MapGeom, MapCoord, pix_tuple_to_idx |
||
15 | from .geom import coordsys_to_frame, skycoord_to_lonlat |
||
16 | from .geom import find_and_read_bands, make_axes |
||
17 | |||
18 | # Not sure if we should expose this in the docs or not: |
||
19 | # HPX_FITS_CONVENTIONS, HpxConv |
||
20 | __all__ = ["HpxGeom"] |
||
21 | |||
22 | # Approximation of the size of HEALPIX pixels (in degrees) for a particular order. |
||
23 | # Used to convert from HEALPIX to WCS-based projections. |
||
24 | HPX_ORDER_TO_PIXSIZE = np.array( |
||
25 | [32.0, 16.0, 8.0, 4.0, 2.0, 1.0, 0.50, 0.25, 0.1, 0.05, 0.025, 0.01, 0.005, 0.002] |
||
26 | ) |
||
27 | |||
28 | |||
29 | class HpxConv(object): |
||
30 | """Data structure to define how a HEALPIX map is stored to FITS.""" |
||
31 | |||
32 | def __init__(self, convname, **kwargs): |
||
33 | self.convname = convname |
||
34 | self.colstring = kwargs.get("colstring", "CHANNEL") |
||
35 | self.firstcol = kwargs.get("firstcol", 1) |
||
36 | self.hduname = kwargs.get("hduname", "SKYMAP") |
||
37 | self.bands_hdu = kwargs.get("bands_hdu", "EBOUNDS") |
||
38 | self.quantity_type = kwargs.get("quantity_type", "integral") |
||
39 | self.coordsys = kwargs.get("coordsys", "COORDSYS") |
||
40 | |||
41 | def colname(self, indx): |
||
42 | return "{}{}".format(self.colstring, indx) |
||
43 | |||
44 | @classmethod |
||
45 | def create(cls, convname="gadf"): |
||
46 | return copy.deepcopy(HPX_FITS_CONVENTIONS[convname]) |
||
47 | |||
48 | |||
49 | HPX_FITS_CONVENTIONS = {} |
||
50 | """Various conventions for storing HEALPIX maps in FITS files""" |
||
51 | HPX_FITS_CONVENTIONS[None] = HpxConv("gadf", bands_hdu="BANDS") |
||
52 | HPX_FITS_CONVENTIONS["gadf"] = HpxConv("gadf", bands_hdu="BANDS") |
||
53 | HPX_FITS_CONVENTIONS["fgst-ccube"] = HpxConv("fgst-ccube") |
||
54 | HPX_FITS_CONVENTIONS["fgst-ltcube"] = HpxConv( |
||
55 | "fgst-ltcube", colstring="COSBINS", hduname="EXPOSURE", bands_hdu="CTHETABOUNDS" |
||
56 | ) |
||
57 | HPX_FITS_CONVENTIONS["fgst-bexpcube"] = HpxConv( |
||
58 | "fgst-bexpcube", colstring="ENERGY", hduname="HPXEXPOSURES", bands_hdu="ENERGIES" |
||
59 | ) |
||
60 | HPX_FITS_CONVENTIONS["fgst-srcmap"] = HpxConv( |
||
61 | "fgst-srcmap", hduname=None, quantity_type="differential" |
||
62 | ) |
||
63 | HPX_FITS_CONVENTIONS["fgst-template"] = HpxConv( |
||
64 | "fgst-template", colstring="ENERGY", bands_hdu="ENERGIES" |
||
65 | ) |
||
66 | HPX_FITS_CONVENTIONS["fgst-srcmap-sparse"] = HpxConv( |
||
67 | "fgst-srcmap-sparse", colstring=None, hduname=None, quantity_type="differential" |
||
68 | ) |
||
69 | HPX_FITS_CONVENTIONS["galprop"] = HpxConv( |
||
70 | "galprop", |
||
71 | colstring="Bin", |
||
72 | hduname="SKYMAP2", |
||
73 | bands_hdu="ENERGIES", |
||
74 | quantity_type="differential", |
||
75 | coordsys="COORDTYPE", |
||
76 | ) |
||
77 | HPX_FITS_CONVENTIONS["galprop2"] = HpxConv( |
||
78 | "galprop", |
||
79 | colstring="Bin", |
||
80 | hduname="SKYMAP2", |
||
81 | bands_hdu="ENERGIES", |
||
82 | quantity_type="differential", |
||
83 | ) |
||
84 | |||
85 | |||
86 | def unravel_hpx_index(idx, npix): |
||
87 | """Convert flattened global map index to an index tuple. |
||
88 | |||
89 | Parameters |
||
90 | ---------- |
||
91 | idx : `~numpy.ndarray` |
||
92 | Flat index. |
||
93 | npix : `~numpy.ndarray` |
||
94 | Number of pixels in each band. |
||
95 | |||
96 | Returns |
||
97 | ------- |
||
98 | idx : tuple of `~numpy.ndarray` |
||
99 | Index array for each dimension of the map. |
||
100 | """ |
||
101 | if npix.size == 1: |
||
102 | return tuple([idx]) |
||
103 | |||
104 | dpix = np.zeros(npix.size, dtype="i") |
||
105 | dpix[1:] = np.cumsum(npix.flat[:-1]) |
||
106 | bidx = np.searchsorted(np.cumsum(npix.flat), idx + 1) |
||
107 | pix = idx - dpix[bidx] |
||
108 | return tuple([pix] + list(np.unravel_index(bidx, npix.shape))) |
||
109 | |||
110 | |||
111 | def ravel_hpx_index(idx, npix): |
||
112 | """Convert map index tuple to a flattened index. |
||
113 | |||
114 | Parameters |
||
115 | ---------- |
||
116 | idx : tuple of `~numpy.ndarray` |
||
117 | |||
118 | Returns |
||
119 | ------- |
||
120 | idx : `~numpy.ndarray` |
||
121 | """ |
||
122 | if len(idx) == 1: |
||
123 | return idx[0] |
||
124 | |||
125 | # TODO: raise exception for indices that are out of bounds |
||
126 | |||
127 | idx0 = idx[0] |
||
128 | idx1 = np.ravel_multi_index(idx[1:], npix.shape, mode="clip") |
||
129 | npix = np.concatenate((np.array([0]), npix.flat[:-1])) |
||
130 | |||
131 | return idx0 + np.cumsum(npix)[idx1] |
||
132 | |||
133 | |||
134 | def coords_to_vec(lon, lat): |
||
135 | """Converts longitude and latitude coordinates to a unit 3-vector. |
||
136 | |||
137 | Returns |
||
138 | ------- |
||
139 | array(3,n) with v_x[i],v_y[i],v_z[i] = directional cosines |
||
140 | """ |
||
141 | phi = np.radians(lon) |
||
142 | theta = (np.pi / 2) - np.radians(lat) |
||
143 | sin_t = np.sin(theta) |
||
144 | cos_t = np.cos(theta) |
||
145 | |||
146 | x = sin_t * np.cos(phi) |
||
147 | y = sin_t * np.sin(phi) |
||
148 | z = cos_t |
||
149 | |||
150 | # Stack them into the output array |
||
151 | out = np.vstack((x, y, z)).swapaxes(0, 1) |
||
152 | return out |
||
153 | |||
154 | |||
155 | def get_nside_from_pix_size(pixsz): |
||
156 | """Get the NSIDE that is closest to the given pixel size. |
||
157 | |||
158 | Parameters |
||
159 | ---------- |
||
160 | pix : `~numpy.ndarray` |
||
161 | Pixel size in degrees. |
||
162 | |||
163 | Returns |
||
164 | ------- |
||
165 | nside : `~numpy.ndarray` |
||
166 | NSIDE parameter. |
||
167 | """ |
||
168 | import healpy as hp |
||
169 | |||
170 | pixsz = np.array(pixsz, ndmin=1) |
||
171 | nside = 2 ** np.linspace(1, 14, 14, dtype=int) |
||
172 | nside_pixsz = np.degrees(hp.nside2resol(nside)) |
||
173 | return nside[np.argmin(np.abs(nside_pixsz - pixsz[..., None]), axis=-1)] |
||
174 | |||
175 | |||
176 | def get_pix_size_from_nside(nside): |
||
177 | """Estimate of the pixel size from the HEALPIX nside coordinate. |
||
178 | |||
179 | This just uses a lookup table to provide a nice round number |
||
180 | for each HEALPIX order. |
||
181 | """ |
||
182 | order = nside_to_order(nside) |
||
183 | if np.any(order < 0) or np.any(order > 13): |
||
184 | raise ValueError("HEALPIX order must be 0 to 13. Got: {!r}".format(order)) |
||
185 | |||
186 | return HPX_ORDER_TO_PIXSIZE[order] |
||
187 | |||
188 | |||
189 | def hpx_to_axes(h, npix): |
||
190 | """Generate a sequence of bin edge vectors corresponding to the axes of a HPX object. |
||
191 | """ |
||
192 | x = h.ebins |
||
193 | z = np.arange(npix[-1] + 1) |
||
194 | return x, z |
||
195 | |||
196 | |||
197 | def hpx_to_coords(h, shape): |
||
198 | """Generate an N x D list of pixel center coordinates. |
||
199 | |||
200 | ``N`` is the number of pixels and ``D`` is the dimensionality of the map. |
||
201 | """ |
||
202 | x, z = hpx_to_axes(h, shape) |
||
203 | |||
204 | x = np.sqrt(x[0:-1] * x[1:]) |
||
205 | z = z[:-1] + 0.5 |
||
206 | |||
207 | x = np.ravel(np.ones(shape) * x[:, np.newaxis]) |
||
208 | z = np.ravel(np.ones(shape) * z[np.newaxis, :]) |
||
209 | |||
210 | return np.vstack((x, z)) |
||
211 | |||
212 | |||
213 | def make_hpx_to_wcs_mapping_centers(hpx, wcs): |
||
214 | """Make the mapping data needed to from from HPX pixelization to a WCS-based array. |
||
215 | |||
216 | Parameters |
||
217 | ---------- |
||
218 | hpx : `~gammapy.maps.HpxGeom` |
||
219 | The HEALPIX geometry. |
||
220 | wcs : `~gammapy.maps.WcsGeom` |
||
221 | The WCS geometry. |
||
222 | |||
223 | Returns |
||
224 | ------- |
||
225 | ipixs : array(nx,ny) |
||
226 | HEALPIX pixel indices for each WCS pixel |
||
227 | -1 indicates the wcs pixel does not contain the center of a HEALPIX pixel |
||
228 | mult_val : array |
||
229 | (nx,ny) of 1. |
||
230 | npix : tuple |
||
231 | (nx,ny) with the shape of the WCS grid |
||
232 | """ |
||
233 | npix = (int(wcs.wcs.crpix[0] * 2), int(wcs.wcs.crpix[1] * 2)) |
||
234 | mult_val = np.ones(npix).T.flatten() |
||
235 | sky_crds = hpx.get_sky_coords() |
||
236 | pix_crds = wcs.wcs_world2pix(sky_crds, 0).astype(int) |
||
237 | ipixs = -1 * np.ones(npix, int).T.flatten() |
||
238 | pix_index = npix[1] * pix_crds[0:, 0] + pix_crds[0:, 1] |
||
239 | |||
240 | if hpx._ipix is None: |
||
241 | for ipix, pix_crd in enumerate(pix_index): |
||
242 | ipixs[pix_crd] = ipix |
||
243 | else: |
||
244 | for pix_crd, ipix in zip(pix_index, hpx._ipix): |
||
245 | ipixs[pix_crd] = ipix |
||
246 | |||
247 | ipixs = ipixs.reshape(npix).T.flatten() |
||
248 | |||
249 | return ipixs, mult_val, npix |
||
250 | |||
251 | |||
252 | def make_hpx_to_wcs_mapping(hpx, wcs): |
||
253 | """Make the pixel mapping from HPX- to a WCS-based geometry. |
||
254 | |||
255 | Parameters |
||
256 | ---------- |
||
257 | hpx : `~gammapy.maps.HpxGeom` |
||
258 | The HEALPIX geometry |
||
259 | wcs : `~gammapy.maps.WcsGeom` |
||
260 | The WCS geometry |
||
261 | |||
262 | Returns |
||
263 | ------- |
||
264 | ipix : `~numpy.ndarray` |
||
265 | array(nx,ny) of HEALPIX pixel indices for each wcs pixel |
||
266 | mult_val : `~numpy.ndarray` |
||
267 | array(nx,ny) of 1./number of WCS pixels pointing at each HEALPIX pixel |
||
268 | npix : tuple |
||
269 | tuple(nx,ny) with the shape of the WCS grid |
||
270 | """ |
||
271 | import healpy as hp |
||
272 | |||
273 | npix = wcs.npix |
||
274 | |||
275 | # FIXME: Calculation of WCS pixel centers should be moved into a |
||
276 | # method of WcsGeom |
||
277 | pix_crds = np.dstack(np.meshgrid(np.arange(npix[0]), np.arange(npix[1]))) |
||
278 | pix_crds = pix_crds.swapaxes(0, 1).reshape((-1, 2)) |
||
279 | sky_crds = wcs.wcs.wcs_pix2world(pix_crds, 0) |
||
280 | sky_crds *= np.radians(1.0) |
||
281 | sky_crds[0:, 1] = (np.pi / 2) - sky_crds[0:, 1] |
||
282 | |||
283 | mask = ~np.any(np.isnan(sky_crds), axis=1) |
||
284 | ipix = -1 * np.ones((len(hpx.nside), int(npix[0] * npix[1])), int) |
||
285 | m = mask[None, :] * np.ones_like(ipix, dtype=bool) |
||
286 | |||
287 | ipix[m] = hp.ang2pix( |
||
288 | hpx.nside[..., None], |
||
289 | sky_crds[:, 1][mask][None, ...], |
||
290 | sky_crds[:, 0][mask][None, ...], |
||
291 | hpx.nest, |
||
292 | ).flatten() |
||
293 | |||
294 | # Here we are counting the number of HEALPIX pixels each WCS pixel |
||
295 | # points to and getting a multiplicative factor that tells use how |
||
296 | # to split up the counts in each HEALPIX pixel (by dividing the |
||
297 | # corresponding WCS pixels by the number of associated HEALPIX |
||
298 | # pixels). |
||
299 | mult_val = np.ones_like(ipix, dtype=float) |
||
300 | for i, t in enumerate(ipix): |
||
301 | count = np.unique(t, return_counts=True) |
||
302 | idx = np.searchsorted(count[0], t) |
||
303 | mult_val[i, ...] = 1.0 / count[1][idx] |
||
304 | |||
305 | if hpx.nside.size == 1: |
||
306 | ipix = np.squeeze(ipix, axis=0) |
||
307 | mult_val = np.squeeze(mult_val, axis=0) |
||
308 | |||
309 | return ipix, mult_val, npix |
||
310 | |||
311 | |||
312 | def match_hpx_pix(nside, nest, nside_pix, ipix_ring): |
||
313 | """TODO |
||
314 | """ |
||
315 | import healpy as hp |
||
316 | |||
317 | ipix_in = np.arange(12 * nside * nside) |
||
318 | vecs = hp.pix2vec(nside, ipix_in, nest) |
||
319 | pix_match = hp.vec2pix(nside_pix, vecs[0], vecs[1], vecs[2]) == ipix_ring |
||
320 | return ipix_in[pix_match] |
||
321 | |||
322 | |||
323 | def parse_hpxregion(region): |
||
324 | """Parse the ``HPX_REG`` header keyword into a list of tokens. |
||
325 | """ |
||
326 | m = re.match(r"([A-Za-z\_]*?)\((.*?)\)", region) |
||
327 | |||
328 | if m is None: |
||
329 | raise ValueError("Failed to parse hpx region string: {!r}".format(region)) |
||
330 | |||
331 | if not m.group(1): |
||
332 | return re.split(",", m.group(2)) |
||
333 | else: |
||
334 | return [m.group(1)] + re.split(",", m.group(2)) |
||
335 | |||
336 | |||
337 | def get_hpxregion_dir(region, coordsys): |
||
338 | """Get the reference direction for a HEALPIX region string. |
||
339 | |||
340 | Parameters |
||
341 | ---------- |
||
342 | region : str |
||
343 | A string describing a HEALPIX region |
||
344 | coordsys : {'CEL', 'GAL'} |
||
345 | Coordinate system |
||
346 | """ |
||
347 | import healpy as hp |
||
348 | |||
349 | frame = coordsys_to_frame(coordsys) |
||
350 | |||
351 | if region is None: |
||
352 | return SkyCoord(0.0, 0.0, frame=frame, unit="deg") |
||
353 | |||
354 | tokens = parse_hpxregion(region) |
||
355 | if tokens[0] in ["DISK", "DISK_INC"]: |
||
356 | lon, lat = float(tokens[1]), float(tokens[2]) |
||
357 | return SkyCoord(lon, lat, frame=frame, unit="deg") |
||
358 | elif tokens[0] == "HPX_PIXEL": |
||
359 | nside_pix = int(tokens[2]) |
||
360 | ipix_pix = int(tokens[3]) |
||
361 | if tokens[1] == "NESTED": |
||
362 | nest_pix = True |
||
363 | elif tokens[1] == "RING": |
||
364 | nest_pix = False |
||
365 | else: |
||
366 | raise ValueError("Invalid ordering scheme: {!r}".format(tokens[1])) |
||
367 | theta, phi = hp.pix2ang(nside_pix, ipix_pix, nest_pix) |
||
368 | lat = np.degrees((np.pi / 2) - theta) |
||
369 | lon = np.degrees(phi) |
||
370 | return SkyCoord(lon, lat, frame=frame, unit="deg") |
||
371 | else: |
||
372 | raise ValueError("Invalid region type: {!r}".format(tokens[0])) |
||
373 | |||
374 | |||
375 | def get_hpxregion_size(region): |
||
376 | """Get the approximate size of region (in degrees) from a HEALPIX region string. |
||
377 | """ |
||
378 | tokens = parse_hpxregion(region) |
||
379 | if tokens[0] in {"DISK", "DISK_INC"}: |
||
380 | return float(tokens[3]) |
||
381 | elif tokens[0] == "HPX_PIXEL": |
||
382 | pix_size = get_pix_size_from_nside(int(tokens[2])) |
||
383 | return 2.0 * pix_size |
||
384 | else: |
||
385 | raise ValueError("Invalid region type: {!r}".format(tokens[0])) |
||
386 | |||
387 | |||
388 | def is_power2(n): |
||
389 | """Check if an integer is a power of 2.""" |
||
390 | return (n > 0) & ((n & (n - 1)) == 0) |
||
391 | |||
392 | |||
393 | def nside_to_order(nside): |
||
394 | """Compute the ORDER given NSIDE. |
||
395 | |||
396 | Returns -1 for NSIDE values that are not a power of 2. |
||
397 | """ |
||
398 | nside = np.array(nside, ndmin=1) |
||
399 | order = -1 * np.ones_like(nside) |
||
400 | m = is_power2(nside) |
||
401 | order[m] = np.log2(nside[m]).astype(int) |
||
402 | return order |
||
403 | |||
404 | |||
405 | def upix_to_pix(upix): |
||
406 | """Get the pixel index and nside from a unique pixel number.""" |
||
407 | nside = np.power(2, np.floor(np.log2(upix / 4)) / 2).astype(int) |
||
408 | pix = upix - 4 * np.power(nside, 2) |
||
409 | return pix, nside |
||
410 | |||
411 | |||
412 | def pix_to_upix(pix, nside): |
||
413 | """Compute the unique pixel number from the pixel number and nside.""" |
||
414 | return pix + 4 * np.power(nside, 2) |
||
415 | |||
416 | |||
417 | def get_superpixels(idx, nside_subpix, nside_superpix, nest=True): |
||
418 | """Compute the indices of superpixels that contain a subpixel. |
||
419 | |||
420 | Parameters |
||
421 | ---------- |
||
422 | idx : `~numpy.ndarray` |
||
423 | Array of HEALPix pixel indices for subpixels of NSIDE |
||
424 | ``nside_subpix``. |
||
425 | nside_subpix : int or `~numpy.ndarray` |
||
426 | NSIDE of subpixel. |
||
427 | nside_superpix : int or `~numpy.ndarray` |
||
428 | NSIDE of superpixel. |
||
429 | nest : bool |
||
430 | If True, assume NESTED pixel ordering, otherwise, RING pixel |
||
431 | ordering. |
||
432 | |||
433 | Returns |
||
434 | ------- |
||
435 | idx_super : `~numpy.ndarray` |
||
436 | Indices of HEALpix pixels of nside ``nside_superpix`` that |
||
437 | contain pixel indices ``idx`` of nside ``nside_subpix``. |
||
438 | """ |
||
439 | |||
440 | import healpy as hp |
||
441 | |||
442 | idx = np.array(idx) |
||
443 | nside_superpix = np.asarray(nside_superpix) |
||
444 | nside_subpix = np.asarray(nside_subpix) |
||
445 | |||
446 | if not nest: |
||
447 | idx = hp.ring2nest(nside_subpix, idx) |
||
448 | |||
449 | if np.any(~is_power2(nside_superpix)) or np.any(~is_power2(nside_subpix)): |
||
450 | raise ValueError("NSIDE must be a power of 2.") |
||
451 | |||
452 | ratio = np.array((nside_subpix // nside_superpix) ** 2, ndmin=1) |
||
453 | idx //= ratio |
||
454 | |||
455 | if not nest: |
||
456 | m = idx == -1 |
||
457 | idx[m] = 0 |
||
458 | idx = hp.nest2ring(nside_superpix, idx) |
||
459 | idx[m] = -1 |
||
460 | |||
461 | return idx |
||
462 | |||
463 | |||
464 | def get_subpixels(idx, nside_superpix, nside_subpix, nest=True): |
||
465 | """Compute the indices of subpixels contained within superpixels. |
||
466 | This function returns an output array with one additional |
||
467 | dimension of size N for subpixel indices where N is the maximum |
||
468 | number of subpixels for any pair of ``nside_superpix`` and |
||
469 | ``nside_subpix``. If the number of subpixels is less than N the |
||
470 | remaining subpixel indices will be set to -1. |
||
471 | |||
472 | Parameters |
||
473 | ---------- |
||
474 | idx : `~numpy.ndarray` |
||
475 | Array of HEALPix pixel indices for superpixels of NSIDE |
||
476 | ``nside_superpix``. |
||
477 | nside_superpix : int or `~numpy.ndarray` |
||
478 | NSIDE of superpixel. |
||
479 | nside_subpix : int or `~numpy.ndarray` |
||
480 | NSIDE of subpixel. |
||
481 | nest : bool |
||
482 | If True, assume NESTED pixel ordering, otherwise, RING pixel |
||
483 | ordering. |
||
484 | |||
485 | Returns |
||
486 | ------- |
||
487 | idx_sub : `~numpy.ndarray` |
||
488 | Indices of HEALpix pixels of nside ``nside_subpix`` contained |
||
489 | within pixel indices ``idx`` of nside ``nside_superpix``. |
||
490 | """ |
||
491 | import healpy as hp |
||
492 | |||
493 | if not nest: |
||
494 | idx = hp.ring2nest(nside_superpix, idx) |
||
495 | |||
496 | idx = np.asarray(idx) |
||
497 | nside_superpix = np.asarray(nside_superpix) |
||
498 | nside_subpix = np.asarray(nside_subpix) |
||
499 | |||
500 | if np.any(~is_power2(nside_superpix)) or np.any(~is_power2(nside_subpix)): |
||
501 | raise ValueError("NSIDE must be a power of 2.") |
||
502 | |||
503 | # number of subpixels in each superpixel |
||
504 | npix = np.array((nside_subpix // nside_superpix) ** 2, ndmin=1) |
||
505 | x = np.arange(np.max(npix), dtype=int) |
||
506 | idx = idx * npix |
||
507 | |||
508 | if not np.all(npix[0] == npix): |
||
509 | x = np.broadcast_to(x, idx.shape + x.shape) |
||
510 | idx = idx[..., None] + x |
||
511 | idx[x >= np.broadcast_to(npix[..., None], x.shape)] = -1 |
||
512 | else: |
||
513 | idx = idx[..., None] + x |
||
514 | |||
515 | if not nest: |
||
516 | m = idx == -1 |
||
517 | idx[m] = 0 |
||
518 | idx = hp.nest2ring(nside_subpix[..., None], idx) |
||
519 | idx[m] = -1 |
||
520 | |||
521 | return idx |
||
522 | |||
523 | |||
524 | class HpxGeom(MapGeom): |
||
525 | """Geometry class for HEALPIX maps. |
||
526 | |||
527 | This class performs mapping between partial-sky indices (pixel |
||
528 | number within a HEALPIX region) and all-sky indices (pixel number |
||
529 | within an all-sky HEALPIX map). Multi-band HEALPIX geometries use |
||
530 | a global indexing scheme that assigns a unique pixel number based |
||
531 | on the all-sky index and band index. In the single-band case the |
||
532 | global index is the same as the HEALPIX index. |
||
533 | |||
534 | By default the constructor will return an all-sky map. |
||
535 | Partial-sky maps can be defined with the ``region`` argument. |
||
536 | |||
537 | Parameters |
||
538 | ---------- |
||
539 | nside : `~numpy.ndarray` |
||
540 | HEALPIX nside parameter, the total number of pixels is |
||
541 | 12*nside*nside. For multi-dimensional maps one can pass |
||
542 | either a single nside value or a vector of nside values |
||
543 | defining the pixel size for each image plane. If nside is not |
||
544 | a scalar then its dimensionality should match that of the |
||
545 | non-spatial axes. |
||
546 | nest : bool |
||
547 | True -> 'NESTED', False -> 'RING' indexing scheme |
||
548 | coordsys : str |
||
549 | Coordinate system, 'CEL' | 'GAL' |
||
550 | region : str or tuple |
||
551 | Spatial geometry for partial-sky maps. If none the map will |
||
552 | encompass the whole sky. String input will be parsed |
||
553 | according to HPX_REG header keyword conventions. Tuple |
||
554 | input can be used to define an explicit list of pixels |
||
555 | encompassed by the geometry. |
||
556 | axes : list |
||
557 | Axes for non-spatial dimensions. |
||
558 | conv : str |
||
559 | Convention for FITS serialization format. |
||
560 | sparse : bool |
||
561 | If True defer allocation of partial- to all-sky index mapping |
||
562 | arrays. This option is only compatible with partial-sky maps |
||
563 | with an analytic geometry (e.g. DISK). |
||
564 | """ |
||
565 | |||
566 | is_hpx = True |
||
567 | |||
568 | def __init__( |
||
569 | self, |
||
570 | nside, |
||
571 | nest=True, |
||
572 | coordsys="CEL", |
||
573 | region=None, |
||
574 | axes=None, |
||
575 | conv="gadf", |
||
576 | sparse=False, |
||
577 | ): |
||
578 | |||
579 | # FIXME: Figure out what to do when sparse=True |
||
580 | # FIXME: Require NSIDE to be power of two when nest=True |
||
581 | |||
582 | self._nside = np.array(nside, ndmin=1) |
||
583 | self._axes = make_axes(axes, conv) |
||
584 | self._shape = tuple([ax.nbin for ax in self._axes]) |
||
585 | if self.nside.size > 1 and self.nside.shape != self._shape: |
||
586 | raise ValueError( |
||
587 | "Wrong dimensionality for nside. nside must " |
||
588 | "be a scalar or have a dimensionality consistent " |
||
589 | "with the axes argument." |
||
590 | ) |
||
591 | |||
592 | self._order = nside_to_order(self._nside) |
||
593 | self._nest = nest |
||
594 | self._coordsys = coordsys |
||
595 | self._maxpix = 12 * self._nside * self._nside |
||
596 | self._maxpix = self._maxpix * np.ones(self._shape, dtype=int) |
||
597 | self._sparse = sparse |
||
598 | |||
599 | self._ipix = None |
||
600 | self._rmap = None |
||
601 | self._region = region |
||
602 | self._create_lookup(region) |
||
603 | |||
604 | if self._ipix is not None: |
||
605 | self._rmap = {} |
||
606 | for i, ipix in enumerate(self._ipix.flat): |
||
607 | self._rmap[ipix] = i |
||
608 | |||
609 | self._npix = self._npix * np.ones(self._shape, dtype=int) |
||
610 | self._conv = conv |
||
611 | self._center_skydir = self._get_ref_dir() |
||
612 | lon, lat, frame = skycoord_to_lonlat(self._center_skydir) |
||
613 | self._center_coord = tuple( |
||
614 | [lon, lat] |
||
615 | + [ax.pix_to_coord((float(ax.nbin) - 1.0) / 2.0) for ax in self.axes] |
||
616 | ) |
||
617 | self._center_pix = self.coord_to_pix(self._center_coord) |
||
618 | |||
619 | @property |
||
620 | def data_shape(self): |
||
621 | """Shape of the Numpy data array matching this geometry.""" |
||
622 | npix_shape = [np.max(self.npix)] |
||
623 | ax_shape = [ax.nbin for ax in self.axes] |
||
624 | return tuple(npix_shape + ax_shape)[::-1] |
||
625 | |||
626 | def _create_lookup(self, region): |
||
627 | """Create local-to-global pixel lookup table.""" |
||
628 | |||
629 | if isinstance(region, six.string_types): |
||
630 | ipix = [ |
||
631 | self.get_index_list(nside, self._nest, region) |
||
632 | for nside in self._nside.flat |
||
633 | ] |
||
634 | self._ibnd = np.concatenate( |
||
635 | [i * np.ones_like(p, dtype="int16") for i, p in enumerate(ipix)] |
||
636 | ) |
||
637 | self._ipix = [ |
||
638 | ravel_hpx_index((p, i * np.ones_like(p)), np.ravel(self._maxpix)) |
||
639 | for i, p in enumerate(ipix) |
||
640 | ] |
||
641 | self._region = region |
||
642 | self._indxschm = "EXPLICIT" |
||
643 | self._npix = np.array([len(t) for t in self._ipix]) |
||
644 | if self.nside.ndim > 1: |
||
645 | self._npix = self._npix.reshape(self.nside.shape) |
||
646 | self._ipix = np.concatenate(self._ipix) |
||
647 | |||
648 | elif isinstance(region, tuple): |
||
649 | |||
650 | region = [np.asarray(t) for t in region] |
||
651 | m = np.any(np.stack([t >= 0 for t in region]), axis=0) |
||
652 | region = [t[m] for t in region] |
||
653 | |||
654 | self._ipix = ravel_hpx_index(region, self._maxpix) |
||
655 | self._ipix = np.unique(self._ipix) |
||
656 | region = unravel_hpx_index(self._ipix, self._maxpix) |
||
657 | self._region = "explicit" |
||
658 | self._indxschm = "EXPLICIT" |
||
659 | if len(region) == 1: |
||
660 | self._npix = np.array([len(region[0])]) |
||
661 | else: |
||
662 | self._npix = np.zeros(self._shape, dtype=int) |
||
663 | idx = np.ravel_multi_index(region[1:], self._shape) |
||
664 | cnt = np.unique(idx, return_counts=True) |
||
665 | self._npix.flat[cnt[0]] = cnt[1] |
||
666 | |||
667 | elif region is None: |
||
668 | self._region = None |
||
669 | self._indxschm = "IMPLICIT" |
||
670 | self._npix = self._maxpix |
||
671 | |||
672 | else: |
||
673 | raise ValueError("Invalid region string: {!r}".format(region)) |
||
674 | |||
675 | def local_to_global(self, idx_local): |
||
676 | """Compute a local index (partial-sky) from a global (all-sky) |
||
677 | index. |
||
678 | |||
679 | Returns |
||
680 | ------- |
||
681 | idx_global : tuple |
||
682 | A tuple of pixel index vectors with global HEALPIX pixel |
||
683 | indices. |
||
684 | """ |
||
685 | if self._ipix is None: |
||
686 | return idx_local |
||
687 | |||
688 | if self.nside.size > 1: |
||
689 | idx = ravel_hpx_index(idx_local, self._npix) |
||
690 | else: |
||
691 | idx_tmp = tuple( |
||
692 | [idx_local[0]] + [np.zeros(t.shape, dtype=int) for t in idx_local[1:]] |
||
693 | ) |
||
694 | idx = ravel_hpx_index(idx_tmp, self._npix) |
||
695 | |||
696 | idx_global = unravel_hpx_index(self._ipix[idx], self._maxpix) |
||
697 | return idx_global[:1] + tuple(idx_local[1:]) |
||
698 | |||
699 | def global_to_local(self, idx_global, ravel=False): |
||
700 | """Compute global (all-sky) index from a local (partial-sky) index. |
||
701 | |||
702 | Parameters |
||
703 | ---------- |
||
704 | idx_global : tuple |
||
705 | A tuple of pixel indices with global HEALPix pixel indices. |
||
706 | ravel : bool |
||
707 | Return a raveled index. |
||
708 | |||
709 | Returns |
||
710 | ------- |
||
711 | idx_local : tuple |
||
712 | A tuple of pixel indices with local HEALPIX pixel indices. |
||
713 | """ |
||
714 | if self.nside.size == 1: |
||
715 | idx = np.array(idx_global[0], ndmin=1) |
||
716 | else: |
||
717 | idx = ravel_hpx_index(idx_global, self._maxpix) |
||
718 | |||
719 | if self._rmap is not None: |
||
720 | retval = np.full(idx.size, -1, "i") |
||
721 | m = np.in1d(idx.flat, self._ipix) |
||
722 | retval[m] = np.searchsorted(self._ipix, idx.flat[m]) |
||
723 | retval = retval.reshape(idx.shape) |
||
724 | else: |
||
725 | retval = idx |
||
726 | |||
727 | if self.nside.size == 1: |
||
728 | idx_local = tuple([retval] + list(idx_global[1:])) |
||
729 | else: |
||
730 | idx_local = unravel_hpx_index(retval, self._npix) |
||
731 | |||
732 | m = np.any(np.stack([t == -1 for t in idx_local]), axis=0) |
||
733 | for i, t in enumerate(idx_local): |
||
734 | idx_local[i][m] = -1 |
||
735 | |||
736 | if not ravel: |
||
737 | return idx_local |
||
738 | else: |
||
739 | return ravel_hpx_index(idx_local, self.npix) |
||
740 | |||
741 | def __getitem__(self, idx_global): |
||
742 | """This implements the global-to-local index lookup. |
||
743 | |||
744 | For all-sky maps it just returns the input array. For |
||
745 | partial-sky maps it returns the local indices corresponding to |
||
746 | the indices in the input array, and -1 for those pixels that |
||
747 | are outside the selected region. For multi-dimensional maps |
||
748 | with a different ``NSIDE`` in each band the global index is an |
||
749 | unrolled index for both HEALPIX pixel number and image slice. |
||
750 | |||
751 | Parameters |
||
752 | ---------- |
||
753 | idx_global : `~numpy.ndarray` |
||
754 | An array of global (all-sky) pixel indices. If this is a |
||
755 | tuple, list, or array of integers it will be interpreted |
||
756 | as a global (raveled) index. If this argument is a tuple |
||
757 | of lists or arrays it will be interpreted as a list of |
||
758 | unraveled index vectors. |
||
759 | |||
760 | Returns |
||
761 | ------- |
||
762 | idx_local : `~numpy.ndarray` |
||
763 | An array of local HEALPIX pixel indices. |
||
764 | """ |
||
765 | # Convert to tuple representation |
||
766 | if ( |
||
767 | isinstance(idx_global, int) |
||
768 | or (isinstance(idx_global, tuple) and isinstance(idx_global[0], int)) |
||
769 | or isinstance(idx_global, np.ndarray) |
||
770 | ): |
||
771 | idx_global = unravel_hpx_index(np.array(idx_global, ndmin=1), self._maxpix) |
||
772 | |||
773 | return self.global_to_local(idx_global, ravel=True) |
||
774 | |||
775 | def coord_to_pix(self, coords): |
||
776 | import healpy as hp |
||
777 | |||
778 | coords = MapCoord.create(coords, coordsys=self.coordsys) |
||
779 | theta, phi = coords.theta, coords.phi |
||
780 | |||
781 | c = self.coord_to_tuple(coords) |
||
782 | |||
783 | if self.axes: |
||
784 | bins = [] |
||
785 | idxs = [] |
||
786 | for i, ax in enumerate(self.axes): |
||
787 | bins += [ax.coord_to_pix(c[i + 2])] |
||
788 | idxs += [ax.coord_to_idx(c[i + 2])] |
||
789 | |||
790 | # FIXME: Figure out how to handle coordinates out of |
||
791 | # bounds of non-spatial dimensions |
||
792 | if self.nside.size > 1: |
||
793 | nside = self.nside[tuple(idxs)] |
||
794 | else: |
||
795 | nside = self.nside |
||
796 | |||
797 | m = ~np.isfinite(theta) |
||
798 | theta[m] = 0.0 |
||
799 | phi[m] = 0.0 |
||
800 | pix = hp.ang2pix(nside, theta, phi, nest=self.nest) |
||
801 | pix = tuple([pix] + bins) |
||
802 | if np.any(m): |
||
803 | for p in pix: |
||
804 | p[m] = -1 |
||
805 | else: |
||
806 | pix = (hp.ang2pix(self.nside, theta, phi, nest=self.nest),) |
||
807 | |||
808 | return pix |
||
809 | |||
810 | def pix_to_coord(self, pix): |
||
811 | import healpy as hp |
||
812 | |||
813 | if self.axes: |
||
814 | bins = [] |
||
815 | vals = [] |
||
816 | for i, ax in enumerate(self.axes): |
||
817 | bins += [pix[1 + i]] |
||
818 | vals += [ax.pix_to_coord(pix[1 + i])] |
||
819 | |||
820 | idxs = pix_tuple_to_idx(bins) |
||
821 | |||
822 | if self.nside.size > 1: |
||
823 | nside = self.nside[idxs] |
||
824 | else: |
||
825 | nside = self.nside |
||
826 | |||
827 | ipix = np.round(pix[0]).astype(int) |
||
828 | m = ipix == -1 |
||
829 | ipix[m] = 0 |
||
830 | theta, phi = hp.pix2ang(nside, ipix, nest=self.nest) |
||
831 | coords = [np.degrees(phi), np.degrees(np.pi / 2.0 - theta)] |
||
832 | coords = tuple(coords + vals) |
||
833 | if np.any(m): |
||
834 | for c in coords: |
||
835 | c[m] = np.nan |
||
836 | else: |
||
837 | ipix = np.round(pix[0]).astype(int) |
||
838 | theta, phi = hp.pix2ang(self.nside, ipix, nest=self.nest) |
||
839 | coords = (np.degrees(phi), np.degrees(np.pi / 2.0 - theta)) |
||
840 | |||
841 | return coords |
||
842 | |||
843 | def pix_to_idx(self, pix, clip=False): |
||
844 | # FIXME: Look for better method to clip HPX indices |
||
845 | # TODO: copy idx to avoid modifying input pix? |
||
846 | # pix_tuple_to_idx seems to always make a copy!? |
||
847 | idx = pix_tuple_to_idx(pix) |
||
848 | idx_local = self.global_to_local(idx) |
||
849 | for i, _ in enumerate(idx): |
||
850 | |||
851 | if clip: |
||
852 | if i > 0: |
||
853 | np.clip(idx[i], 0, self.axes[i - 1].nbin - 1, out=idx[i]) |
||
854 | else: |
||
855 | np.clip(idx[i], 0, None, out=idx[i]) |
||
856 | else: |
||
857 | if i > 0: |
||
858 | mask = (idx[i] < 0) | (idx[i] >= self.axes[i - 1].nbin) |
||
859 | np.putmask(idx[i], mask, -1) |
||
860 | else: |
||
861 | mask = (idx_local[i] < 0) | (idx[i] < 0) |
||
862 | np.putmask(idx[i], mask, -1) |
||
863 | |||
864 | return tuple(idx) |
||
865 | |||
866 | def to_slice(self, slices, drop_axes=True): |
||
867 | if len(slices) == 0 and self.ndim == 2: |
||
868 | return copy.deepcopy(self) |
||
869 | |||
870 | if len(slices) != self.ndim - 2: |
||
871 | raise ValueError() |
||
872 | |||
873 | nside = np.ones(self.shape, dtype=int) * self.nside |
||
874 | nside = np.squeeze(nside[slices]) |
||
875 | |||
876 | axes = [ax.slice(s) for ax, s in zip(self.axes, slices)] |
||
877 | if drop_axes: |
||
878 | axes = [ax for ax in axes if ax.nbin > 1] |
||
879 | slice_dims = [0] + [i + 1 for i, ax in enumerate(axes) if ax.nbin > 1] |
||
880 | else: |
||
881 | slice_dims = np.arange(self.ndim) |
||
882 | |||
883 | if self.region == "explicit": |
||
884 | idx = self.get_idx() |
||
885 | slices = (slice(None),) + slices |
||
886 | idx = [p[slices[::-1]] for p in idx] |
||
887 | idx = [p[p != -1] for p in idx] |
||
888 | if drop_axes: |
||
889 | idx = [idx[i] for i in range(len(idx)) if i in slice_dims] |
||
890 | region = tuple(idx) |
||
891 | else: |
||
892 | region = self.region |
||
893 | |||
894 | return self.__class__(nside, self.nest, self.coordsys, region, axes, self.conv) |
||
895 | |||
896 | @property |
||
897 | def axes(self): |
||
898 | """List of non-spatial axes.""" |
||
899 | return self._axes |
||
900 | |||
901 | @property |
||
902 | def shape(self): |
||
903 | """Shape of non-spatial axes.""" |
||
904 | return self._shape |
||
905 | |||
906 | @property |
||
907 | def ndim(self): |
||
908 | """Number of dimensions (int).""" |
||
909 | return len(self._axes) + 2 |
||
910 | |||
911 | @property |
||
912 | def ordering(self): |
||
913 | """HEALPix ordering ('NESTED' or 'RING').""" |
||
914 | return "NESTED" if self.nest else "RING" |
||
915 | |||
916 | @property |
||
917 | def nside(self): |
||
918 | """NSIDE in each band.""" |
||
919 | return self._nside |
||
920 | |||
921 | @property |
||
922 | def order(self): |
||
923 | """ORDER in each band (NSIDE = 2**ORDER). Set to -1 for bands with |
||
924 | NSIDE that is not a power of 2. |
||
925 | """ |
||
926 | return self._order |
||
927 | |||
928 | @property |
||
929 | def nest(self): |
||
930 | """Is HEALPix order nested? (bool).""" |
||
931 | return self._nest |
||
932 | |||
933 | @property |
||
934 | def npix(self): |
||
935 | """Number of pixels in each band. |
||
936 | |||
937 | For partial-sky geometries this can |
||
938 | be less than the number of pixels for the band NSIDE. |
||
939 | """ |
||
940 | return self._npix |
||
941 | |||
942 | @property |
||
943 | def conv(self): |
||
944 | """Name of default FITS convention associated with this geometry.""" |
||
945 | return self._conv |
||
946 | |||
947 | @property |
||
948 | def hpx_conv(self): |
||
949 | return HPX_FITS_CONVENTIONS[self.conv] |
||
950 | |||
951 | @property |
||
952 | def coordsys(self): |
||
953 | return self._coordsys |
||
954 | |||
955 | @property |
||
956 | def projection(self): |
||
957 | """Map projection.""" |
||
958 | return "HPX" |
||
959 | |||
960 | @property |
||
961 | def region(self): |
||
962 | """Region string.""" |
||
963 | return self._region |
||
964 | |||
965 | @property |
||
966 | def is_allsky(self): |
||
967 | """Flag for all-sky maps.""" |
||
968 | if self._region is None: |
||
969 | return True |
||
970 | else: |
||
971 | return False |
||
972 | |||
973 | @property |
||
974 | def is_regular(self): |
||
975 | """Flag identifying whether this geometry is regular in non-spatial dimensions. |
||
976 | |||
977 | False for multi-resolution or irregular geometries. |
||
978 | If true all image planes have the same pixel geometry. |
||
979 | """ |
||
980 | if self.nside.size > 1 or self.region == "explicit": |
||
981 | return False |
||
982 | else: |
||
983 | return True |
||
984 | |||
985 | @property |
||
986 | def center_coord(self): |
||
987 | """Map coordinates of the center of the geometry (tuple).""" |
||
988 | return self._center_coord |
||
989 | |||
990 | @property |
||
991 | def center_pix(self): |
||
992 | """Pixel coordinates of the center of the geometry (tuple).""" |
||
993 | return self._center_pix |
||
994 | |||
995 | @property |
||
996 | def center_skydir(self): |
||
997 | """Sky coordinate of the center of the geometry. |
||
998 | |||
999 | Returns |
||
1000 | ------- |
||
1001 | pix : `~astropy.coordinates.SkyCoord` |
||
1002 | """ |
||
1003 | return self._center_skydir |
||
1004 | |||
1005 | @property |
||
1006 | def ipix(self): |
||
1007 | """HEALPIX pixel and band indices for every pixel in the map.""" |
||
1008 | return self.get_idx() |
||
1009 | |||
1010 | def to_ud_graded(self, order): |
||
1011 | """Upgrade or downgrade the resolution of this geometry to the given |
||
1012 | order. This method does not preserve the geometry footprint. |
||
1013 | |||
1014 | Returns |
||
1015 | ------- |
||
1016 | geom : `~HpxGeom` |
||
1017 | A HEALPix geometry object. |
||
1018 | """ |
||
1019 | if np.any(self.order < 0): |
||
1020 | raise ValueError("Upgrade and degrade only implemented for standard maps") |
||
1021 | |||
1022 | axes = copy.deepcopy(self.axes) |
||
1023 | return self.__class__( |
||
1024 | 2 ** order, |
||
1025 | self.nest, |
||
1026 | coordsys=self.coordsys, |
||
1027 | region=self.region, |
||
1028 | axes=axes, |
||
1029 | conv=self.conv, |
||
1030 | ) |
||
1031 | |||
1032 | def to_swapped(self): |
||
1033 | """Make a copy of this geometry with a swapped ORDERING. (NEST->RING |
||
1034 | or vice versa) |
||
1035 | |||
1036 | Returns |
||
1037 | ------- |
||
1038 | geom : `~HpxGeom` |
||
1039 | A HEALPix geometry object. |
||
1040 | """ |
||
1041 | axes = copy.deepcopy(self.axes) |
||
1042 | return self.__class__( |
||
1043 | self.nside, |
||
1044 | not self.nest, |
||
1045 | coordsys=self.coordsys, |
||
1046 | region=self.region, |
||
1047 | axes=axes, |
||
1048 | conv=self.conv, |
||
1049 | ) |
||
1050 | |||
1051 | def to_image(self): |
||
1052 | return self.__class__( |
||
1053 | np.max(self.nside), |
||
1054 | self.nest, |
||
1055 | coordsys=self.coordsys, |
||
1056 | region=self.region, |
||
1057 | conv=self.conv, |
||
1058 | ) |
||
1059 | |||
1060 | def to_cube(self, axes): |
||
1061 | axes = copy.deepcopy(self.axes) + axes |
||
1062 | return self.__class__( |
||
1063 | np.max(self.nside), |
||
1064 | self.nest, |
||
1065 | coordsys=self.coordsys, |
||
1066 | region=self.region, |
||
1067 | conv=self.conv, |
||
1068 | axes=axes, |
||
1069 | ) |
||
1070 | |||
1071 | def _get_neighbors(self, idx): |
||
1072 | import healpy as hp |
||
1073 | |||
1074 | nside = self._get_nside(idx) |
||
1075 | idx_nb = (hp.get_all_neighbours(nside, idx[0], nest=self.nest),) |
||
1076 | idx_nb += tuple([t[None, ...] * np.ones_like(idx_nb[0]) for t in idx[1:]]) |
||
1077 | |||
1078 | return idx_nb |
||
1079 | |||
1080 | def pad(self, pad_width): |
||
1081 | if self.is_allsky: |
||
1082 | raise ValueError("Cannot pad an all-sky map.") |
||
1083 | |||
1084 | idx = self.get_idx(flat=True) |
||
1085 | idx_r = ravel_hpx_index(idx, self._maxpix) |
||
1086 | |||
1087 | # TODO: Pre-filter indices to find those close to the edge |
||
1088 | idx_nb = self._get_neighbors(idx) |
||
1089 | idx_nb = ravel_hpx_index(idx_nb, self._maxpix) |
||
1090 | |||
1091 | for _ in range(pad_width): |
||
1092 | # Mask of neighbors that are not contained in the geometry |
||
1093 | # TODO: change this to numpy.isin when we require Numpy 1.13+ |
||
1094 | # Here and everywhere in Gamampy -> search for "isin" |
||
1095 | # see https://github.com/gammapy/gammapy/pull/1371 |
||
1096 | mask_edge = np.in1d(idx_nb, idx_r, invert=True).reshape(idx_nb.shape) |
||
1097 | idx_edge = idx_nb[mask_edge] |
||
1098 | idx_edge = np.unique(idx_edge) |
||
1099 | idx_r = np.sort(np.concatenate((idx_r, idx_edge))) |
||
1100 | idx_nb = unravel_hpx_index(idx_edge, self._maxpix) |
||
1101 | idx_nb = self._get_neighbors(idx_nb) |
||
1102 | idx_nb = ravel_hpx_index(idx_nb, self._maxpix) |
||
1103 | |||
1104 | idx = unravel_hpx_index(idx_r, self._maxpix) |
||
1105 | return self.__class__( |
||
1106 | self.nside.copy(), |
||
1107 | self.nest, |
||
1108 | coordsys=self.coordsys, |
||
1109 | region=idx, |
||
1110 | conv=self.conv, |
||
1111 | axes=copy.deepcopy(self.axes), |
||
1112 | ) |
||
1113 | |||
1114 | def crop(self, crop_width): |
||
1115 | if self.is_allsky: |
||
1116 | raise ValueError("Cannot crop an all-sky map.") |
||
1117 | |||
1118 | idx = self.get_idx(flat=True) |
||
1119 | idx_r = ravel_hpx_index(idx, self._maxpix) |
||
1120 | |||
1121 | # TODO: Pre-filter indices to find those close to the edge |
||
1122 | idx_nb = self._get_neighbors(idx) |
||
1123 | idx_nb = ravel_hpx_index(idx_nb, self._maxpix) |
||
1124 | |||
1125 | for _ in range(crop_width): |
||
1126 | # Mask of pixels that have at least one neighbor not |
||
1127 | # contained in the geometry |
||
1128 | mask_edge = np.in1d(idx_nb, idx_r, invert=True) |
||
1129 | mask_edge = np.any(mask_edge, axis=0) |
||
1130 | idx_r = idx_r[~mask_edge] |
||
1131 | idx_nb = idx_nb[:, ~mask_edge] |
||
1132 | |||
1133 | idx = unravel_hpx_index(idx_r, self._maxpix) |
||
1134 | return self.__class__( |
||
1135 | self.nside.copy(), |
||
1136 | self.nest, |
||
1137 | coordsys=self.coordsys, |
||
1138 | region=idx, |
||
1139 | conv=self.conv, |
||
1140 | axes=copy.deepcopy(self.axes), |
||
1141 | ) |
||
1142 | |||
1143 | def upsample(self, factor): |
||
1144 | if not is_power2(factor): |
||
1145 | raise ValueError("Upsample factor must be a power of 2.") |
||
1146 | |||
1147 | if self.is_allsky: |
||
1148 | return self.__class__( |
||
1149 | self.nside * factor, |
||
1150 | self.nest, |
||
1151 | coordsys=self.coordsys, |
||
1152 | region=self.region, |
||
1153 | conv=self.conv, |
||
1154 | axes=copy.deepcopy(self.axes), |
||
1155 | ) |
||
1156 | |||
1157 | idx = list(self.get_idx(flat=True)) |
||
1158 | nside = self._get_nside(idx) |
||
1159 | |||
1160 | idx_new = get_subpixels(idx[0], nside, nside * factor, nest=self.nest) |
||
1161 | for i in range(1, len(idx)): |
||
1162 | idx[i] = idx[i][..., None] * np.ones(idx_new.shape, dtype=int) |
||
1163 | |||
1164 | idx[0] = idx_new |
||
1165 | return self.__class__( |
||
1166 | self.nside * factor, |
||
1167 | self.nest, |
||
1168 | coordsys=self.coordsys, |
||
1169 | region=tuple(idx), |
||
1170 | conv=self.conv, |
||
1171 | axes=copy.deepcopy(self.axes), |
||
1172 | ) |
||
1173 | |||
1174 | def downsample(self, factor): |
||
1175 | if not is_power2(factor): |
||
1176 | raise ValueError("Downsample factor must be a power of 2.") |
||
1177 | |||
1178 | if self.is_allsky: |
||
1179 | return self.__class__( |
||
1180 | self.nside // factor, |
||
1181 | self.nest, |
||
1182 | coordsys=self.coordsys, |
||
1183 | region=self.region, |
||
1184 | conv=self.conv, |
||
1185 | axes=copy.deepcopy(self.axes), |
||
1186 | ) |
||
1187 | |||
1188 | idx = list(self.get_idx(flat=True)) |
||
1189 | nside = self._get_nside(idx) |
||
1190 | idx_new = get_superpixels(idx[0], nside, nside // factor, nest=self.nest) |
||
1191 | idx[0] = idx_new |
||
1192 | return self.__class__( |
||
1193 | self.nside // factor, |
||
1194 | self.nest, |
||
1195 | coordsys=self.coordsys, |
||
1196 | region=tuple(idx), |
||
1197 | conv=self.conv, |
||
1198 | axes=copy.deepcopy(self.axes), |
||
1199 | ) |
||
1200 | |||
1201 | @classmethod |
||
1202 | def create( |
||
1203 | cls, |
||
1204 | nside=None, |
||
1205 | binsz=None, |
||
1206 | nest=True, |
||
1207 | coordsys="CEL", |
||
1208 | region=None, |
||
1209 | axes=None, |
||
1210 | conv="gadf", |
||
1211 | skydir=None, |
||
1212 | width=None, |
||
1213 | ): |
||
1214 | """Create an HpxGeom object. |
||
1215 | |||
1216 | Parameters |
||
1217 | ---------- |
||
1218 | nside : int or `~numpy.ndarray` |
||
1219 | HEALPix NSIDE parameter. This parameter sets the size of |
||
1220 | the spatial pixels in the map. |
||
1221 | binsz : float or `~numpy.ndarray` |
||
1222 | Approximate pixel size in degrees. An NSIDE will be |
||
1223 | chosen that correponds to a pixel size closest to this |
||
1224 | value. This option is superseded by nside. |
||
1225 | nest : bool |
||
1226 | True for HEALPIX "NESTED" indexing scheme, False for "RING" scheme |
||
1227 | coordsys : {'CEL', 'GAL'}, optional |
||
1228 | Coordinate system, either Galactic ('GAL') or Equatorial ('CEL'). |
||
1229 | skydir : tuple or `~astropy.coordinates.SkyCoord` |
||
1230 | Sky position of map center. Can be either a SkyCoord |
||
1231 | object or a tuple of longitude and latitude in deg in the |
||
1232 | coordinate system of the map. |
||
1233 | region : str |
||
1234 | HPX region string. Allows for partial-sky maps. |
||
1235 | width : float |
||
1236 | Diameter of the map in degrees. If set the map will |
||
1237 | encompass all pixels within a circular region centered on |
||
1238 | ``skydir``. |
||
1239 | axes : list |
||
1240 | List of axes for non-spatial dimensions. |
||
1241 | conv : str |
||
1242 | Convention for FITS file format. |
||
1243 | |||
1244 | Returns |
||
1245 | ------- |
||
1246 | geom : `~HpxGeom` |
||
1247 | A HEALPix geometry object. |
||
1248 | |||
1249 | Examples |
||
1250 | -------- |
||
1251 | >>> from gammapy.maps import HpxGeom, MapAxis |
||
1252 | >>> axis = MapAxis.from_bounds(0,1,2) |
||
1253 | >>> geom = HpxGeom.create(nside=16) |
||
1254 | >>> geom = HpxGeom.create(binsz=0.1, width=10.0) |
||
1255 | >>> geom = HpxGeom.create(nside=64, width=10.0, axes=[axis]) |
||
1256 | >>> geom = HpxGeom.create(nside=[32,64], width=10.0, axes=[axis]) |
||
1257 | """ |
||
1258 | if nside is None and binsz is None: |
||
1259 | raise ValueError("Either nside or binsz must be defined.") |
||
1260 | |||
1261 | if nside is None and binsz is not None: |
||
1262 | nside = get_nside_from_pix_size(binsz) |
||
1263 | |||
1264 | if skydir is None: |
||
1265 | lon, lat = (0.0, 0.0) |
||
1266 | elif isinstance(skydir, tuple): |
||
1267 | lon, lat = skydir |
||
1268 | elif isinstance(skydir, SkyCoord): |
||
1269 | lon, lat, frame = skycoord_to_lonlat(skydir, coordsys=coordsys) |
||
1270 | else: |
||
1271 | raise ValueError("Invalid type for skydir: {!r}".format(type(skydir))) |
||
1272 | |||
1273 | if region is None and width is not None: |
||
1274 | region = "DISK({:f},{:f},{:f})".format(lon, lat, width / 2.0) |
||
1275 | |||
1276 | return cls( |
||
1277 | nside, nest=nest, coordsys=coordsys, region=region, axes=axes, conv=conv |
||
1278 | ) |
||
1279 | |||
1280 | @staticmethod |
||
1281 | def identify_hpx_convention(header): |
||
1282 | """Identify the convention used to write this file.""" |
||
1283 | # Hopefully the file contains the HPX_CONV keyword specifying |
||
1284 | # the convention used |
||
1285 | if "HPX_CONV" in header: |
||
1286 | return header["HPX_CONV"].lower() |
||
1287 | |||
1288 | # Try based on the EXTNAME keyword |
||
1289 | hduname = header.get("EXTNAME", None) |
||
1290 | if hduname == "HPXEXPOSURES": |
||
1291 | return "fgst-bexpcube" |
||
1292 | elif hduname == "SKYMAP2": |
||
1293 | if "COORDTYPE" in header.keys(): |
||
1294 | return "galprop" |
||
1295 | else: |
||
1296 | return "galprop2" |
||
1297 | |||
1298 | # Check the name of the first column |
||
1299 | colname = header["TTYPE1"] |
||
1300 | if colname == "PIX": |
||
1301 | colname = header["TTYPE2"] |
||
1302 | |||
1303 | if colname == "KEY": |
||
1304 | return "fgst-srcmap-sparse" |
||
1305 | elif colname == "ENERGY1": |
||
1306 | return "fgst-template" |
||
1307 | elif colname == "COSBINS": |
||
1308 | return "fgst-ltcube" |
||
1309 | elif colname == "Bin0": |
||
1310 | return "galprop" |
||
1311 | elif colname == "CHANNEL1" or colname == "CHANNEL0": |
||
1312 | if hduname == "SKYMAP": |
||
1313 | return "fgst-ccube" |
||
1314 | else: |
||
1315 | return "fgst-srcmap" |
||
1316 | else: |
||
1317 | raise ValueError("Could not identify HEALPIX convention") |
||
1318 | |||
1319 | @classmethod |
||
1320 | def from_header(cls, header, hdu_bands=None, pix=None): |
||
1321 | """Create an HPX object from a FITS header. |
||
1322 | |||
1323 | Parameters |
||
1324 | ---------- |
||
1325 | header : `~astropy.io.fits.Header` |
||
1326 | The FITS header |
||
1327 | hdu_bands : `~astropy.io.fits.BinTableHDU` |
||
1328 | The BANDS table HDU. |
||
1329 | pix : tuple |
||
1330 | List of pixel index vectors defining the pixels |
||
1331 | encompassed by the geometry. For EXPLICIT geometries with |
||
1332 | HPX_REG undefined this tuple defines the geometry. |
||
1333 | |||
1334 | Returns |
||
1335 | ------- |
||
1336 | hpx : `~HpxGeom` |
||
1337 | HEALPix geometry. |
||
1338 | """ |
||
1339 | convname = HpxGeom.identify_hpx_convention(header) |
||
1340 | conv = HPX_FITS_CONVENTIONS[convname] |
||
1341 | |||
1342 | axes = find_and_read_bands(hdu_bands) |
||
1343 | shape = [ax.nbin for ax in axes] |
||
1344 | |||
1345 | if header["PIXTYPE"] != "HEALPIX": |
||
1346 | raise ValueError( |
||
1347 | "Header PIXTYPE must be 'HEALPIX'. Got: {}".format(header["PIXTYPE"]) |
||
1348 | ) |
||
1349 | |||
1350 | if header["ORDERING"] == "RING": |
||
1351 | nest = False |
||
1352 | elif header["ORDERING"] == "NESTED": |
||
1353 | nest = True |
||
1354 | else: |
||
1355 | raise ValueError( |
||
1356 | "Header ORDERING must be RING or NESTED. Got: {}".format( |
||
1357 | header["ORDERING"] |
||
1358 | ) |
||
1359 | ) |
||
1360 | |||
1361 | if hdu_bands is not None and "NSIDE" in hdu_bands.columns.names: |
||
1362 | nside = hdu_bands.data.field("NSIDE").reshape(shape).astype(int) |
||
1363 | elif "NSIDE" in header: |
||
1364 | nside = header["NSIDE"] |
||
1365 | elif "ORDER" in header: |
||
1366 | nside = 2 ** header["ORDER"] |
||
1367 | else: |
||
1368 | raise ValueError("Failed to extract NSIDE or ORDER.") |
||
1369 | |||
1370 | try: |
||
1371 | coordsys = header[conv.coordsys] |
||
1372 | except KeyError: |
||
1373 | coordsys = header.get("COORDSYS", "CEL") |
||
1374 | |||
1375 | try: |
||
1376 | region = header["HPX_REG"] |
||
1377 | except KeyError: |
||
1378 | try: |
||
1379 | region = header["HPXREGION"] |
||
1380 | except KeyError: |
||
1381 | region = None |
||
1382 | |||
1383 | return cls( |
||
1384 | nside, nest, coordsys=coordsys, region=region, axes=axes, conv=convname |
||
1385 | ) |
||
1386 | |||
1387 | @classmethod |
||
1388 | def from_hdu(cls, hdu, hdu_bands=None): |
||
1389 | """Create an HPX object from a BinTable HDU. |
||
1390 | |||
1391 | Parameters |
||
1392 | ---------- |
||
1393 | hdu : `~astropy.io.fits.BinTableHDU` |
||
1394 | The FITS HDU |
||
1395 | hdu_bands : `~astropy.io.fits.BinTableHDU` |
||
1396 | The BANDS table HDU |
||
1397 | |||
1398 | Returns |
||
1399 | ------- |
||
1400 | hpx : `~HpxGeom` |
||
1401 | HEALPix geometry. |
||
1402 | """ |
||
1403 | # FIXME: Need correct handling of IMPLICIT and EXPLICIT maps |
||
1404 | |||
1405 | # if HPX region is not defined then geometry is defined by |
||
1406 | # the set of all pixels in the table |
||
1407 | if "HPX_REG" not in hdu.header: |
||
1408 | pix = (hdu.data.field("PIX"), hdu.data.field("CHANNEL")) |
||
1409 | else: |
||
1410 | pix = None |
||
1411 | |||
1412 | return cls.from_header(hdu.header, hdu_bands=hdu_bands, pix=pix) |
||
1413 | |||
1414 | def make_header(self, **kwargs): |
||
1415 | """"Build and return FITS header for this HEALPIX map.""" |
||
1416 | header = fits.Header() |
||
1417 | conv = kwargs.get("conv", HPX_FITS_CONVENTIONS[self.conv]) |
||
1418 | |||
1419 | # FIXME: For some sparse maps we may want to allow EXPLICIT |
||
1420 | # with an empty region string |
||
1421 | indxschm = kwargs.get("indxschm", None) |
||
1422 | |||
1423 | if indxschm is None: |
||
1424 | if self._region is None: |
||
1425 | indxschm = "IMPLICIT" |
||
1426 | elif self.is_regular == 1: |
||
1427 | indxschm = "EXPLICIT" |
||
1428 | else: |
||
1429 | indxschm = "LOCAL" |
||
1430 | |||
1431 | if "FGST" in conv.convname.upper(): |
||
1432 | header["TELESCOP"] = "GLAST" |
||
1433 | header["INSTRUME"] = "LAT" |
||
1434 | |||
1435 | header[conv.coordsys] = self.coordsys |
||
1436 | header["PIXTYPE"] = "HEALPIX" |
||
1437 | header["ORDERING"] = self.ordering |
||
1438 | header["INDXSCHM"] = indxschm |
||
1439 | header["ORDER"] = np.max(self._order) |
||
1440 | header["NSIDE"] = np.max(self._nside) |
||
1441 | header["FIRSTPIX"] = 0 |
||
1442 | header["LASTPIX"] = np.max(self._maxpix) - 1 |
||
1443 | header["HPX_CONV"] = conv.convname.upper() |
||
1444 | |||
1445 | if self.coordsys == "CEL": |
||
1446 | header["EQUINOX"] = (2000.0, "Equinox of RA & DEC specifications") |
||
1447 | |||
1448 | if self.region: |
||
1449 | header["HPX_REG"] = self._region |
||
1450 | |||
1451 | return header |
||
1452 | |||
1453 | def _make_bands_cols(self): |
||
1454 | cols = [] |
||
1455 | if self.nside.size > 1: |
||
1456 | cols += [fits.Column("NSIDE", "I", array=np.ravel(self.nside))] |
||
1457 | return cols |
||
1458 | |||
1459 | @staticmethod |
||
1460 | def get_index_list(nside, nest, region): |
||
1461 | """Get list of pixels indices for all the pixels in a region. |
||
1462 | |||
1463 | Parameters |
||
1464 | ---------- |
||
1465 | nside : int |
||
1466 | HEALPIX nside parameter |
||
1467 | nest : bool |
||
1468 | True for 'NESTED', False = 'RING' |
||
1469 | region : str |
||
1470 | HEALPIX region string |
||
1471 | |||
1472 | Returns |
||
1473 | ------- |
||
1474 | ilist : `~numpy.ndarray` |
||
1475 | List of pixel indices. |
||
1476 | """ |
||
1477 | import healpy as hp |
||
1478 | |||
1479 | # TODO: this should return something more friendly than a tuple |
||
1480 | # e.g. a namedtuple or a dict |
||
1481 | tokens = parse_hpxregion(region) |
||
1482 | |||
1483 | reg_type = tokens[0] |
||
1484 | if reg_type == "DISK": |
||
1485 | lon, lat = float(tokens[1]), float(tokens[2]) |
||
1486 | radius = np.radians(float(tokens[3])) |
||
1487 | vec = coords_to_vec(lon, lat)[0] |
||
1488 | ilist = hp.query_disc(nside, vec, radius, inclusive=False, nest=nest) |
||
1489 | elif reg_type == "DISK_INC": |
||
1490 | lon, lat = float(tokens[1]), float(tokens[2]) |
||
1491 | radius = np.radians(float(tokens[3])) |
||
1492 | vec = coords_to_vec(lon, lat)[0] |
||
1493 | fact = int(tokens[4]) |
||
1494 | ilist = hp.query_disc( |
||
1495 | nside, vec, radius, inclusive=True, nest=nest, fact=fact |
||
1496 | ) |
||
1497 | elif reg_type == "HPX_PIXEL": |
||
1498 | nside_pix = int(tokens[2]) |
||
1499 | if tokens[1] == "NESTED": |
||
1500 | ipix_ring = hp.nest2ring(nside_pix, int(tokens[3])) |
||
1501 | elif tokens[1] == "RING": |
||
1502 | ipix_ring = int(tokens[3]) |
||
1503 | else: |
||
1504 | raise ValueError("Invalid ordering scheme: {!r}".format(tokens[1])) |
||
1505 | ilist = match_hpx_pix(nside, nest, nside_pix, ipix_ring) |
||
1506 | else: |
||
1507 | raise ValueError("Invalid region type: {!r}".format(reg_type)) |
||
1508 | |||
1509 | return ilist |
||
1510 | |||
1511 | def _get_ref_dir(self): |
||
1512 | """Compute the reference direction for this geometry.""" |
||
1513 | import healpy as hp |
||
1514 | |||
1515 | frame = coordsys_to_frame(self.coordsys) |
||
1516 | |||
1517 | if self.region == "explicit": |
||
1518 | idx = unravel_hpx_index(self._ipix, self._maxpix) |
||
1519 | nside = self._get_nside(idx) |
||
1520 | vec = hp.pix2vec(nside, idx[0], nest=self.nest) |
||
1521 | vec = np.array([np.mean(t) for t in vec]) |
||
1522 | lonlat = hp.vec2ang(vec, lonlat=True) |
||
1523 | return SkyCoord(lonlat[0], lonlat[1], frame=frame, unit="deg") |
||
1524 | |||
1525 | return get_hpxregion_dir(self.region, self.coordsys) |
||
1526 | |||
1527 | def _get_region_size(self): |
||
1528 | import healpy as hp |
||
1529 | |||
1530 | if self.region is None: |
||
1531 | return 180.0 |
||
1532 | if self.region == "explicit": |
||
1533 | idx = unravel_hpx_index(self._ipix, self._maxpix) |
||
1534 | nside = self._get_nside(idx) |
||
1535 | ang = hp.pix2ang(nside, idx[0], nest=self.nest, lonlat=True) |
||
1536 | frame = coordsys_to_frame(self.coordsys) |
||
1537 | dirs = SkyCoord(ang[0], ang[1], unit="deg", frame=frame) |
||
1538 | return np.max(dirs.separation(self.center_skydir).deg) |
||
1539 | |||
1540 | return get_hpxregion_size(self.region) |
||
1541 | |||
1542 | def _get_nside(self, idx): |
||
1543 | if self.nside.size > 1: |
||
1544 | return self.nside[tuple(idx[1:])] |
||
1545 | else: |
||
1546 | return self.nside |
||
1547 | |||
1548 | def make_wcs(self, proj="AIT", oversample=2, drop_axes=True, width_pix=None): |
||
1549 | """Make a WCS projection appropriate for this HPX pixelization. |
||
1550 | |||
1551 | Parameters |
||
1552 | ---------- |
||
1553 | drop_axes : bool |
||
1554 | Drop non-spatial axes from the |
||
1555 | HEALPIX geometry. If False then all dimensions of the |
||
1556 | HEALPIX geometry will be copied to the WCS geometry. |
||
1557 | proj : str |
||
1558 | Projection type of WCS geometry. |
||
1559 | oversample : float |
||
1560 | Oversampling factor for WCS map. This will be the |
||
1561 | approximate ratio of the width of a HPX pixel to a WCS |
||
1562 | pixel. If this parameter is None then the width will be |
||
1563 | set from ``width_pix``. |
||
1564 | width_pix : int |
||
1565 | Width of the WCS geometry in pixels. The pixel size will |
||
1566 | be set to the number of pixels satisfying ``oversample`` |
||
1567 | or ``width_pix`` whichever is smaller. If this parameter |
||
1568 | is None then the width will be set from ``oversample``. |
||
1569 | |||
1570 | Returns |
||
1571 | ------- |
||
1572 | wcs : `~gammapy.maps.WcsGeom` |
||
1573 | WCS geometry |
||
1574 | """ |
||
1575 | pix_size = get_pix_size_from_nside(self.nside) |
||
1576 | binsz = np.min(pix_size) / oversample |
||
1577 | width = 2.0 * self._get_region_size() + np.max(pix_size) |
||
1578 | |||
1579 | if width_pix is not None and int(width / binsz) > width_pix: |
||
1580 | binsz = width / width_pix |
||
1581 | |||
1582 | if width > 90.0: |
||
1583 | width = min(360.0, width), min(180.0, width) |
||
1584 | |||
1585 | if drop_axes: |
||
1586 | axes = None |
||
1587 | else: |
||
1588 | axes = copy.deepcopy(self.axes) |
||
1589 | |||
1590 | return WcsGeom.create( |
||
1591 | width=width, |
||
1592 | binsz=binsz, |
||
1593 | coordsys=self.coordsys, |
||
1594 | axes=axes, |
||
1595 | skydir=self.center_skydir, |
||
1596 | proj=proj, |
||
1597 | ) |
||
1598 | |||
1599 | def get_idx(self, idx=None, local=False, flat=False): |
||
1600 | if idx is not None and np.any(np.array(idx) >= np.array(self._shape)): |
||
1601 | raise ValueError("Image index out of range: {!r}".format(idx)) |
||
1602 | |||
1603 | # Regular all- and partial-sky maps |
||
1604 | if self.is_regular: |
||
1605 | |||
1606 | pix = [np.arange(np.max(self._npix))] |
||
1607 | if idx is None: |
||
1608 | pix += [np.arange(ax.nbin, dtype=int) for ax in self.axes] |
||
1609 | else: |
||
1610 | pix += [t for t in idx] |
||
1611 | pix = np.meshgrid(*pix[::-1], indexing="ij", sparse=False)[::-1] |
||
1612 | pix = self.local_to_global(pix) |
||
1613 | |||
1614 | # Non-regular all-sky |
||
1615 | elif self.is_allsky and not self.is_regular: |
||
1616 | |||
1617 | shape = (np.max(self.npix),) |
||
1618 | if idx is None: |
||
1619 | shape = shape + self.shape |
||
1620 | else: |
||
1621 | shape = shape + (1,) * len(self.axes) |
||
1622 | pix = [np.full(shape, -1, dtype=int) for i in range(1 + len(self.axes))] |
||
1623 | for idx_img in np.ndindex(self.shape): |
||
1624 | |||
1625 | if idx is not None and idx_img != idx: |
||
1626 | continue |
||
1627 | |||
1628 | npix = self._npix[idx_img] |
||
1629 | if idx is None: |
||
1630 | s_img = (slice(0, npix),) + idx_img |
||
1631 | else: |
||
1632 | s_img = (slice(0, npix),) + (0,) * len(self.axes) |
||
1633 | |||
1634 | pix[0][s_img] = np.arange(self._npix[idx_img]) |
||
1635 | for j in range(len(self.axes)): |
||
1636 | pix[j + 1][s_img] = idx_img[j] |
||
1637 | pix = [p.T for p in pix] |
||
1638 | |||
1639 | # Explicit pixel indices |
||
1640 | else: |
||
1641 | |||
1642 | if idx is not None: |
||
1643 | npix_sum = np.concatenate(([0], np.cumsum(self._npix))) |
||
1644 | idx_ravel = np.ravel_multi_index(idx, self._shape) |
||
1645 | s = slice(npix_sum[idx_ravel], npix_sum[idx_ravel + 1]) |
||
1646 | else: |
||
1647 | s = slice(None) |
||
1648 | pix_flat = unravel_hpx_index(self._ipix[s], self._maxpix) |
||
1649 | |||
1650 | shape = (np.max(self.npix),) |
||
1651 | if idx is None: |
||
1652 | shape = shape + self.shape |
||
1653 | else: |
||
1654 | shape = shape + (1,) * len(self.axes) |
||
1655 | pix = [np.full(shape, -1, dtype=int) for _ in range(1 + len(self.axes))] |
||
1656 | |||
1657 | for idx_img in np.ndindex(self.shape): |
||
1658 | |||
1659 | if idx is not None and idx_img != idx: |
||
1660 | continue |
||
1661 | |||
1662 | npix = int(self._npix[idx_img]) |
||
1663 | if idx is None: |
||
1664 | s_img = (slice(0, npix),) + idx_img |
||
1665 | else: |
||
1666 | s_img = (slice(0, npix),) + (0,) * len(self.axes) |
||
1667 | |||
1668 | if self.axes: |
||
1669 | m = np.all( |
||
1670 | np.stack([pix_flat[i + 1] == t for i, t in enumerate(idx_img)]), |
||
1671 | axis=0, |
||
1672 | ) |
||
1673 | pix[0][s_img] = pix_flat[0][m] |
||
1674 | else: |
||
1675 | pix[0][s_img] = pix_flat[0] |
||
1676 | |||
1677 | for j in range(len(self.axes)): |
||
1678 | pix[j + 1][s_img] = idx_img[j] |
||
1679 | |||
1680 | pix = [p.T for p in pix] |
||
1681 | |||
1682 | if local: |
||
1683 | pix = self.global_to_local(pix) |
||
1684 | |||
1685 | if flat: |
||
1686 | pix = tuple([p[p != -1] for p in pix]) |
||
1687 | |||
1688 | return pix |
||
1689 | |||
1690 | def get_coord(self, idx=None, flat=False): |
||
1691 | pix = self.get_idx(idx=idx, flat=flat) |
||
1692 | coords = self.pix_to_coord(pix) |
||
1693 | cdict = OrderedDict([("lon", coords[0]), ("lat", coords[1])]) |
||
1694 | |||
1695 | for i, axis in enumerate(self.axes): |
||
1696 | cdict[axis.name] = coords[i + 2] |
||
1697 | |||
1698 | return MapCoord.create(cdict, coordsys=self.coordsys) |
||
1699 | |||
1700 | def contains(self, coords): |
||
1701 | idx = self.coord_to_idx(coords) |
||
1702 | return np.all(np.stack([t != -1 for t in idx]), axis=0) |
||
1703 | |||
1704 | def get_skydirs(self): |
||
1705 | """Get the sky coordinates of all the pixels in this geometry. """ |
||
1706 | coords = self.get_coord() |
||
1707 | frame = "galactic" if self.coordsys == "GAL" else "icrs" |
||
1708 | return SkyCoord(coords[0], coords[1], unit="deg", frame=frame) |
||
1709 | |||
1710 | def solid_angle(self): |
||
1711 | """Solid angle array (`~astropy.units.Quantity` in ``sr``). |
||
1712 | |||
1713 | The array has the same dimensionality as ``map.nside`` |
||
1714 | since all pixels have the same solid angle. |
||
1715 | """ |
||
1716 | import healpy as hp |
||
1717 | |||
1718 | return Quantity(hp.nside2pixarea(self.nside), "sr") |
||
1719 | |||
1720 | def __repr__(self): |
||
1721 | str_ = self.__class__.__name__ |
||
1722 | str_ += "\n\n" |
||
1723 | axes = ["skycoord"] + [_.name for _ in self.axes] |
||
1724 | str_ += "\taxes : {}\n".format(", ".join(axes)) |
||
1725 | str_ += "\tshape : {}\n".format(self.data_shape[::-1]) |
||
1726 | str_ += "\tndim : {}\n".format(self.ndim) |
||
1727 | str_ += "\tnside : {nside[0]}\n".format(nside=self.nside) |
||
1728 | str_ += "\tnested : {}\n".format(self.nest) |
||
1729 | str_ += "\tcoordsys : {}\n".format(self.coordsys) |
||
1730 | str_ += "\tprojection : {}\n".format(self.projection) |
||
1731 | lon, lat = self.center_skydir.data.lon.deg, self.center_skydir.data.lat.deg |
||
1732 | str_ += "\tcenter : {lon:.1f} deg, {lat:.1f} deg\n".format(lon=lon, lat=lat) |
||
1733 | return str_ |
||
1734 | |||
1735 | def __eq__(self, other): |
||
1736 | """Test equality between two `~gammapy.maps.HpxGeom`""" |
||
1737 | if not isinstance(other, self.__class__): |
||
1738 | return NotImplemented |
||
1739 | |||
1740 | if self._sparse or other._sparse: |
||
1741 | return NotImplemented |
||
1742 | if self.is_allsky and other.is_allsky is False: |
||
1743 | return NotImplemented |
||
1744 | |||
1745 | # check overall shape and axes compatibility |
||
1746 | if self.data_shape != other.data_shape: |
||
1747 | return False |
||
1748 | |||
1749 | for axis, otheraxis in zip(self.axes, other.axes): |
||
1750 | if axis != otheraxis: |
||
1751 | return False |
||
1752 | |||
1753 | return ( |
||
1754 | self.nside == other.nside |
||
1755 | and self.coordsys == other.coordsys |
||
1756 | and self.order == other.order |
||
1757 | and self.nest == other.nest |
||
1758 | ) |
||
1759 | |||
1760 | def __ne__(self, other): |
||
1761 | return not self.__eq__(other) |
||
1762 | |||
1763 | |||
1764 | class HpxToWcsMapping(object): |
||
1765 | """Stores the indices need to convert from HEALPIX to WCS. |
||
1766 | |||
1767 | Parameters |
||
1768 | ---------- |
||
1769 | hpx : `~HpxGeom` |
||
1770 | HEALPix geometry object. |
||
1771 | wcs : `~gammapy.maps.WcsGeom` |
||
1772 | WCS geometry object. |
||
1773 | """ |
||
1774 | |||
1775 | def __init__(self, hpx, wcs, ipix, mult_val, npix): |
||
1776 | self._hpx = hpx |
||
1777 | self._wcs = wcs |
||
1778 | self._ipix = ipix |
||
1779 | self._mult_val = mult_val |
||
1780 | self._npix = npix |
||
1781 | self._lmap = self._hpx[self._ipix] |
||
1782 | self._valid = self._lmap >= 0 |
||
1783 | |||
1784 | @property |
||
1785 | def hpx(self): |
||
1786 | """The HEALPIX projection""" |
||
1787 | return self._hpx |
||
1788 | |||
1789 | @property |
||
1790 | def wcs(self): |
||
1791 | """The WCS projection""" |
||
1792 | return self._wcs |
||
1793 | |||
1794 | @property |
||
1795 | def ipix(self): |
||
1796 | """An array(nx,ny) of the global HEALPIX pixel indices for each WCS pixel""" |
||
1797 | return self._ipix |
||
1798 | |||
1799 | @property |
||
1800 | def mult_val(self): |
||
1801 | """An array(nx,ny) of 1/number of WCS pixels pointing at each HEALPIX pixel""" |
||
1802 | return self._mult_val |
||
1803 | |||
1804 | @property |
||
1805 | def npix(self): |
||
1806 | """A tuple(nx,ny) of the shape of the WCS grid""" |
||
1807 | return self._npix |
||
1808 | |||
1809 | @property |
||
1810 | def lmap(self): |
||
1811 | """An array(nx,ny) giving the mapping of the local HEALPIX pixel |
||
1812 | indices for each WCS pixel""" |
||
1813 | return self._lmap |
||
1814 | |||
1815 | @property |
||
1816 | def valid(self): |
||
1817 | """An array(nx,ny) of bools giving if each WCS pixel in inside the |
||
1818 | HEALPIX region""" |
||
1819 | return self._valid |
||
1820 | |||
1821 | @classmethod |
||
1822 | def create(cls, hpx, wcs): |
||
1823 | """Create an object that maps pixels from HEALPix geometry ``hpx`` to |
||
1824 | WCS geometry ``wcs``. |
||
1825 | |||
1826 | Parameters |
||
1827 | ---------- |
||
1828 | hpx : `~HpxGeom` |
||
1829 | HEALPix geometry object. |
||
1830 | wcs : `~gammapy.maps.WcsGeom` |
||
1831 | WCS geometry object. |
||
1832 | |||
1833 | Returns |
||
1834 | ------- |
||
1835 | hpx2wcs : `~HpxToWcsMapping` |
||
1836 | |||
1837 | """ |
||
1838 | ipix, mult_val, npix = make_hpx_to_wcs_mapping(hpx, wcs) |
||
1839 | return cls(hpx, wcs, ipix, mult_val, npix) |
||
1840 | |||
1841 | def fill_wcs_map_from_hpx_data( |
||
1842 | self, hpx_data, wcs_data, normalize=True, fill_nan=True |
||
1843 | ): |
||
1844 | """Fill the WCS map from the hpx data using the pre-calculated mappings. |
||
1845 | |||
1846 | Parameters |
||
1847 | ---------- |
||
1848 | hpx_data : `~numpy.ndarray` |
||
1849 | The input HEALPIX data |
||
1850 | wcs_data : `~numpy.ndarray` |
||
1851 | The data array being filled |
||
1852 | normalize : bool |
||
1853 | True -> preserve integral by splitting HEALPIX values between bins |
||
1854 | fill_nan : bool |
||
1855 | Fill pixels outside the HPX geometry with NaN. |
||
1856 | """ |
||
1857 | # FIXME: Do we want to flatten mapping arrays? |
||
1858 | |||
1859 | shape = tuple([t.flat[0] for t in self._npix]) |
||
1860 | if self._valid.ndim != 1: |
||
1861 | shape = hpx_data.shape[:-1] + shape |
||
1862 | |||
1863 | valid = np.where(self._valid.reshape(shape)) |
||
1864 | lmap = self._lmap[self._valid] |
||
1865 | mult_val = self._mult_val[self._valid] |
||
1866 | |||
1867 | wcs_slice = [slice(None) for _ in range(wcs_data.ndim - 2)] |
||
1868 | wcs_slice = tuple(wcs_slice + list(valid)[::-1][:2]) |
||
1869 | |||
1870 | hpx_slice = [slice(None) for _ in range(wcs_data.ndim - 2)] |
||
1871 | hpx_slice = tuple(hpx_slice + [lmap]) |
||
1872 | |||
1873 | if normalize: |
||
1874 | wcs_data[wcs_slice] = mult_val * hpx_data[hpx_slice] |
||
1875 | else: |
||
1876 | wcs_data[wcs_slice] = hpx_data[hpx_slice] |
||
1877 | |||
1878 | if fill_nan: |
||
1879 | valid = np.swapaxes(self._valid.reshape(shape), -1, -2) |
||
1880 | valid = valid * np.ones_like(wcs_data, dtype=bool) |
||
1881 | wcs_data[~valid] = np.nan |
||
1882 | |||
1883 | return wcs_data |
||
1884 | |||
1885 | def make_wcs_data_from_hpx_data(self, hpx_data, wcs, normalize=True): |
||
1886 | """Create and fill a WCS map from the HEALPIX data using the pre-calculated mappings. |
||
1887 | |||
1888 | Parameters |
||
1889 | ---------- |
||
1890 | hpx_data : TODO |
||
1891 | The input HEALPIX data |
||
1892 | wcs : TODO |
||
1893 | The WCS object |
||
1894 | normalize : bool |
||
1895 | True -> preserve integral by splitting HEALPIX values between bins |
||
1896 | """ |
||
1897 | wcs_data = np.zeros(wcs.npix) |
||
1898 | self.fill_wcs_map_from_hpx_data(hpx_data, wcs_data, normalize) |
||
1899 | return wcs_data |
||
1900 |