# Licensed under a 3-clause BSD style license - see LICENSE.rstimportcopyimportloggingfromfunctoolsimportlru_cacheimportnumpyasnpfromastropyimportunitsasufromastropy.coordinatesimportAngle,SkyCoordfromastropy.ioimportfitsfromastropy.tableimportQTable,Tablefromastropy.utilsimportlazypropertyfromastropy.visualization.wcsaxesimportWCSAxesfromastropy.wcs.utilsimport(proj_plane_pixel_area,proj_plane_pixel_scales,wcs_to_celestial_frame,)fromregionsimport(CompoundSkyRegion,PixCoord,PointSkyRegion,RectanglePixelRegion,Regions,SkyRegion,)importmatplotlib.pyplotaspltfromgammapy.utils.regionsimport(compound_region_center,compound_region_to_regions,regions_to_compound_region,)fromgammapy.visualization.utilsimportARTIST_TO_LINE_PROPERTIESfrom..axesimportMapAxesfrom..coordimportMapCoordfrom..coreimportMapfrom..geomimportGeom,pix_tuple_to_idxfrom..utilsimport_check_widthfrom..wcsimportWcsGeom__all__=["RegionGeom"]
[docs]classRegionGeom(Geom):"""Map geometry representing a region on the sky. The spatial component of the geometry is made up of a single pixel with an arbitrary shape and size. It can also have any number of non-spatial dimensions. This class represents the geometry for the `RegionNDMap` maps. Parameters ---------- region : `~regions.SkyRegion` Region object. axes : list of `MapAxis` Non-spatial data axes. wcs : `~astropy.wcs.WCS` Optional WCS object to project the region if needed. binsz_wcs : `~astropy.units.Quantity` Angular bin size of the underlying `~WcsGeom` used to evaluate quantities in the region. Default is "0.1 deg". This default value is adequate for the majority of use cases. If a WCS object is provided, the input of ``binsz_wcs`` is overridden. """is_regular=Trueis_allsky=Falseis_hpx=Falseis_region=True_slice_spatial_axes=slice(0,2)_slice_non_spatial_axes=slice(2,None)projection="TAN"def__init__(self,region,axes=None,wcs=None,binsz_wcs="0.1 deg"):self._region=regionself._axes=MapAxes.from_default(axes,n_spatial_axes=2)self._binsz_wcs=u.Quantity(binsz_wcs)ifwcsisNoneandregionisnotNone:ifisinstance(region,CompoundSkyRegion):self._center=compound_region_center(region)else:self._center=region.centerwcs=WcsGeom.create(binsz=binsz_wcs,skydir=self._center,proj=self.projection,frame=self._center.frame.name,).wcsself._wcs=wcs# TODO : can we get the width before defining the wcs ?wcs=WcsGeom.create(binsz=binsz_wcs,width=tuple(self.width),skydir=self._center,proj=self.projection,frame=self._center.frame.name,).wcsself._wcs=wcsself.ndim=len(self.data_shape)# define cached methodsself.get_wcs_coord_and_weights=lru_cache()(self.get_wcs_coord_and_weights)def__setstate__(self,state):forkey,valueinstate.items():ifkeyin["get_wcs_coord_and_weights"]:state[key]=lru_cache()(value)self.__dict__=state@propertydefframe(self):"""Coordinate system, either Galactic ("galactic") or Equatorial ("icrs")."""ifself.regionisNone:return"icrs"try:returnself.region.center.frame.nameexceptAttributeError:returnwcs_to_celestial_frame(self.wcs).name@propertydefbinsz_wcs(self):"""Angular bin size of the underlying `~WcsGeom`. Returns ------- binsz_wcs: `~astropy.coordinates.Angle` Angular bin size. """returnAngle(proj_plane_pixel_scales(self.wcs),unit="deg")@lazypropertydef_rectangle_bbox(self):ifself.regionisNone:raiseValueError("Region definition required.")regions=compound_region_to_regions(self.region)regions_pix=[_.to_pixel(self.wcs)for_inregions]bbox=regions_pix[0].bounding_boxforregion_pixinregions_pix[1:]:bbox=bbox.union(region_pix.bounding_box)try:rectangle_pix=bbox.to_region()exceptValueError:rectangle_pix=RectanglePixelRegion(center=PixCoord(*bbox.center[::-1]),width=1,height=1)returnrectangle_pix.to_sky(self.wcs)@propertydefwidth(self):"""Width of bounding box of the region. Returns ------- width : `~astropy.units.Quantity` Dimensions of the region in both spatial dimensions. Units: ``deg`` """rectangle=self._rectangle_bboxreturnu.Quantity([rectangle.width.to("deg"),rectangle.height.to("deg")])@propertydefregion(self):"""The spatial component of the region geometry as a `~regions.SkyRegion`."""returnself._region@propertydefis_all_point_sky_regions(self):"""Whether regions are all point regions."""regions=compound_region_to_regions(self.region)returnnp.all([isinstance(_,PointSkyRegion)for_inregions])@propertydefaxes(self):"""List of non-spatial axes."""returnself._axes@propertydefaxes_names(self):"""All axes names."""return["lon","lat"]+self.axes.names@propertydefwcs(self):"""WCS projection object."""returnself._wcs@propertydefcenter_coord(self):"""Center coordinate of the region as a `astropy.coordinates.SkyCoord`."""returnself.pix_to_coord(self.center_pix)@propertydefcenter_pix(self):"""Pixel values corresponding to the center of the region."""returntuple((np.array(self.data_shape)-1.0)/2)[::-1]@lazypropertydefcenter_skydir(self):"""Sky coordinate of the center of the region."""ifself.regionisNone:returnSkyCoord(np.nan*u.deg,np.nan*u.deg)returnself._rectangle_bbox.center@propertydefnpix(self):"""Number of spatial pixels."""return([1],[1])
[docs]defcontains(self,coords):"""Check if a given map coordinate is contained in the region. Requires the `.region` attribute to be set. For `PointSkyRegion` the method always returns True. Parameters ---------- coords : tuple, dict, `MapCoord` or `~astropy.coordinates.SkyCoord` Object containing coordinate arrays we wish to check for inclusion in the region. Returns ------- mask : `~numpy.ndarray` Boolean array with the same shape as the input that indicates which coordinates are inside the region. """ifself.regionisNone:raiseValueError("Region definition required.")coords=MapCoord.create(coords,frame=self.frame,axis_names=self.axes.names)ifself.is_all_point_sky_regions:returnnp.ones(coords.skycoord.shape,dtype=bool)returnself.region.contains(coords.skycoord,self.wcs)
[docs]defcontains_wcs_pix(self,pix):"""Check if a given WCS pixel coordinate is contained in the region. For `PointSkyRegion` the method always returns True. Parameters ---------- pix : tuple Tuple of pixel coordinates. Returns ------- containment : `~numpy.ndarray` Boolean array. """ifself.is_all_point_sky_regions:returnnp.ones(pix[0].shape,dtype=bool)region_pix=self.region.to_pixel(self.wcs)returnregion_pix.contains(PixCoord(pix[0],pix[1]))
[docs]defseparation(self,position):"""Angular distance between the center of the region and the given position. Parameters ---------- position : `astropy.coordinates.SkyCoord` Sky coordinate we want the angular distance to. Returns ------- sep : `~astropy.coordinates.Angle` The on-sky separation between the given coordinate and the region center. """returnself.center_skydir.separation(position)
@propertydefdata_shape(self):"""Shape of the `~numpy.ndarray` matching this geometry."""returnself._shape[::-1]@propertydefdata_shape_axes(self):"""Shape of data of the non-spatial axes and unit spatial axes."""returnself.axes.shape[::-1]+(1,1)@propertydef_shape(self):"""Number of bins in each dimension. The spatial dimension is always (1, 1), as a `RegionGeom` is not pixelized further. """returntuple((1,1)+self.axes.shape)
[docs]defget_coord(self,mode="center",frame=None,sparse=False,axis_name=None):"""Get map coordinates from the geometry. Parameters ---------- mode : {'center', 'edges'}, optional Get center or edge coordinates for the non-spatial axes. Default is "center". frame : str or `~astropy.coordinates.Frame`, optional Coordinate frame. Default is None. sparse : bool, optional Compute sparse coordinates. Default is False. axis_name : str, optional If mode = "edges", the edges will be returned for this axis only. Default is None. Returns ------- coord : `~MapCoord` Map coordinate object. """ifmode=="edges"andaxis_nameisNone:raiseValueError("Mode 'edges' requires axis name")coords=self.axes.get_coord(mode=mode,axis_name=axis_name)coords["skycoord"]=self.center_skydir.reshape((1,1))ifframeisNone:frame=self.framecoords=MapCoord.create(coords,frame=self.frame).to_frame(frame)ifnotsparse:coords=coords.broadcastedreturncoords
def_pad_spatial(self,pad_width):raiseNotImplementedError("Spatial padding of `RegionGeom` not supported")
[docs]defcrop(self):raiseNotImplementedError("Cropping of `RegionGeom` not supported")
[docs]defsolid_angle(self):"""Get solid angle of the region. Returns ------- angle : `~astropy.units.Quantity` Solid angle of the region in steradians. Units: ``sr`` """ifself.regionisNone:raiseValueError("Region definition required.")# compound regions do not implement area()# so we use the mask representation and estimate the area# from the pixels in the mask using oversamplingifisinstance(self.region,CompoundSkyRegion):# oversample by a factor of tenoversampling=10.0wcs=self.to_binsz_wcs(self.binsz_wcs/oversampling).wcspixel_region=self.region.to_pixel(wcs)mask=pixel_region.to_mask()area=np.count_nonzero(mask)/oversampling**2else:# all other types of regions should implement areaarea=self.region.to_pixel(self.wcs).areasolid_angle=area*proj_plane_pixel_area(self.wcs)*u.deg**2returnsolid_angle.to("sr")
[docs]defbin_volume(self):"""If the `RegionGeom` has a non-spatial axis, it returns the volume of the region. If not, it returns the solid angle size. Returns ------- volume : `~astropy.units.Quantity` Volume of the region. """bin_volume=self.solid_angle()*np.ones(self.data_shape)foridx,axinenumerate(self.axes):shape=self.ndim*[1]shape[-(idx+3)]=-1bin_volume=bin_volume*ax.bin_width.reshape(tuple(shape))returnbin_volume
[docs]defto_wcs_geom(self,width_min=None):"""Get the minimal equivalent geometry which contains the region. Parameters ---------- width_min : `~astropy.quantity.Quantity`, optional Minimum width for the resulting geometry. Can be a single number or two, for different minimum widths in each spatial dimension. Default is None. Returns ------- wcs_geom : `~WcsGeom` A WCS geometry object. """ifwidth_minisnotNone:width=np.max([self.width.to_value("deg"),_check_width(width_min)],axis=0)else:width=self.widthwcs_geom_region=WcsGeom(wcs=self.wcs)wcs_geom=wcs_geom_region.cutout(position=self.center_skydir,width=width)wcs_geom=wcs_geom.to_cube(self.axes)returnwcs_geom
[docs]defto_binsz_wcs(self,binsz):"""Change the bin size of the underlying WCS geometry. Parameters ---------- binsz : float, str or `~astropy.quantity.Quantity` Bin size. Returns ------- region : `~RegionGeom` Region geometry with the same axes and region as the input, but different WCS pixelization. """new_geom=RegionGeom(self.region,axes=self.axes,binsz_wcs=binsz)returnnew_geom
[docs]defget_wcs_coord_and_weights(self,factor=10):"""Get the array of spatial coordinates and corresponding weights. The coordinates are the center of a pixel that intersects the region and the weights that represent which fraction of the pixel is contained in the region. Parameters ---------- factor : int, optional Oversampling factor to compute the weights. Default is 10. Returns ------- region_coord : `~MapCoord` MapCoord object with the coordinates inside the region. weights : `~numpy.ndarray` Weights representing the fraction of each pixel contained in the region. """wcs_geom=self.to_wcs_geom()weights=wcs_geom.to_image().region_weights(regions=[self.region],oversampling_factor=factor)mask=weights.data>0weights=weights.data[mask]# Get coordinatescoords=wcs_geom.get_coord(sparse=True).apply_mask(mask)returncoords,weights
[docs]defto_cube(self,axes):"""Append non-spatial axes to create a higher-dimensional geometry. Returns ------- region : `~RegionGeom` Region geometry with the added axes. """axes=copy.deepcopy(self.axes)+axesreturnself._init_copy(region=self.region,wcs=self.wcs,axes=axes)
[docs]defto_image(self):"""Remove non-spatial axes to create a 2D region. Returns ------- region : `~RegionGeom` Region geometry without any non-spatial axes. """returnself._init_copy(region=self.region,wcs=self.wcs,axes=None)
[docs]defupsample(self,factor,axis_name=None):"""Upsample a non-spatial dimension of the region by a given factor. Returns ------- region : `~RegionGeom` Region geometry with the upsampled axis. """axes=self.axes.upsample(factor=factor,axis_name=axis_name)returnself._init_copy(region=self.region,wcs=self.wcs,axes=axes)
[docs]defdownsample(self,factor,axis_name):"""Downsample a non-spatial dimension of the region by a given factor. Returns ------- region : `~RegionGeom` Region geometry with the downsampled axis. """axes=self.axes.downsample(factor=factor,axis_name=axis_name)returnself._init_copy(region=self.region,wcs=self.wcs,axes=axes)
[docs]@classmethoddefcreate(cls,region,**kwargs):"""Create region geometry. The input region can be passed in the form of a ds9 string and will be parsed internally by `~regions.Regions.parse`. See: * https://astropy-regions.readthedocs.io/en/stable/region_io.html * http://ds9.si.edu/doc/ref/region.html Parameters ---------- region : str or `~regions.SkyRegion` Region definition. **kwargs : dict Keyword arguments passed to `RegionGeom.__init__`. Returns ------- geom : `RegionGeom` Region geometry. """returncls.from_regions(regions=region,**kwargs)
[docs]defis_allclose(self,other,rtol_axes=1e-6,atol_axes=1e-6):"""Compare two data IRFs for equivalency. Parameters ---------- other : `RegionGeom` Region geometry to compare against. rtol_axes : float, optional Relative tolerance for the axes comparison. Default is 1e-6. atol_axes : float, optional Relative tolerance for the axes comparison. Default is 1e-6. Returns ------- is_allclose : bool Whether the geometry is all close. """ifnotisinstance(other,self.__class__):returnTypeError(f"Cannot compare {type(self)} and {type(other)}")ifself.data_shape!=other.data_shape:returnFalseaxes_eq=self.axes.is_allclose(other.axes,rtol=rtol_axes,atol=atol_axes)# TODO: compare regions based on masks...regions_eq=Truereturnaxes_eqandregions_eq
def__eq__(self,other):ifnotisinstance(other,self.__class__):returnFalsereturnself.is_allclose(other=other)def_to_region_table(self):"""Export region to a FITS region table."""ifself.regionisNone:raiseValueError("Region definition required.")region_list=compound_region_to_regions(self.region)pixel_region_list=[]forreginregion_list:pixel_region_list.append(reg.to_pixel(self.wcs))table=Regions(pixel_region_list).serialize(format="fits")header=WcsGeom(wcs=self.wcs).to_header()table.meta.update(header)returntable
[docs]defto_hdulist(self,format="ogip",hdu_bands=None,hdu_region=None):"""Convert geometry to HDU list. Parameters ---------- format : {"ogip", "gadf", "ogip-sherpa"} HDU format. Default is "ogip". hdu_bands : str, optional Name or index of the HDU with the BANDS table. Default is None. hdu_region : str, optional Name or index of the HDU with the region table. Not used for the "gadf" format. Default is None. Returns ------- hdulist : `~astropy.io.fits.HDUList` HDU list. """ifhdu_bandsisNone:hdu_bands="HDU_BANDS"ifhdu_regionisNone:hdu_region="HDU_REGION"ifformat!="gadf":hdu_region="REGION"hdulist=fits.HDUList()hdulist.append(self.axes.to_table_hdu(hdu_bands=hdu_bands,format=format))# region HDUifself.region:region_table=self._to_region_table()region_hdu=fits.BinTableHDU(region_table,name=hdu_region)hdulist.append(region_hdu)returnhdulist
[docs]@classmethoddeffrom_regions(cls,regions,**kwargs):"""Create region geometry from list of regions. The regions are combined with union to a compound region. Parameters ---------- regions : list of `~regions.SkyRegion` or str Regions. **kwargs: dict Keyword arguments forwarded to `RegionGeom`. Returns ------- geom : `RegionGeom` Region map geometry. """ifisinstance(regions,str):regions=Regions.parse(data=regions,format="ds9")elifisinstance(regions,SkyRegion):regions=[regions]elifisinstance(regions,SkyCoord):regions=[PointSkyRegion(center=regions)]elifisinstance(regions,list)andlen(regions)==0:regions=Noneifregions:regions=regions_to_compound_region(regions)returncls(region=regions,**kwargs)
[docs]@classmethoddeffrom_hdulist(cls,hdulist,format="ogip",hdu=None):"""Read region table and convert it to region list. Parameters ---------- hdulist : `~astropy.io.fits.HDUList` HDU list. format : {"ogip", "ogip-arf", "gadf"} HDU format. Default is "ogip". hdu : str, optional Name of the HDU. Default is None. Returns ------- geom : `RegionGeom` Region map geometry. """region_hdu="REGION"ifformat=="gadf"andhdu:region_hdu=hdu+"_"+region_hduifregion_hduinhdulist:try:region_table=QTable.read(hdulist[region_hdu])regions_pix=Regions.parse(data=region_table,format="fits")exceptTypeError:region_table=Table.read(hdulist[region_hdu])regions_pix=Regions.parse(data=region_table,format="fits")wcs=WcsGeom.from_header(region_table.meta).wcsregions=[]forregion_pixinregions_pix:# TODO: remove workaround once regions issue with fits serialization is sorted out# see https://github.com/astropy/regions/issues/400region_pix.meta["include"]=Trueregions.append(region_pix.to_sky(wcs))region=regions_to_compound_region(regions)else:region,wcs=None,Noneifformat=="ogip":hdu_bands="EBOUNDS"elifformat=="ogip-arf":hdu_bands="SPECRESP"elifformat=="gadf":hdu_bands=hdu+"_BANDS"else:raiseValueError(f"Unknown format {format}")axes=MapAxes.from_table_hdu(hdulist[hdu_bands],format=format)returncls(region=region,wcs=wcs,axes=axes)
[docs]defunion(self,other):"""Stack a `RegionGeom` by making the union."""ifnotself==other:raiseValueError("Can only make union if extra axes are equivalent.")ifother.region:ifself.region:self._region=self.region.union(other.region)else:self._region=other.region
[docs]defplot_region(self,ax=None,kwargs_point=None,path_effect=None,**kwargs):"""Plot region in the sky. Parameters ---------- ax : `~astropy.visualization.WCSAxes`, optional Axes to plot on. If no axes are given, the region is shown using the minimal equivalent WCS geometry. Default is None. kwargs_point : dict, optional Keyword arguments passed to `~matplotlib.lines.Line2D` for plotting of point sources. Default is None. path_effect : `~matplotlib.patheffects.PathEffect`, optional Path effect applied to artists and lines. Default is None. **kwargs : dict Keyword arguments forwarded to `~regions.PixelRegion.as_artist`. Returns ------- ax : `~astropy.visualization.WCSAxes` Axes to plot on. """ifself.region:kwargs_point=kwargs_pointor{}ifaxisNone:ax=plt.gca()ifnotisinstance(ax,WCSAxes):ax.remove()wcs_geom=self.to_wcs_geom()m=Map.from_geom(geom=wcs_geom.to_image())ax=m.plot(add_cbar=False,vmin=-1,vmax=0)kwargs.setdefault("facecolor","None")kwargs.setdefault("edgecolor","tab:blue")kwargs_point.setdefault("marker","*")forkey,valueinkwargs.items():key_point=ARTIST_TO_LINE_PROPERTIES.get(key,None)ifkey_point:kwargs_point[key_point]=valueforregionincompound_region_to_regions(self.region):region_pix=region.to_pixel(wcs=ax.wcs)ifisinstance(region,PointSkyRegion):artist=region_pix.as_artist(**kwargs_point)else:artist=region_pix.as_artist(**kwargs)ifpath_effect:artist.add_path_effect(path_effect)ax.add_artist(artist)returnaxelse:logging.info("Region definition required.")