1
|
|
|
import copy |
2
|
|
|
from functools import lru_cache |
3
|
|
|
import numpy as np |
4
|
|
|
from astropy import units as u |
5
|
|
|
from astropy.coordinates import Angle, SkyCoord |
6
|
|
|
from astropy.io import fits |
7
|
|
|
from astropy.table import QTable, Table |
8
|
|
|
from astropy.utils import lazyproperty |
9
|
|
|
from astropy.visualization.wcsaxes import WCSAxes |
10
|
|
|
from astropy.wcs.utils import ( |
11
|
|
|
proj_plane_pixel_area, |
12
|
|
|
proj_plane_pixel_scales, |
13
|
|
|
wcs_to_celestial_frame, |
14
|
|
|
) |
15
|
|
|
from regions import ( |
16
|
|
|
CompoundSkyRegion, |
17
|
|
|
PixCoord, |
18
|
|
|
PointSkyRegion, |
19
|
|
|
RectanglePixelRegion, |
20
|
|
|
RectangleSkyRegion, |
21
|
|
|
Regions, |
22
|
|
|
SkyRegion, |
23
|
|
|
) |
24
|
|
|
import matplotlib.pyplot as plt |
25
|
|
|
from gammapy.utils.regions import ( |
26
|
|
|
compound_region_center, |
27
|
|
|
compound_region_to_regions, |
28
|
|
|
regions_to_compound_region, |
29
|
|
|
) |
30
|
|
|
from gammapy.visualization.utils import ARTIST_TO_LINE_PROPERTIES |
31
|
|
|
from ..axes import MapAxes |
32
|
|
|
from ..coord import MapCoord |
33
|
|
|
from ..core import Map |
34
|
|
|
from ..geom import Geom, pix_tuple_to_idx |
35
|
|
|
from ..utils import _check_width |
36
|
|
|
from ..wcs import WcsGeom |
37
|
|
|
|
38
|
|
|
__all__ = ["RegionGeom"] |
39
|
|
|
|
40
|
|
|
|
41
|
|
|
class RegionGeom(Geom): |
42
|
|
|
"""Map geometry representing a region on the sky. |
43
|
|
|
The spatial component of the geometry is made up of a |
44
|
|
|
single pixel with an arbitrary shape and size. It can |
45
|
|
|
also have any number of non-spatial dimensions. This class |
46
|
|
|
represents the geometry for the `RegionNDMap` maps. |
47
|
|
|
|
48
|
|
|
Parameters |
49
|
|
|
---------- |
50
|
|
|
region : `~regions.SkyRegion` |
51
|
|
|
Region object. |
52
|
|
|
axes : list of `MapAxis` |
53
|
|
|
Non-spatial data axes. |
54
|
|
|
wcs : `~astropy.wcs.WCS` |
55
|
|
|
Optional wcs object to project the region if needed. |
56
|
|
|
binsz_wcs : `float` |
57
|
|
|
Angular bin size of the underlying `~WcsGeom` used to evaluate |
58
|
|
|
quantities in the region. Default size is 0.01 deg. This default |
59
|
|
|
value is adequate for the majority of use cases. If a wcs object |
60
|
|
|
is provided, the input of binsz_wcs is overridden. |
61
|
|
|
""" |
62
|
|
|
|
63
|
|
|
is_regular = True |
64
|
|
|
is_allsky = False |
65
|
|
|
is_hpx = False |
66
|
|
|
is_region = True |
67
|
|
|
|
68
|
|
|
_slice_spatial_axes = slice(0, 2) |
69
|
|
|
_slice_non_spatial_axes = slice(2, None) |
70
|
|
|
projection = "TAN" |
71
|
|
|
|
72
|
|
|
def __init__(self, region, axes=None, wcs=None, binsz_wcs="0.1 deg"): |
73
|
|
|
self._region = region |
74
|
|
|
self._axes = MapAxes.from_default(axes, n_spatial_axes=2) |
75
|
|
|
self._binsz_wcs = u.Quantity(binsz_wcs) |
76
|
|
|
|
77
|
|
|
if wcs is None and region is not None: |
78
|
|
|
if isinstance(region, CompoundSkyRegion): |
79
|
|
|
self._center = compound_region_center(region) |
80
|
|
|
else: |
81
|
|
|
self._center = region.center |
82
|
|
|
|
83
|
|
|
wcs = WcsGeom.create( |
84
|
|
|
binsz=binsz_wcs, |
85
|
|
|
skydir=self._center, |
86
|
|
|
proj=self.projection, |
87
|
|
|
frame=self._center.frame.name, |
88
|
|
|
).wcs |
89
|
|
|
|
90
|
|
|
self._wcs = wcs |
91
|
|
|
self.ndim = len(self.data_shape) |
92
|
|
|
|
93
|
|
|
# define cached methods |
94
|
|
|
self.get_wcs_coord_and_weights = lru_cache()(self.get_wcs_coord_and_weights) |
95
|
|
|
|
96
|
|
|
def __setstate__(self, state): |
97
|
|
|
for key, value in state.items(): |
98
|
|
|
if key in ["get_wcs_coord_and_weights"]: |
99
|
|
|
state[key] = lru_cache()(value) |
100
|
|
|
self.__dict__ = state |
101
|
|
|
|
102
|
|
|
@property |
103
|
|
|
def frame(self): |
104
|
|
|
"""Coordinate system, either Galactic ("galactic") or Equatorial |
105
|
|
|
("icrs").""" |
106
|
|
|
if self.region is None: |
107
|
|
|
return "icrs" |
108
|
|
|
try: |
109
|
|
|
return self.region.center.frame.name |
110
|
|
|
except AttributeError: |
111
|
|
|
return wcs_to_celestial_frame(self.wcs).name |
112
|
|
|
|
113
|
|
|
@property |
114
|
|
|
def binsz_wcs(self): |
115
|
|
|
"""Angular bin size of the underlying `~WcsGeom` |
116
|
|
|
|
117
|
|
|
Returns |
118
|
|
|
------- |
119
|
|
|
binsz_wcs: `~astropy.coordinates.Angle` |
120
|
|
|
""" |
121
|
|
|
return Angle(proj_plane_pixel_scales(self.wcs), unit="deg") |
122
|
|
|
|
123
|
|
|
@lazyproperty |
124
|
|
|
def _rectangle_bbox(self): |
125
|
|
|
if self.region is None: |
126
|
|
|
raise ValueError("Region definition required.") |
127
|
|
|
|
128
|
|
|
regions = compound_region_to_regions(self.region) |
129
|
|
|
regions_pix = [_.to_pixel(self.wcs) for _ in regions] |
130
|
|
|
|
131
|
|
|
bbox = regions_pix[0].bounding_box |
132
|
|
|
|
133
|
|
|
for region_pix in regions_pix[1:]: |
134
|
|
|
bbox = bbox.union(region_pix.bounding_box) |
135
|
|
|
|
136
|
|
|
try: |
137
|
|
|
rectangle_pix = bbox.to_region() |
138
|
|
|
except ValueError: |
139
|
|
|
rectangle_pix = RectanglePixelRegion( |
140
|
|
|
center=PixCoord(*bbox.center[::-1]), width=1, height=1 |
141
|
|
|
) |
142
|
|
|
return rectangle_pix.to_sky(self.wcs) |
143
|
|
|
|
144
|
|
|
@property |
145
|
|
|
def width(self): |
146
|
|
|
"""Width of bounding box of the region. |
147
|
|
|
|
148
|
|
|
Returns |
149
|
|
|
------- |
150
|
|
|
width : `~astropy.units.Quantity` |
151
|
|
|
Dimensions of the region in both spatial dimensions. |
152
|
|
|
Units: ``deg`` |
153
|
|
|
""" |
154
|
|
|
rectangle = self._rectangle_bbox |
155
|
|
|
return u.Quantity([rectangle.width.to("deg"), rectangle.height.to("deg")]) |
156
|
|
|
|
157
|
|
|
@property |
158
|
|
|
def region(self): |
159
|
|
|
"""`~regions.SkyRegion` object that defines the spatial component |
160
|
|
|
of the region geometry""" |
161
|
|
|
return self._region |
162
|
|
|
|
163
|
|
|
@property |
164
|
|
|
def is_all_point_sky_regions(self): |
165
|
|
|
"""Whether regions are all point regions""" |
166
|
|
|
regions = compound_region_to_regions(self.region) |
167
|
|
|
return np.all([isinstance(_, PointSkyRegion) for _ in regions]) |
168
|
|
|
|
169
|
|
|
@property |
170
|
|
|
def axes(self): |
171
|
|
|
"""List of non-spatial axes.""" |
172
|
|
|
return self._axes |
173
|
|
|
|
174
|
|
|
@property |
175
|
|
|
def axes_names(self): |
176
|
|
|
"""All axes names""" |
177
|
|
|
return ["lon", "lat"] + self.axes.names |
178
|
|
|
|
179
|
|
|
@property |
180
|
|
|
def wcs(self): |
181
|
|
|
"""WCS projection object.""" |
182
|
|
|
return self._wcs |
183
|
|
|
|
184
|
|
|
@property |
185
|
|
|
def center_coord(self): |
186
|
|
|
"""(`astropy.coordinates.SkyCoord`)""" |
187
|
|
|
return self.pix_to_coord(self.center_pix) |
188
|
|
|
|
189
|
|
|
@property |
190
|
|
|
def center_pix(self): |
191
|
|
|
"""Pixel values corresponding to the center of the region""" |
192
|
|
|
return tuple((np.array(self.data_shape) - 1.0) / 2)[::-1] |
193
|
|
|
|
194
|
|
|
@lazyproperty |
195
|
|
|
def center_skydir(self): |
196
|
|
|
"""Sky coordinate of the center of the region""" |
197
|
|
|
if self.region is None: |
198
|
|
|
return SkyCoord(np.nan * u.deg, np.nan * u.deg) |
199
|
|
|
|
200
|
|
|
return self._rectangle_bbox.center |
201
|
|
|
|
202
|
|
|
@property |
203
|
|
|
def npix(self): |
204
|
|
|
"""Number of spatial pixels""" |
205
|
|
|
return (1, 1) |
206
|
|
|
|
207
|
|
|
def contains(self, coords): |
208
|
|
|
"""Check if a given map coordinate is contained in the region. |
209
|
|
|
Requires the `.region` attribute to be set. |
210
|
|
|
|
211
|
|
|
Parameters |
212
|
|
|
---------- |
213
|
|
|
coords : tuple, dict, `MapCoord` or `~astropy.coordinates.SkyCoord` |
214
|
|
|
Object containing coordinate arrays we wish to check for inclusion |
215
|
|
|
in the region. |
216
|
|
|
|
217
|
|
|
Returns |
218
|
|
|
------- |
219
|
|
|
mask : `~numpy.ndarray` |
220
|
|
|
Boolean Numpy array with the same shape as the input that indicates |
221
|
|
|
which coordinates are inside the region. |
222
|
|
|
""" |
223
|
|
|
if self.region is None: |
224
|
|
|
raise ValueError("Region definition required.") |
225
|
|
|
|
226
|
|
|
coords = MapCoord.create(coords, frame=self.frame, axis_names=self.axes.names) |
227
|
|
|
|
228
|
|
|
if isinstance(self.region, PointSkyRegion): |
229
|
|
|
# the point region is approximated by 1x1 pixel here |
230
|
|
|
region = RectangleSkyRegion( |
231
|
|
|
center=self.center_skydir, |
232
|
|
|
width=self.binsz_wcs[0], |
233
|
|
|
height=self.binsz_wcs[1], |
234
|
|
|
) |
235
|
|
|
else: |
236
|
|
|
region = self.region |
237
|
|
|
|
238
|
|
|
return region.contains(coords.skycoord, self.wcs) |
239
|
|
|
|
240
|
|
|
def contains_wcs_pix(self, pix): |
241
|
|
|
"""Check if a given wcs pixel coordinate is contained in the region. |
242
|
|
|
|
243
|
|
|
Parameters |
244
|
|
|
---------- |
245
|
|
|
pix : tuple |
246
|
|
|
Tuple of pixel coordinates. |
247
|
|
|
|
248
|
|
|
Returns |
249
|
|
|
------- |
250
|
|
|
containment : `~numpy.ndarray` |
251
|
|
|
Bool array. |
252
|
|
|
""" |
253
|
|
|
region_pix = self.region.to_pixel(self.wcs) |
254
|
|
|
return region_pix.contains(PixCoord(pix[0], pix[1])) |
255
|
|
|
|
256
|
|
|
def separation(self, position): |
257
|
|
|
"""Angular distance between the center of the region and the given position. |
258
|
|
|
|
259
|
|
|
Parameters |
260
|
|
|
---------- |
261
|
|
|
position : `astropy.coordinates.SkyCoord` |
262
|
|
|
Sky coordinate we want the angular distance to. |
263
|
|
|
|
264
|
|
|
Returns |
265
|
|
|
------- |
266
|
|
|
sep : `~astropy.coordinates.Angle` |
267
|
|
|
The on-sky separation between the given coordinate and the region center. |
268
|
|
|
""" |
269
|
|
|
return self.center_skydir.separation(position) |
270
|
|
|
|
271
|
|
|
@property |
272
|
|
|
def data_shape(self): |
273
|
|
|
"""Shape of the Numpy data array matching this geometry.""" |
274
|
|
|
return self._shape[::-1] |
275
|
|
|
|
276
|
|
|
@property |
277
|
|
|
def data_shape_axes(self): |
278
|
|
|
"""Shape of data of the non-spatial axes and unit spatial axes.""" |
279
|
|
|
return self.axes.shape[::-1] + (1, 1) |
280
|
|
|
|
281
|
|
|
@property |
282
|
|
|
def _shape(self): |
283
|
|
|
"""Number of bins in each dimension. |
284
|
|
|
The spatial dimension is always (1, 1), as a |
285
|
|
|
`RegionGeom` is not pixelized further |
286
|
|
|
""" |
287
|
|
|
return tuple((1, 1) + self.axes.shape) |
288
|
|
|
|
289
|
|
|
def get_coord(self, mode="center", frame=None, sparse=False, axis_name=None): |
290
|
|
|
"""Get map coordinates from the geometry. |
291
|
|
|
|
292
|
|
|
Parameters |
293
|
|
|
---------- |
294
|
|
|
mode : {'center', 'edges'} |
295
|
|
|
Get center or edge coordinates for the non-spatial axes. |
296
|
|
|
frame : str or `~astropy.coordinates.Frame` |
297
|
|
|
Coordinate frame |
298
|
|
|
sparse : bool |
299
|
|
|
Compute sparse coordinates |
300
|
|
|
axis_name : str |
301
|
|
|
If mode = "edges", the edges will be returned for this axis only. |
302
|
|
|
|
303
|
|
|
|
304
|
|
|
Returns |
305
|
|
|
------- |
306
|
|
|
coord : `~MapCoord` |
307
|
|
|
Map coordinate object. |
308
|
|
|
""" |
309
|
|
|
if mode == "edges" and axis_name is None: |
310
|
|
|
raise ValueError("Mode 'edges' requires axis name") |
311
|
|
|
|
312
|
|
|
coords = self.axes.get_coord(mode=mode, axis_name=axis_name) |
313
|
|
|
coords["skycoord"] = self.center_skydir.reshape((1, 1)) |
314
|
|
|
|
315
|
|
|
if frame is None: |
316
|
|
|
frame = self.frame |
317
|
|
|
|
318
|
|
|
return MapCoord.create(coords, frame=self.frame).to_frame(frame) |
319
|
|
|
|
320
|
|
|
def _pad_spatial(self, pad_width): |
321
|
|
|
raise NotImplementedError("Spatial padding of `RegionGeom` not supported") |
322
|
|
|
|
323
|
|
|
def crop(self): |
324
|
|
|
raise NotImplementedError("Cropping of `RegionGeom` not supported") |
325
|
|
|
|
326
|
|
|
def solid_angle(self): |
327
|
|
|
"""Get solid angle of the region. |
328
|
|
|
|
329
|
|
|
Returns |
330
|
|
|
------- |
331
|
|
|
angle : `~astropy.units.Quantity` |
332
|
|
|
Solid angle of the region. In sr. |
333
|
|
|
Units: ``sr`` |
334
|
|
|
""" |
335
|
|
|
if self.region is None: |
336
|
|
|
raise ValueError("Region definition required.") |
337
|
|
|
|
338
|
|
|
# compound regions do not implement area() |
339
|
|
|
# so we use the mask representation and estimate the area |
340
|
|
|
# from the pixels in the mask using oversampling |
341
|
|
|
if isinstance(self.region, CompoundSkyRegion): |
342
|
|
|
# oversample by a factor of ten |
343
|
|
|
oversampling = 10.0 |
344
|
|
|
wcs = self.to_binsz_wcs(self.binsz_wcs / oversampling).wcs |
345
|
|
|
pixel_region = self.region.to_pixel(wcs) |
346
|
|
|
mask = pixel_region.to_mask() |
347
|
|
|
area = np.count_nonzero(mask) / oversampling**2 |
348
|
|
|
else: |
349
|
|
|
# all other types of regions should implement area |
350
|
|
|
area = self.region.to_pixel(self.wcs).area |
351
|
|
|
|
352
|
|
|
solid_angle = area * proj_plane_pixel_area(self.wcs) * u.deg**2 |
353
|
|
|
return solid_angle.to("sr") |
354
|
|
|
|
355
|
|
|
def bin_volume(self): |
356
|
|
|
"""If the RegionGeom has a non-spatial axis, it |
357
|
|
|
returns the volume of the region. If not, it |
358
|
|
|
just returns the solid angle size. |
359
|
|
|
|
360
|
|
|
Returns |
361
|
|
|
------- |
362
|
|
|
volume : `~astropy.units.Quantity` |
363
|
|
|
Volume of the region. |
364
|
|
|
""" |
365
|
|
|
bin_volume = self.solid_angle() * np.ones(self.data_shape) |
366
|
|
|
|
367
|
|
|
for idx, ax in enumerate(self.axes): |
368
|
|
|
shape = self.ndim * [1] |
369
|
|
|
shape[-(idx + 3)] = -1 |
370
|
|
|
bin_volume = bin_volume * ax.bin_width.reshape(tuple(shape)) |
371
|
|
|
|
372
|
|
|
return bin_volume |
373
|
|
|
|
374
|
|
|
def to_wcs_geom(self, width_min=None): |
375
|
|
|
"""Get the minimal equivalent geometry |
376
|
|
|
which contains the region. |
377
|
|
|
|
378
|
|
|
Parameters |
379
|
|
|
---------- |
380
|
|
|
width_min : `~astropy.quantity.Quantity` |
381
|
|
|
Minimal width for the resulting geometry. |
382
|
|
|
Can be a single number or two, for |
383
|
|
|
different minimum widths in each spatial dimension. |
384
|
|
|
|
385
|
|
|
Returns |
386
|
|
|
------- |
387
|
|
|
wcs_geom : `~WcsGeom` |
388
|
|
|
A WCS geometry object. |
389
|
|
|
""" |
390
|
|
|
if width_min is not None: |
391
|
|
|
width = np.max( |
392
|
|
|
[self.width.to_value("deg"), _check_width(width_min)], axis=0 |
393
|
|
|
) |
394
|
|
|
else: |
395
|
|
|
width = self.width |
396
|
|
|
wcs_geom_region = WcsGeom(wcs=self.wcs, npix=self.wcs.array_shape) |
397
|
|
|
wcs_geom = wcs_geom_region.cutout(position=self.center_skydir, width=width) |
398
|
|
|
wcs_geom = wcs_geom.to_cube(self.axes) |
399
|
|
|
return wcs_geom |
400
|
|
|
|
401
|
|
|
def to_binsz_wcs(self, binsz): |
402
|
|
|
|
403
|
|
|
"""Change the bin size of the underlying WCS geometry. |
404
|
|
|
|
405
|
|
|
Parameters |
406
|
|
|
---------- |
407
|
|
|
binzs : float, string or `~astropy.quantity.Quantity` |
408
|
|
|
|
409
|
|
|
Returns |
410
|
|
|
------- |
411
|
|
|
region : `~RegionGeom` |
412
|
|
|
A RegionGeom with the same axes and region as the input, |
413
|
|
|
but different wcs pixelization. |
414
|
|
|
""" |
415
|
|
|
new_geom = RegionGeom(self.region, axes=self.axes, binsz_wcs=binsz) |
416
|
|
|
return new_geom |
417
|
|
|
|
418
|
|
|
def get_wcs_coord_and_weights(self, factor=10): |
419
|
|
|
"""Get the array of spatial coordinates and corresponding weights |
420
|
|
|
|
421
|
|
|
The coordinates are the center of a pixel that intersects the region and |
422
|
|
|
the weights that represent which fraction of the pixel is contained |
423
|
|
|
in the region. |
424
|
|
|
|
425
|
|
|
Parameters |
426
|
|
|
---------- |
427
|
|
|
factor : int |
428
|
|
|
Oversampling factor to compute the weights |
429
|
|
|
|
430
|
|
|
Returns |
431
|
|
|
------- |
432
|
|
|
region_coord : `~MapCoord` |
433
|
|
|
MapCoord object with the coordinates inside |
434
|
|
|
the region. |
435
|
|
|
weights : `~np.array` |
436
|
|
|
Weights representing the fraction of each pixel |
437
|
|
|
contained in the region. |
438
|
|
|
""" |
439
|
|
|
wcs_geom = self.to_wcs_geom() |
440
|
|
|
|
441
|
|
|
weights = wcs_geom.to_image().region_weights( |
442
|
|
|
regions=[self.region], oversampling_factor=factor |
443
|
|
|
) |
444
|
|
|
|
445
|
|
|
mask = weights.data > 0 |
446
|
|
|
weights = weights.data[mask] |
447
|
|
|
|
448
|
|
|
# Get coordinates |
449
|
|
|
coords = wcs_geom.get_coord(sparse=True).apply_mask(mask) |
450
|
|
|
return coords, weights |
451
|
|
|
|
452
|
|
|
def to_binsz(self, binsz): |
453
|
|
|
"""Returns self""" |
454
|
|
|
return self |
455
|
|
|
|
456
|
|
|
def to_cube(self, axes): |
457
|
|
|
"""Append non-spatial axes to create a higher-dimensional geometry. |
458
|
|
|
|
459
|
|
|
Returns |
460
|
|
|
------- |
461
|
|
|
region : `~RegionGeom` |
462
|
|
|
RegionGeom with the added axes. |
463
|
|
|
""" |
464
|
|
|
axes = copy.deepcopy(self.axes) + axes |
465
|
|
|
return self._init_copy(axes=axes) |
466
|
|
|
|
467
|
|
|
def to_image(self): |
468
|
|
|
"""Remove non-spatial axes to create a 2D region. |
469
|
|
|
|
470
|
|
|
Returns |
471
|
|
|
------- |
472
|
|
|
region : `~RegionGeom` |
473
|
|
|
RegionGeom without any non-spatial axes. |
474
|
|
|
""" |
475
|
|
|
return self._init_copy(axes=None) |
476
|
|
|
|
477
|
|
|
def upsample(self, factor, axis_name=None): |
478
|
|
|
"""Upsample a non-spatial dimension of the region by a given factor. |
479
|
|
|
|
480
|
|
|
Returns |
481
|
|
|
------- |
482
|
|
|
region : `~RegionGeom` |
483
|
|
|
RegionGeom with the upsampled axis. |
484
|
|
|
""" |
485
|
|
|
axes = self.axes.upsample(factor=factor, axis_name=axis_name) |
486
|
|
|
return self._init_copy(axes=axes) |
487
|
|
|
|
488
|
|
|
def downsample(self, factor, axis_name): |
489
|
|
|
"""Downsample a non-spatial dimension of the region by a given factor. |
490
|
|
|
|
491
|
|
|
Returns |
492
|
|
|
------- |
493
|
|
|
region : `~RegionGeom` |
494
|
|
|
RegionGeom with the downsampled axis. |
495
|
|
|
""" |
496
|
|
|
axes = self.axes.downsample(factor=factor, axis_name=axis_name) |
497
|
|
|
return self._init_copy(axes=axes) |
498
|
|
|
|
499
|
|
|
def pix_to_coord(self, pix): |
500
|
|
|
lon = np.where( |
501
|
|
|
(-0.5 < pix[0]) & (pix[0] < 0.5), |
502
|
|
|
self.center_skydir.data.lon, |
503
|
|
|
np.nan * u.deg, |
504
|
|
|
) |
505
|
|
|
lat = np.where( |
506
|
|
|
(-0.5 < pix[1]) & (pix[1] < 0.5), |
507
|
|
|
self.center_skydir.data.lat, |
508
|
|
|
np.nan * u.deg, |
509
|
|
|
) |
510
|
|
|
coords = (lon, lat) |
511
|
|
|
|
512
|
|
|
for p, ax in zip(pix[self._slice_non_spatial_axes], self.axes): |
513
|
|
|
coords += (ax.pix_to_coord(p),) |
514
|
|
|
|
515
|
|
|
return coords |
516
|
|
|
|
517
|
|
|
def pix_to_idx(self, pix, clip=False): |
518
|
|
|
idxs = list(pix_tuple_to_idx(pix)) |
519
|
|
|
|
520
|
|
|
for i, idx in enumerate(idxs[self._slice_non_spatial_axes]): |
521
|
|
|
if clip: |
522
|
|
|
np.clip(idx, 0, self.axes[i].nbin - 1, out=idx) |
523
|
|
|
else: |
524
|
|
|
np.putmask(idx, (idx < 0) | (idx >= self.axes[i].nbin), -1) |
525
|
|
|
|
526
|
|
|
return tuple(idxs) |
527
|
|
|
|
528
|
|
|
def coord_to_pix(self, coords): |
529
|
|
|
# inherited docstring |
530
|
|
|
if isinstance(coords, tuple) and len(coords) == len(self.axes): |
531
|
|
|
skydir = self.center_skydir.transform_to(self.frame) |
532
|
|
|
coords = (skydir.data.lon, skydir.data.lat) + coords |
533
|
|
|
elif isinstance(coords, dict): |
534
|
|
|
valid_keys = ["lon", "lat", "skycoord"] |
535
|
|
|
if not any([_ in coords for _ in valid_keys]): |
536
|
|
|
coords.setdefault("skycoord", self.center_skydir) |
537
|
|
|
|
538
|
|
|
coords = MapCoord.create(coords, frame=self.frame, axis_names=self.axes.names) |
539
|
|
|
|
540
|
|
|
if self.region is None: |
541
|
|
|
pix = (0, 0) |
542
|
|
|
else: |
543
|
|
|
# TODO: remove once fix is available in regions |
544
|
|
|
if isinstance(self.region, PointSkyRegion): |
545
|
|
|
point_region = self.region.to_pixel(self.wcs) |
546
|
|
|
point_region.meta["include"] = False |
547
|
|
|
pix_coord = PixCoord.from_sky(coords.skycoord, self.wcs) |
548
|
|
|
in_region = point_region.contains(pix_coord) |
549
|
|
|
else: |
550
|
|
|
in_region = self.region.contains(coords.skycoord, wcs=self.wcs) |
551
|
|
|
|
552
|
|
|
x = np.zeros(coords.skycoord.shape) |
553
|
|
|
x[~in_region] = np.nan |
554
|
|
|
|
555
|
|
|
y = np.zeros(coords.skycoord.shape) |
556
|
|
|
y[~in_region] = np.nan |
557
|
|
|
|
558
|
|
|
pix = (x, y) |
559
|
|
|
|
560
|
|
|
pix += self.axes.coord_to_pix(coords) |
561
|
|
|
return pix |
562
|
|
|
|
563
|
|
|
def get_idx(self): |
564
|
|
|
idxs = [np.arange(n, dtype=float) for n in self.data_shape[::-1]] |
565
|
|
|
return np.meshgrid(*idxs[::-1], indexing="ij")[::-1] |
566
|
|
|
|
567
|
|
|
def _make_bands_cols(self): |
568
|
|
|
return [] |
569
|
|
|
|
570
|
|
|
@classmethod |
571
|
|
|
def create(cls, region, **kwargs): |
572
|
|
|
"""Create region geometry. |
573
|
|
|
|
574
|
|
|
The input region can be passed in the form of a ds9 string and will be parsed |
575
|
|
|
internally by `~regions.Regions.parse`. See: |
576
|
|
|
|
577
|
|
|
* https://astropy-regions.readthedocs.io/en/stable/region_io.html |
578
|
|
|
* http://ds9.si.edu/doc/ref/region.html |
579
|
|
|
|
580
|
|
|
Parameters |
581
|
|
|
---------- |
582
|
|
|
region : str or `~regions.SkyRegion` |
583
|
|
|
Region definition |
584
|
|
|
**kwargs : dict |
585
|
|
|
Keyword arguments passed to `RegionGeom.__init__` |
586
|
|
|
|
587
|
|
|
Returns |
588
|
|
|
------- |
589
|
|
|
geom : `RegionGeom` |
590
|
|
|
Region geometry |
591
|
|
|
""" |
592
|
|
|
return cls.from_regions(regions=region, **kwargs) |
593
|
|
|
|
594
|
|
|
def __repr__(self): |
595
|
|
|
axes = ["lon", "lat"] + [_.name for _ in self.axes] |
596
|
|
|
try: |
597
|
|
|
frame = self.center_skydir.frame.name |
598
|
|
|
lon = self.center_skydir.data.lon.deg |
599
|
|
|
lat = self.center_skydir.data.lat.deg |
600
|
|
|
except AttributeError: |
601
|
|
|
frame, lon, lat = "", np.nan, np.nan |
602
|
|
|
|
603
|
|
|
return ( |
604
|
|
|
f"{self.__class__.__name__}\n\n" |
605
|
|
|
f"\tregion : {self.region.__class__.__name__}\n" |
606
|
|
|
f"\taxes : {axes}\n" |
607
|
|
|
f"\tshape : {self.data_shape[::-1]}\n" |
608
|
|
|
f"\tndim : {self.ndim}\n" |
609
|
|
|
f"\tframe : {frame}\n" |
610
|
|
|
f"\tcenter : {lon:.1f} deg, {lat:.1f} deg\n" |
611
|
|
|
) |
612
|
|
|
|
613
|
|
|
def is_allclose(self, other, rtol_axes=1e-6, atol_axes=1e-6): |
614
|
|
|
"""Compare two data IRFs for equivalency |
615
|
|
|
|
616
|
|
|
Parameters |
617
|
|
|
---------- |
618
|
|
|
other : `RegionGeom` |
619
|
|
|
Geom to compare against. |
620
|
|
|
rtol_axes : float |
621
|
|
|
Relative tolerance for the axes comparison. |
622
|
|
|
atol_axes : float |
623
|
|
|
Relative tolerance for the axes comparison. |
624
|
|
|
|
625
|
|
|
Returns |
626
|
|
|
------- |
627
|
|
|
is_allclose : bool |
628
|
|
|
Whether the geometry is all close. |
629
|
|
|
""" |
630
|
|
|
if not isinstance(other, self.__class__): |
631
|
|
|
return TypeError(f"Cannot compare {type(self)} and {type(other)}") |
632
|
|
|
|
633
|
|
|
if self.data_shape != other.data_shape: |
634
|
|
|
return False |
635
|
|
|
|
636
|
|
|
axes_eq = self.axes.is_allclose(other.axes, rtol=rtol_axes, atol=atol_axes) |
637
|
|
|
# TODO: compare regions based on masks... |
638
|
|
|
regions_eq = True |
639
|
|
|
return axes_eq and regions_eq |
640
|
|
|
|
641
|
|
|
def __eq__(self, other): |
642
|
|
|
if not isinstance(other, self.__class__): |
643
|
|
|
return False |
644
|
|
|
|
645
|
|
|
return self.is_allclose(other=other) |
646
|
|
|
|
647
|
|
|
def _to_region_table(self): |
648
|
|
|
"""Export region to a FITS region table.""" |
649
|
|
|
if self.region is None: |
650
|
|
|
raise ValueError("Region definition required.") |
651
|
|
|
|
652
|
|
|
region_list = compound_region_to_regions(self.region) |
653
|
|
|
|
654
|
|
|
pixel_region_list = [] |
655
|
|
|
|
656
|
|
|
for reg in region_list: |
657
|
|
|
pixel_region_list.append(reg.to_pixel(self.wcs)) |
658
|
|
|
|
659
|
|
|
table = Regions(pixel_region_list).serialize(format="fits") |
660
|
|
|
|
661
|
|
|
header = WcsGeom(wcs=self.wcs, npix=self.wcs.array_shape).to_header() |
662
|
|
|
table.meta.update(header) |
663
|
|
|
return table |
664
|
|
|
|
665
|
|
|
def to_hdulist(self, format="ogip", hdu_bands=None, hdu_region=None): |
666
|
|
|
"""Convert geom to hdulist |
667
|
|
|
|
668
|
|
|
Parameters |
669
|
|
|
---------- |
670
|
|
|
format : {"gadf", "ogip", "ogip-sherpa"} |
671
|
|
|
HDU format |
672
|
|
|
hdu : str |
673
|
|
|
Name of the HDU with the map data. |
674
|
|
|
|
675
|
|
|
Returns |
676
|
|
|
------- |
677
|
|
|
hdulist : `~astropy.io.fits.HDUList` |
678
|
|
|
HDU list |
679
|
|
|
|
680
|
|
|
""" |
681
|
|
|
if hdu_bands is None: |
682
|
|
|
hdu_bands = "HDU_BANDS" |
683
|
|
|
if hdu_region is None: |
684
|
|
|
hdu_region = "HDU_REGION" |
685
|
|
|
if format != "gadf": |
686
|
|
|
hdu_region = "REGION" |
687
|
|
|
|
688
|
|
|
hdulist = fits.HDUList() |
689
|
|
|
|
690
|
|
|
hdulist.append(self.axes.to_table_hdu(hdu_bands=hdu_bands, format=format)) |
691
|
|
|
|
692
|
|
|
# region HDU |
693
|
|
|
if self.region: |
694
|
|
|
region_table = self._to_region_table() |
695
|
|
|
|
696
|
|
|
region_hdu = fits.BinTableHDU(region_table, name=hdu_region) |
697
|
|
|
hdulist.append(region_hdu) |
698
|
|
|
|
699
|
|
|
return hdulist |
700
|
|
|
|
701
|
|
|
@classmethod |
702
|
|
|
def from_regions(cls, regions, **kwargs): |
703
|
|
|
"""Create region geom from list of regions |
704
|
|
|
|
705
|
|
|
The regions are combined with union to a compound region. |
706
|
|
|
|
707
|
|
|
Parameters |
708
|
|
|
---------- |
709
|
|
|
regions : list of `~regions.SkyRegion` or str |
710
|
|
|
Regions |
711
|
|
|
**kwargs: dict |
712
|
|
|
Keyword arguments forwarded to `RegionGeom` |
713
|
|
|
|
714
|
|
|
Returns |
715
|
|
|
------- |
716
|
|
|
geom : `RegionGeom` |
717
|
|
|
Region map geometry |
718
|
|
|
""" |
719
|
|
|
if isinstance(regions, str): |
720
|
|
|
regions = Regions.parse(data=regions, format="ds9") |
721
|
|
|
elif isinstance(regions, SkyRegion): |
722
|
|
|
regions = [regions] |
723
|
|
|
elif isinstance(regions, SkyCoord): |
724
|
|
|
regions = [PointSkyRegion(center=regions)] |
725
|
|
|
|
726
|
|
|
if regions: |
727
|
|
|
regions = regions_to_compound_region(regions) |
728
|
|
|
|
729
|
|
|
return cls(region=regions, **kwargs) |
730
|
|
|
|
731
|
|
|
@classmethod |
732
|
|
|
def from_hdulist(cls, hdulist, format="ogip", hdu=None): |
733
|
|
|
"""Read region table and convert it to region list. |
734
|
|
|
|
735
|
|
|
Parameters |
736
|
|
|
---------- |
737
|
|
|
hdulist : `~astropy.io.fits.HDUList` |
738
|
|
|
HDU list |
739
|
|
|
format : {"ogip", "ogip-arf", "gadf"} |
740
|
|
|
HDU format |
741
|
|
|
|
742
|
|
|
Returns |
743
|
|
|
------- |
744
|
|
|
geom : `RegionGeom` |
745
|
|
|
Region map geometry |
746
|
|
|
|
747
|
|
|
""" |
748
|
|
|
region_hdu = "REGION" |
749
|
|
|
|
750
|
|
|
if format == "gadf" and hdu: |
751
|
|
|
region_hdu = hdu + "_" + region_hdu |
752
|
|
|
|
753
|
|
|
if region_hdu in hdulist: |
754
|
|
|
try: |
755
|
|
|
region_table = QTable.read(hdulist[region_hdu]) |
756
|
|
|
regions_pix = Regions.parse(data=region_table, format="fits") |
757
|
|
|
except TypeError: |
758
|
|
|
# TODO: this is needed to support regions=0.5 |
759
|
|
|
region_table = Table.read(hdulist[region_hdu]) |
760
|
|
|
regions_pix = Regions.parse(data=region_table, format="fits") |
761
|
|
|
|
762
|
|
|
wcs = WcsGeom.from_header(region_table.meta).wcs |
763
|
|
|
regions = [] |
764
|
|
|
|
765
|
|
|
for region_pix in regions_pix: |
766
|
|
|
# TODO: remove workaround once regions issue with fits serialization is sorted out |
767
|
|
|
# see https://github.com/astropy/regions/issues/400 |
768
|
|
|
region_pix.meta["include"] = True |
769
|
|
|
regions.append(region_pix.to_sky(wcs)) |
770
|
|
|
|
771
|
|
|
region = regions_to_compound_region(regions) |
772
|
|
|
else: |
773
|
|
|
region, wcs = None, None |
774
|
|
|
|
775
|
|
|
if format == "ogip": |
776
|
|
|
hdu_bands = "EBOUNDS" |
777
|
|
|
elif format == "ogip-arf": |
778
|
|
|
hdu_bands = "SPECRESP" |
779
|
|
|
elif format == "gadf": |
780
|
|
|
hdu_bands = hdu + "_BANDS" |
781
|
|
|
else: |
782
|
|
|
raise ValueError(f"Unknown format {format}") |
783
|
|
|
|
784
|
|
|
axes = MapAxes.from_table_hdu(hdulist[hdu_bands], format=format) |
785
|
|
|
return cls(region=region, wcs=wcs, axes=axes) |
786
|
|
|
|
787
|
|
|
def union(self, other): |
788
|
|
|
"""Stack a RegionGeom by making the union""" |
789
|
|
|
if not self == other: |
790
|
|
|
raise ValueError("Can only make union if extra axes are equivalent.") |
791
|
|
|
if other.region: |
792
|
|
|
if self.region: |
793
|
|
|
self._region = self.region.union(other.region) |
794
|
|
|
else: |
795
|
|
|
self._region = other.region |
796
|
|
|
|
797
|
|
|
def plot_region(self, ax=None, kwargs_point=None, path_effect=None, **kwargs): |
798
|
|
|
"""Plot region in the sky. |
799
|
|
|
|
800
|
|
|
Parameters |
801
|
|
|
---------- |
802
|
|
|
ax : `~astropy.visualization.WCSAxes` |
803
|
|
|
Axes to plot on. If no axes are given, |
804
|
|
|
the region is shown using the minimal |
805
|
|
|
equivalent WCS geometry. |
806
|
|
|
kwargs_point : dict |
807
|
|
|
Keyword arguments passed to `~matplotlib.lines.Line2D` for plotting |
808
|
|
|
of point sources |
809
|
|
|
path_effect : `~matplotlib.patheffects.PathEffect` |
810
|
|
|
Path effect applied to artists and lines. |
811
|
|
|
**kwargs : dict |
812
|
|
|
Keyword arguments forwarded to `~regions.PixelRegion.as_artist` |
813
|
|
|
|
814
|
|
|
Returns |
815
|
|
|
------- |
816
|
|
|
ax : `~astropy.visualization.WCSAxes` |
817
|
|
|
Axes to plot on. |
818
|
|
|
""" |
819
|
|
|
kwargs_point = kwargs_point or {} |
820
|
|
|
|
821
|
|
|
if ax is None: |
822
|
|
|
ax = plt.gca() |
823
|
|
|
|
824
|
|
|
if not isinstance(ax, WCSAxes): |
825
|
|
|
ax.remove() |
826
|
|
|
wcs_geom = self.to_wcs_geom() |
827
|
|
|
m = Map.from_geom(geom=wcs_geom.to_image()) |
828
|
|
|
ax = m.plot(add_cbar=False, vmin=-1, vmax=0) |
829
|
|
|
|
830
|
|
|
kwargs.setdefault("facecolor", "None") |
831
|
|
|
kwargs.setdefault("edgecolor", "tab:blue") |
832
|
|
|
kwargs_point.setdefault("marker", "*") |
833
|
|
|
|
834
|
|
|
for key, value in kwargs.items(): |
835
|
|
|
key_point = ARTIST_TO_LINE_PROPERTIES.get(key, None) |
836
|
|
|
if key_point: |
837
|
|
|
kwargs_point[key_point] = value |
838
|
|
|
|
839
|
|
|
for region in compound_region_to_regions(self.region): |
840
|
|
|
region_pix = region.to_pixel(wcs=ax.wcs) |
841
|
|
|
|
842
|
|
|
if isinstance(region, PointSkyRegion): |
843
|
|
|
artist = region_pix.as_artist(**kwargs_point) |
844
|
|
|
else: |
845
|
|
|
artist = region_pix.as_artist(**kwargs) |
846
|
|
|
|
847
|
|
|
if path_effect: |
848
|
|
|
artist.add_path_effect(path_effect) |
849
|
|
|
|
850
|
|
|
ax.add_artist(artist) |
851
|
|
|
|
852
|
|
|
return ax |
853
|
|
|
|