# Licensed under a 3-clause BSD style license - see LICENSE.rst"""Tools to create profiles (i.e. 1D "slices" from 2D images)."""importnumpyasnpimportscipy.ndimagefromastropyimportunitsasufromastropy.convolutionimportBox1DKernel,Gaussian1DKernelfromastropy.coordinatesimportAnglefromastropy.tableimportTableimportmatplotlib.pyplotaspltfromgammapy.maps.axesimportUNIT_STRING_FORMATfrom.coreimportEstimator__all__=["ImageProfile","ImageProfileEstimator"]# TODO: implement measuring profile along arbitrary directions# TODO: think better about error handling. e.g. MC based methods
[docs]classImageProfileEstimator(Estimator):"""Estimate profile from image. Parameters ---------- x_edges : `~astropy.coordinates.Angle`, optional Coordinate edges to define a custom measurement grid. method : {'sum', 'mean'} Compute sum or mean within profile bins. Default is 'sum'. axis : {'lon', 'lat', 'radial'} Along which axis to estimate the profile. Default is 'lon'. center : `~astropy.coordinates.SkyCoord`, optional Center coordinate for the radial profile option. Examples -------- This example shows how to compute a counts profile for the Fermi galactic center region:: import matplotlib.pyplot as plt from gammapy.estimators import ImageProfileEstimator from gammapy.maps import Map from astropy import units as u # load example data filename = '$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-counts.fits.gz' fermi_cts = Map.read(filename) # set up profile estimator and run p = ImageProfileEstimator(axis='lon', method='sum') profile = p.run(fermi_cts) # smooth profile and plot smoothed = profile.smooth(kernel='gauss') smoothed.peek() plt.show() """tag="ImageProfileEstimator"def__init__(self,x_edges=None,method="sum",axis="lon",center=None):ifmethodnotin["sum","mean"]:raiseValueError("Not a valid method, choose either 'sum' or 'mean'")ifaxisnotin["lon","lat","radial"]:raiseValueError("Not a valid axis, choose either 'lon', 'lat' or 'radial'")ifaxis=="radial"andcenterisNone:raiseValueError("Please provide center coordinate for radial profiles")self._x_edges=x_edgesself.parameters={"method":method,"axis":axis,"center":center}def_get_x_edges(self,image):ifself._x_edgesisnotNone:returnself._x_edgesp=self.parameterscoordinates=image.geom.get_coord(mode="edges").skycoordifp["axis"]=="lat":x_edges=coordinates[:,0].data.latelifp["axis"]=="lon":lon=coordinates[0,:].data.lonx_edges=lon.wrap_at("180d")elifp["axis"]=="radial":rad_step=image.geom.pixel_scales.mean()corners=[0,0,-1,-1],[0,-1,0,-1]rad_max=coordinates[corners].separation(p["center"]).max()x_edges=Angle(np.arange(0,rad_max.deg,rad_step.deg),unit="deg")returnx_edgesdef_estimate_profile(self,image,image_err,mask):p=self.parameterslabels=self._label_image(image,mask)profile_err=Noneindex=np.arange(1,len(self._get_x_edges(image)))ifp["method"]=="sum":profile=scipy.ndimage.sum(image.data,labels.data,index)ifimage.unit.is_equivalent("counts"):profile_err=np.sqrt(profile)elifimage_err:# gaussian error propagationerr_sum=scipy.ndimage.sum(image_err.data**2,labels.data,index)profile_err=np.sqrt(err_sum)elifp["method"]=="mean":# gaussian error propagationprofile=scipy.ndimage.mean(image.data,labels.data,index)ifimage_err:N=scipy.ndimage.sum(~np.isnan(image_err.data),labels.data,index)err_sum=scipy.ndimage.sum(image_err.data**2,labels.data,index)profile_err=np.sqrt(err_sum)/Nreturnprofile,profile_errdef_label_image(self,image,mask=None):p=self.parameterscoordinates=image.geom.get_coord().skycoordx_edges=self._get_x_edges(image)ifp["axis"]=="lon":lon=coordinates.data.lon.wrap_at("180d")data=np.digitize(lon.degree,x_edges.deg)elifp["axis"]=="lat":lat=coordinates.data.latdata=np.digitize(lat.degree,x_edges.deg)elifp["axis"]=="radial":separation=coordinates.separation(p["center"])data=np.digitize(separation.degree,x_edges.deg)ifmaskisnotNone:# assign masked values to backgrounddata[mask.data]=0returnimage.copy(data=data)
[docs]defrun(self,image,image_err=None,mask=None):"""Run image profile estimator. Parameters ---------- image : `~gammapy.maps.Map` Input image to run profile estimator on. image_err : `~gammapy.maps.Map`, optional Input error image to run profile estimator on. Default is None. mask : `~gammapy.maps.Map` Optional mask to exclude regions from the measurement. Returns ------- profile : `ImageProfile` Result image profile object. """p=self.parametersifimage.unit.is_equivalent("count"):image_err=image.copy(data=np.sqrt(image.data))profile,profile_err=self._estimate_profile(image,image_err,mask)result=Table()x_edges=self._get_x_edges(image)result["x_min"]=x_edges[:-1]result["x_max"]=x_edges[1:]result["x_ref"]=(x_edges[:-1]+x_edges[1:])/2result["profile"]=profile*image.unitifprofile_errisnotNone:result["profile_err"]=profile_err*image.unitresult.meta["PROFILE_TYPE"]=p["axis"]returnImageProfile(result)
[docs]classImageProfile:"""Image profile class. The image profile data is stored in `~astropy.table.Table` object, with the following columns: * `x_ref` Coordinate bin center (required). * `x_min` Coordinate bin minimum (optional). * `x_max` Coordinate bin maximum (optional). * `profile` Image profile data (required). * `profile_err` Image profile data error (optional). Parameters ---------- table : `~astropy.table.Table` Table instance with the columns specified as above. """def__init__(self,table):self.table=table
[docs]defsmooth(self,kernel="box",radius="0.1 deg",**kwargs):r"""Smooth profile with error propagation. Smoothing is described by a convolution: .. math:: x_j = \sum_i x_{(j - i)} h_i Where :math:`h_i` are the coefficients of the convolution kernel. The corresponding error on :math:`x_j` is then estimated using Gaussian error propagation, neglecting correlations between the individual :math:`x_{(j - i)}`: .. math:: \Delta x_j = \sqrt{\sum_i \Delta x^{2}_{(j - i)} h^{2}_i} Parameters ---------- kernel : {'gauss', 'box'} Kernel shape. Default is "box". radius : `~astropy.units.Quantity`, str or float Smoothing width given as quantity or float. If a float is given it is interpreted as smoothing width in pixels. If an (angular) quantity is given it is converted to pixels using `xref[1] - x_ref[0]`. Default is "0.1 deg". kwargs : dict, optional Keyword arguments passed to `~scipy.ndimage.uniform_filter` ('box') and `~scipy.ndimage.gaussian_filter` ('gauss'). Returns ------- profile : `ImageProfile` Smoothed image profile. """table=self.table.copy()profile=table["profile"]radius=u.Quantity(radius)radius=np.abs(radius/np.diff(self.x_ref))[0]width=2*radius.value+1ifkernel=="box":smoothed=scipy.ndimage.uniform_filter(profile.astype("float"),width,**kwargs)# renormalize dataiftable["profile"].unit.is_equivalent("count"):smoothed*=int(width)smoothed_err=np.sqrt(smoothed)elif"profile_err"intable.colnames:profile_err=table["profile_err"]# use gaussian error propagationbox=Box1DKernel(width)err_sum=scipy.ndimage.convolve(profile_err**2,box.array**2)smoothed_err=np.sqrt(err_sum)elifkernel=="gauss":smoothed=scipy.ndimage.gaussian_filter(profile.astype("float"),width,**kwargs)# use gaussian error propagationif"profile_err"intable.colnames:profile_err=table["profile_err"]gauss=Gaussian1DKernel(width)err_sum=scipy.ndimage.convolve(profile_err**2,gauss.array**2)smoothed_err=np.sqrt(err_sum)else:raiseValueError("Not valid kernel choose either 'box' or 'gauss'")table["profile"]=smoothed*self.table["profile"].unitif"profile_err"intable.colnames:table["profile_err"]=smoothed_err*self.table["profile"].unitreturnself.__class__(table)
@propertydefx_ref(self):"""Reference x coordinates."""returnself.table["x_ref"].quantity@propertydefx_min(self):"""Min. x coordinates."""returnself.table["x_min"].quantity@propertydefx_max(self):"""Max. x coordinates."""returnself.table["x_max"].quantity@propertydefprofile(self):"""Image profile quantity."""returnself.table["profile"].quantity@propertydefprofile_err(self):"""Image profile error quantity."""try:returnself.table["profile_err"].quantityexceptKeyError:returnNone
[docs]defpeek(self,figsize=(8,4.5),**kwargs):"""Show image profile and error. Parameters ---------- figsize : tuple Size of the figure. Default is (8, 4.5). **kwargs : dict, optional Keyword arguments passed to `ImageProfile.plot_profile()`. Returns ------- ax : `~matplotlib.axes.Axes` Axes object. """fig=plt.figure(figsize=figsize)ax=fig.add_axes([0.1,0.1,0.8,0.8])ax=self.plot(ax,**kwargs)if"profile_err"inself.table.colnames:ax=self.plot_err(ax,color=kwargs.get("c"))returnax
[docs]defnormalize(self,mode="peak"):"""Normalize profile to peak value or integral. Parameters ---------- mode : ['integral', 'peak'] Normalize image profile so that it integrates to unity ('integral') or the maximum value corresponds to one ('peak'). Default is "peak". Returns ------- profile : `ImageProfile` Normalized image profile. """table=self.table.copy()profile=self.table["profile"]ifmode=="peak":norm=np.nanmax(profile)elifmode=="integral":norm=np.nansum(profile)else:raiseValueError(f"Invalid normalization mode: {mode!r}")table["profile"]/=normif"profile_err"intable.colnames:table["profile_err"]/=normreturnself.__class__(table)