# Licensed under a 3-clause BSD style license - see LICENSE.rstimportnumpyasnpfromastropy.coordinatesimportSkyCoordfromregionsimportPolygonSkyRegionimportmatplotlibasmplimportmatplotlib.pyplotasplt
[docs]defcontainment_region(map_,fraction=0.393,apply_union=True):"""Find the iso-contours region corresponding to a given containment for a map of integral quantities with a flat geometry. Parameters ---------- map_ : `~gammapy.maps.WcsNDMap` Map of integral quantities. fraction : float Containment fraction. Default is 0.393. apply_union : bool It True return a compound region otherwise return a list of polygon regions. Default is True. Note that compound regions cannot be written in ds9 format but can always be saved with numpy.savez. Returns ------- regions : list of ~regions.PolygonSkyRegion` or `~regions.CompoundSkyRegion` regions from iso-contours matching containment fraction. """from.importWcsNDMapifnotisinstance(map_,WcsNDMap):raiseTypeError(f"This function is only supported for WcsNDMap, got {type(map_)} instead.")fmax=np.nanmax(map_.data)iffmax>0.0:ordered=np.sort(map_.data.flatten())[::-1]cumsum=np.nancumsum(ordered)ind=np.argmin(np.abs(cumsum/cumsum.max()-fraction))fval=ordered[ind]plt.ioff()fig=plt.figure()cs=plt.contour(map_.data.squeeze(),[fval])plt.close(fig)plt.ion()regions_pieces=[]try:paths_all=cs.get_paths()[0]starts=np.where(paths_all.codes==1)[0]stops=np.where(paths_all.codes==79)[0]+1paths=[]forstart,stopinzip(starts,stops):paths.append(mpl.path.Path(paths_all.vertices[start:stop],codes=paths_all.codes[start:stop],))exceptAttributeError:paths=cs.collections[0].get_paths()forppinpaths:vertices=[]forvinpp.vertices:v_coord=map_.geom.pix_to_coord(v)vertices.append([v_coord[0],v_coord[1]])vertices=SkyCoord(vertices,frame=map_.geom.frame)regions_pieces.append(PolygonSkyRegion(vertices))ifapply_union:regions_union=regions_pieces[0]forregioninregions_pieces[1:]:regions_union=regions_union.union(region)returnregions_unionelse:returnregions_pieceselse:raiseValueError("No positive values in the map.")
[docs]defcontainment_radius(map_,fraction=0.393,position=None):"""Compute containment radius from a position in a map with integral quantities and flat geometry. Parameters ---------- map_ : `~gammapy.maps.WcsNDMap` Map of integral quantities. fraction : float Containment fraction. Default is 0.393. position : `~astropy.coordinates.SkyCoord` Position from where the containment is computed. Default is the center of the Map. Returns ------- radius : `~astropy.coordinates.Angle` Minimal radius required to reach the given containement fraction. """from.importWcsNDMapifnotisinstance(map_,WcsNDMap):raiseTypeError(f"This function is only supported for WcsNDMap, got {type(map_)} instead.")ifpositionisNone:position=map_.geom.center_skydirfmax=np.nanmax(map_.data)iffmax>0.0:sep=map_.geom.separation(position).flatten()order=np.argsort(sep)ordered=map_.data.flatten()[order]cumsum=np.nancumsum(ordered)ind=np.argmin(np.abs(cumsum/cumsum.max()-fraction))else:raiseValueError("No positive values in the map.")returnsep[order][ind]