# Licensed under a 3-clause BSD style license - see LICENSE.rstimportloggingimportnumpyasnpimportscipy.interpolateimportscipy.ndimageasndiimportscipy.signalimportastropy.unitsasufromastropy.convolutionimportTophat2DKernelfromastropy.coordinatesimportSkyCoordfromastropy.ioimportfitsfromastropy.nddataimportblock_reducefromregionsimportPixCoord,PointPixelRegion,PointSkyRegion,SkyRegionimportmatplotlib.pyplotaspltfromgammapy.utils.interpolationimportScaledRegularGridInterpolatorfromgammapy.utils.unitsimportunit_from_fits_image_hdufrom..geomimportpix_tuple_to_idxfrom..utilsimportINVALID_INDEXfrom.coreimportWcsMapfrom.geomimportWcsGeom__all__=["WcsNDMap"]log=logging.getLogger(__name__)
[docs]classWcsNDMap(WcsMap):"""WCS map with any number of non-spatial dimensions. This class uses an ND numpy array to store map values. For maps with non-spatial dimensions and variable pixel size it will allocate an array with dimensions commensurate with the largest image plane. Parameters ---------- geom : `~gammapy.maps.WcsGeom` WCS geometry object. data : `~numpy.ndarray` Data array. If none then an empty array will be allocated. dtype : str, optional Data type, default is float32 meta : `dict` Dictionary to store meta data. unit : str or `~astropy.units.Unit` The map unit """def__init__(self,geom,data=None,dtype="float32",meta=None,unit=""):# TODO: Figure out how to mask pixels for integer data typesdata_shape=geom.data_shapeifdataisNone:data=self._make_default_data(geom,data_shape,dtype)super().__init__(geom,data,meta,unit)@staticmethoddef_make_default_data(geom,shape_np,dtype):# Check whether corners of each image plane are validdata=np.zeros(shape_np,dtype=dtype)ifnotgeom.is_regularorgeom.is_allsky:coords=geom.get_coord()is_nan=np.isnan(coords.lon)data[is_nan]=np.nanreturndata
[docs]@classmethoddeffrom_hdu(cls,hdu,hdu_bands=None,format=None):"""Make a WcsNDMap object from a FITS HDU. Parameters ---------- hdu : `~astropy.io.fits.BinTableHDU` or `~astropy.io.fits.ImageHDU` The map FITS HDU. hdu_bands : `~astropy.io.fits.BinTableHDU` The BANDS table HDU. format : {'gadf', 'fgst-ccube','fgst-template'} FITS format convention. Returns ------- map : `WcsNDMap` Wcs map """geom=WcsGeom.from_header(hdu.header,hdu_bands,format=format)shape=geom.axes.shapeshape_wcs=tuple([np.max(geom.npix[0]),np.max(geom.npix[1])])meta=cls._get_meta_from_header(hdu.header)unit=unit_from_fits_image_hdu(hdu.header)# TODO: Should we support extracting slices?ifisinstance(hdu,fits.BinTableHDU):map_out=cls(geom,meta=meta,unit=unit)pix=hdu.data.field("PIX")pix=np.unravel_index(pix,shape_wcs[::-1])vals=hdu.data.field("VALUE")if"CHANNEL"inhdu.data.columns.namesandshape:chan=hdu.data.field("CHANNEL")chan=np.unravel_index(chan,shape[::-1])idx=chan+pixelse:idx=pixmap_out.set_by_idx(idx[::-1],vals)else:ifany(xinhdu.name.lower()forxin["mask","is_ul","success"]):data=hdu.data.astype(bool)else:data=hdu.datamap_out=cls(geom=geom,meta=meta,data=data,unit=unit)returnmap_out
[docs]definterp_by_coord(self,coords,method="linear",fill_value=None,values_scale="lin"):"""Interpolate map values at the given map coordinates. Parameters ---------- coords : tuple, dict or `~gammapy.maps.MapCoord` Coordinate arrays for each dimension of the map. Tuple should be ordered as (lon, lat, x_0, ..., x_n) where x_i are coordinates for non-spatial dimensions of the map. "lon" and "lat" are optional and will be taken at the center of the region by default. method : {"linear", "nearest"} Method to interpolate data values. By default linear interpolation is performed. fill_value : None or float value The value to use for points outside of the interpolation domain. If None, values outside the domain are extrapolated. values_scale : {"lin", "log", "sqrt"} Optional value scaling. Default is "lin". Returns ------- vals : `~numpy.ndarray` Interpolated pixel values. """ifself.geom.is_regular:pix=self.geom.coord_to_pix(coords)returnself.interp_by_pix(pix,method=method,fill_value=fill_value,values_scale=values_scale)else:returnself._interp_by_coord_griddata(coords,method=method)
[docs]definterp_by_pix(self,pix,method="linear",fill_value=None,values_scale="lin"):ifnotself.geom.is_regular:raiseValueError("interp_by_pix only supported for regular geom.")grid_pix=[np.arange(n,dtype=float)forninself.data.shape[::-1]]ifnp.any(np.isfinite(self.data)):data=self.data.copy().Tdata[~np.isfinite(data)]=0.0else:data=self.data.Tfn=ScaledRegularGridInterpolator(grid_pix,data,fill_value=None,bounds_error=False,method=method,values_scale=values_scale,)interp_data=fn(tuple(pix),clip=False)iffill_valueisnotNone:idxs=self.geom.pix_to_idx(pix,clip=False)invalid=np.broadcast_arrays(*[idx==-1foridxinidxs])mask=np.any(invalid,axis=0)ifnotinterp_data.shape:mask=mask.squeeze()interp_data[mask]=fill_valueinterp_data[~np.isfinite(interp_data)]=fill_valuereturninterp_data
def_pad_spatial(self,pad_width,axis_name=None,mode="constant",cval=0,method="linear"):ifaxis_nameisNone:ifnp.isscalar(pad_width):pad_width=(pad_width,pad_width)iflen(pad_width)==2:pad_width+=(0,)*(self.geom.ndim-2)geom=self.geom._pad_spatial(pad_width[:2])ifself.geom.is_regularandmode!="interp":returnself._pad_np(geom,pad_width,mode,cval)else:returnself._pad_coadd(geom,pad_width,mode,cval,method)def_pad_np(self,geom,pad_width,mode,cval):"""Pad a map using ``numpy.pad``. This method only works for regular geometries but should be more efficient when working with large maps. """kwargs={}ifmode=="constant":kwargs["constant_values"]=cvalpad_width=[(t,t)fortinpad_width]data=np.pad(self.data,pad_width[::-1],mode,**kwargs)returnself._init_copy(geom=geom,data=data)def_pad_coadd(self,geom,pad_width,mode,cval,method):"""Pad a map manually by coadding the original map with the new map."""idx_in=self.geom.get_idx(flat=True)idx_in=tuple([t+wfort,winzip(idx_in,pad_width)])[::-1]idx_out=geom.get_idx(flat=True)[::-1]map_out=self._init_copy(geom=geom,data=None)map_out.coadd(self)ifmode=="constant":pad_msk=np.zeros_like(map_out.data,dtype=bool)pad_msk[idx_out]=Truepad_msk[idx_in]=Falsemap_out.data[pad_msk]=cvalelifmode=="interp":coords=geom.pix_to_coord(idx_out[::-1])m=self.geom.contains(coords)coords=tuple([c[~m]forcincoords])vals=self.interp_by_coord(coords,method=method)map_out.set_by_coord(coords,vals)else:raiseValueError(f"Invalid mode: {mode!r}")returnmap_out
[docs]defcrop(self,crop_width):ifnp.isscalar(crop_width):crop_width=(crop_width,crop_width)geom=self.geom.crop(crop_width)ifself.geom.is_regular:slices=[slice(None)]*len(self.geom.axes)slices+=[slice(crop_width[1],int(self.geom.npix[1]-crop_width[1])),slice(crop_width[0],int(self.geom.npix[0]-crop_width[0])),]data=self.data[tuple(slices)]map_out=self._init_copy(geom=geom,data=data)else:# FIXME: This could be done more efficiently by# constructing the appropriate slices for each image planemap_out=self._init_copy(geom=geom,data=None)map_out.coadd(self)returnmap_out
[docs]defplot(self,ax=None,fig=None,add_cbar=False,stretch="linear",**kwargs):""" Plot image on matplotlib WCS axes. Parameters ---------- ax : `~astropy.visualization.wcsaxes.WCSAxes`, optional WCS axis object to plot on. fig : `~matplotlib.figure.Figure` Figure object. add_cbar : bool Add color bar? stretch : str Passed to `astropy.visualization.simple_norm`. **kwargs : dict Keyword arguments passed to `~matplotlib.pyplot.imshow`. Returns ------- ax : `~astropy.visualization.wcsaxes.WCSAxes` WCS axis object """fromastropy.visualizationimportsimple_normifnotself.geom.is_flat:raiseTypeError("Use .plot_interactive() for Map dimension > 2")ax=self._plot_default_axes(ax=ax)iffigisNone:fig=plt.gcf()ifself.geom.is_image:data=self.data.astype(float)else:axis=tuple(np.arange(len(self.geom.axes)))data=np.squeeze(self.data,axis=axis).astype(float)kwargs.setdefault("interpolation","nearest")kwargs.setdefault("origin","lower")kwargs.setdefault("cmap","afmhot")mask=np.isfinite(data)ifmask.any():min_cut,max_cut=kwargs.pop("vmin",None),kwargs.pop("vmax",None)norm=simple_norm(data[mask],stretch,min_cut=min_cut,max_cut=max_cut)kwargs.setdefault("norm",norm)im=ax.imshow(data,**kwargs)ifadd_cbar:fig.colorbar(im,ax=ax,label=str(self.unit))ifself.geom.is_allsky:ax=self._plot_format_allsky(ax)else:ax=self._plot_format(ax)# without this the axis limits are changed when calling scatterax.autoscale(enable=False)returnax
[docs]defplot_mask(self,ax=None,**kwargs):"""Plot the mask as a shaded area Parameters ---------- ax : `~astropy.visualization.wcsaxes.WCSAxes`, optional WCS axis object to plot on. **kwargs : dict Keyword arguments passed to `~matplotlib.pyplot.contourf` Returns ------- ax : `~astropy.visualization.wcsaxes.WCSAxes`, optional WCS axis object to plot on. """ifnotself.geom.is_flat:raiseTypeError("Use .plot_interactive() for Map dimension > 2")ifnotself.is_mask:raiseValueError("`.plot_mask()` only supports maps containing boolean values.")ax=self._plot_default_axes(ax=ax)kwargs.setdefault("alpha",0.5)kwargs.setdefault("colors","w")data=np.squeeze(self.data).astype(float)ax.contourf(data,levels=[0,0.5],**kwargs)ifself.geom.is_allsky:ax=self._plot_format_allsky(ax)else:ax=self._plot_format(ax)# without this the axis limits are changed when calling scatterax.autoscale(enable=False)returnax
def_plot_default_axes(self,ax):fromastropy.visualization.wcsaxes.frameimportEllipticalFrameifaxisNone:fig=plt.gcf()ifself.geom.projectionin["AIT"]:ax=fig.add_subplot(1,1,1,projection=self.geom.wcs,frame_class=EllipticalFrame)else:ax=fig.add_subplot(1,1,1,projection=self.geom.wcs)returnax@staticmethoddef_plot_format(ax):try:ax.coords["glon"].set_axislabel("Galactic Longitude")ax.coords["glat"].set_axislabel("Galactic Latitude")exceptKeyError:ax.coords["ra"].set_axislabel("Right Ascension")ax.coords["dec"].set_axislabel("Declination")exceptAttributeError:log.info("Can't set coordinate axes. No WCS information available.")returnaxdef_plot_format_allsky(self,ax):# Remove frameax.coords.frame.set_linewidth(0)# Set plot axis limitsxmin,_=self.geom.to_image().coord_to_pix({"lon":180,"lat":0})xmax,_=self.geom.to_image().coord_to_pix({"lon":-180,"lat":0})_,ymin=self.geom.to_image().coord_to_pix({"lon":0,"lat":-90})_,ymax=self.geom.to_image().coord_to_pix({"lon":0,"lat":90})ax.set_xlim(xmin,xmax)ax.set_ylim(ymin,ymax)ax.text(0,ymax,self.geom.frame+" coords")# Grid and ticksglon_spacing,glat_spacing=45,15lon,lat=ax.coordslon.set_ticks(spacing=glon_spacing*u.deg,color="w",alpha=0.8)lat.set_ticks(spacing=glat_spacing*u.deg)lon.set_ticks_visible(False)lon.set_major_formatter("d")lat.set_major_formatter("d")lon.set_ticklabel(color="w",alpha=0.8)lon.grid(alpha=0.2,linestyle="solid",color="w")lat.grid(alpha=0.2,linestyle="solid",color="w")returnax
[docs]defcutout_and_mask_region(self,region=None):"""Compute cutout and mask for a given region of the map. The function will estimate the minimal size of the cutout, which encloses the region. Parameters ---------- region: `~regions.Region` Extended region Returns ------- cutout, mask : tuple of `WcsNDMap` Cutout and mask map """fromgammapy.mapsimportRegionGeomifregionisNone:region=self.geom.footprint_rectangle_sky_regiongeom=RegionGeom.from_regions(regions=region,wcs=self.geom.wcs)cutout=self.cutout(position=geom.center_skydir,width=geom.width)mask=cutout.geom.to_image().region_mask([region])returnself.__class__(data=cutout.data,geom=cutout.geom,unit=self.unit),mask
[docs]defto_region_nd_map(self,region=None,func=np.nansum,weights=None,method="nearest"):"""Get region ND map in a given region. By default the whole map region is considered. Parameters ---------- region: `~regions.Region` or `~astropy.coordinates.SkyCoord` Region. func : numpy.func Function to reduce the data. Default is np.nansum. For boolean Map, use np.any or np.all. weights : `WcsNDMap` Array to be used as weights. The geometry must be equivalent. method : {"nearest", "linear"} How to interpolate if a position is given. Returns ------- spectrum : `~gammapy.maps.RegionNDMap` Spectrum in the given region. """fromgammapy.mapsimportRegionGeom,RegionNDMapifregionisNone:region=self.geom.footprint_rectangle_sky_regionifweightsisnotNone:ifnotself.geom==weights.geom:raiseValueError("Incompatible spatial geoms between map and weights")geom=RegionGeom.from_regions(regions=region,axes=self.geom.axes,wcs=self.geom.wcs)ifgeom.is_all_point_sky_regions:coords=geom.get_coord()data=self.interp_by_coord(coords=coords,method=method)ifweightsisnotNone:data*=weights.interp_by_coord(coords=coords,method=method)# Casting needed as interp_by_coord transforms booleandata=data.astype(self.data.dtype)else:cutout,mask=self.cutout_and_mask_region(region=region)ifweightsisnotNone:weights_cutout=weights.cutout(position=geom.center_skydir,width=geom.width)cutout.data*=weights_cutout.dataidx_y,idx_x=np.where(mask)data=func(cutout.data[...,idx_y,idx_x],axis=-1)returnRegionNDMap(geom=geom,data=data,unit=self.unit,meta=self.meta.copy())
[docs]defto_region_nd_map_histogram(self,region=None,bins_axis=None,nbin=100,density=False):"""Convert map into region map by histogramming. By default it creates a linearly spaced axis with 100 bins between (-max(abs(data)), max(abs(data))) within the given region. Parameters ---------- region: `~regions.Region` Region to histogram over. bins_axis : `MapAxis` Binning of the histogram. nbin : int Number of bins to use if no bins_axis is given. density : bool Normalize integral of the histogram to 1. Examples -------- This is how to use the method to create energy dependent histograms: :: from gammapy.maps import MapAxis, Map import numpy as np random_state = np.random.RandomState(seed=0) energy_axis = MapAxis.from_energy_bounds("1 TeV", "10 TeV", nbin=3) data = Map.create(axes=[energy_axis], width=10, unit="cm2 s-1", binsz=0.02) data.data = random_state.normal( size=data.data.shape, loc=0, scale=np.array([1.0, 2.0, 3.0]).reshape((-1, 1, 1)) ) hist = data.to_region_nd_map_histogram() hist.plot(axis_name="bins") Returns ------- region_map : `RegionNDMap` Region map with histogram. """fromgammapy.mapsimportMapAxis,RegionGeom,RegionNDMapifisinstance(region,(PointSkyRegion,SkyCoord)):raiseValueError("Histogram method not supported for point regions")cutout,mask=self.cutout_and_mask_region(region=region)idx_y,idx_x=np.where(mask)quantity=cutout.quantity[...,idx_y,idx_x]value=np.abs(quantity).max()ifbins_axisisNone:bins_axis=MapAxis.from_bounds(-value,value,nbin=nbin,interp="lin",unit=self.unit,name="bins",)ifnotbins_axis.unit.is_equivalent(self.unit):raiseValueError("Unit of bins_axis must be equivalent to unit of map.")axes=[bins_axis]+list(self.geom.axes)geom_hist=RegionGeom(region=region,axes=axes,wcs=self.geom.wcs)# This is likely not the most efficient way to do thisdata=np.apply_along_axis(lambdaa:np.histogram(a,bins=bins_axis.edges,density=density)[0],axis=-1,arr=quantity,)ifdensity:unit=1.0/bins_axis.unitdata=data.valueelse:unit=""returnRegionNDMap.from_geom(geom=geom_hist,data=data,unit=unit)
[docs]defmask_contains_region(self,region):"""Check if input region is contained in a boolean mask map. Parameters ---------- region: `~regions.SkyRegion` or `~regions.PixRegion` Region or list of Regions (pixel or sky regions accepted). Returns ------- contained : bool Whether region is contained in the mask """ifnotself.is_mask:raiseValueError("mask_contains_region is only supported for boolean masks")ifnotself.geom.is_image:raiseValueError("Method only supported for 2D images")ifisinstance(region,SkyRegion):region=region.to_pixel(self.geom.wcs)ifisinstance(region,PointPixelRegion):lon,lat=region.center.x,region.center.ycontains=self.get_by_pix((lon,lat))else:idx=self.geom.get_idx()coords_pix=PixCoord(idx[0][self.data],idx[1][self.data])contains=region.contains(coords_pix)returnnp.any(contains)
[docs]defbinary_erode(self,width,kernel="disk",use_fft=True):"""Binary erosion of boolean mask removing a given margin Parameters ---------- width : `~astropy.units.Quantity`, str or float If a float is given it interpreted as width in pixels. If an (angular) quantity is given it converted to pixels using ``geom.wcs.wcs.cdelt``. The width corresponds to radius in case of a disk kernel, and the side length in case of a box kernel. kernel : {'disk', 'box'} Kernel shape use_fft : bool Use `scipy.signal.fftconvolve` if True (default) and `scipy.ndimage.binary_erosion` otherwise. Returns ------- map : `WcsNDMap` Eroded mask map """ifnotself.is_mask:raiseValueError("Binary operations only supported for boolean masks")structure=self.geom.binary_structure(width=width,kernel=kernel)ifuse_fft:returnself.convolve(structure.squeeze(),method="fft")>(structure.sum()-1)data=ndi.binary_erosion(self.data,structure=structure)returnself._init_copy(data=data)
[docs]defbinary_dilate(self,width,kernel="disk",use_fft=True):"""Binary dilation of boolean mask adding a given margin Parameters ---------- width : tuple of `~astropy.units.Quantity` Angular sizes of the margin in (lon, lat) in that specific order. If only one value is passed, the same margin is applied in (lon, lat). kernel : {'disk', 'box'} Kernel shape use_fft : bool Use `scipy.signal.fftconvolve` if True (default) and `scipy.ndimage.binary_dilation` otherwise. Returns ------- map : `WcsNDMap` Dilated mask map """ifnotself.is_mask:raiseValueError("Binary operations only supported for boolean masks")structure=self.geom.binary_structure(width=width,kernel=kernel)ifuse_fft:returnself.convolve(structure.squeeze(),method="fft")>1data=ndi.binary_dilation(self.data,structure=structure)returnself._init_copy(data=data)
[docs]defconvolve(self,kernel,method="fft",mode="same"):"""Convolve map with a kernel. If the kernel is two dimensional, it is applied to all image planes likewise. If the kernel is higher dimensional, it should either match the map in the number of dimensions or the map must be an image (no non-spatial axes). In that case, the corresponding kernel is selected and applied to every image plane or to the single input image respectively. Parameters ---------- kernel : `~gammapy.irf.PSFKernel` or `numpy.ndarray` Convolution kernel. method : str The method used by `~scipy.signal.convolve`. Default is 'fft'. mode : str The convolution mode used by `~scipy.signal.convolve`. Default is 'same'. Returns ------- map : `WcsNDMap` Convolved map. """fromgammapy.irfimportPSFKernelconvolve=scipy.signal.convolveifself.geom.is_imageandnotisinstance(kernel,PSFKernel):ifkernel.ndim>2:raiseValueError("Image convolution with 3D kernel requires a PSFKernel object")geom=self.geom.copy()ifisinstance(kernel,PSFKernel):kmap=kernel.psf_kernel_mapifnotnp.allclose(self.geom.pixel_scales.deg,kmap.geom.pixel_scales.deg,rtol=1e-5):raiseValueError("Pixel size of kernel and map not compatible.")kernel=kmap.data.astype(np.float32)ifself.geom.is_image:geom=geom.to_cube(kmap.geom.axes)ifmode=="full":pad_width=[0.5*(width-1)forwidthinkernel.shape[-2:]]geom=geom.pad(pad_width,axis_name=None)elifmode=="valid":raiseNotImplementedError("WcsNDMap.convolve: mode='valid' is not supported.")data=np.empty(geom.data_shape,dtype=np.float32)shape_axes_kernel=kernel.shape[slice(0,-2)]iflen(shape_axes_kernel)>0:ifnotgeom.shape_axes==shape_axes_kernel:raiseValueError(f"Incompatible shape between data {geom.shape_axes}"" and kernel {shape_axes_kernel}")ifself.geom.is_imageandkernel.ndim==3:foridxinrange(kernel.shape[0]):data[idx]=convolve(self.data.astype(np.float32),kernel[idx],method=method,mode=mode)else:forimg,idxinself.iter_by_image_data():ikern=Ellipsisifkernel.ndim==2elseidxdata[idx]=convolve(img.astype(np.float32),kernel[ikern],method=method,mode=mode)returnself._init_copy(data=data,geom=geom)
[docs]defsmooth(self,width,kernel="gauss",**kwargs):"""Smooth the map. Iterates over 2D image planes, processing one at a time. Parameters ---------- width : `~astropy.units.Quantity`, str or float Smoothing width given as quantity or float. If a float is given it interpreted as smoothing width in pixels. If an (angular) quantity is given it converted to pixels using ``geom.wcs.wcs.cdelt``. It corresponds to the standard deviation in case of a Gaussian kernel, the radius in case of a disk kernel, and the side length in case of a box kernel. kernel : {'gauss', 'disk', 'box'} Kernel shape kwargs : dict Keyword arguments passed to `~ndi.uniform_filter` ('box'), `~ndi.gaussian_filter` ('gauss') or `~ndi.convolve` ('disk'). Returns ------- image : `WcsNDMap` Smoothed image (a copy, the original object is unchanged). """ifisinstance(width,(u.Quantity,str)):width=u.Quantity(width)/self.geom.pixel_scales.mean()width=width.to_value("")smoothed_data=np.empty(self.data.shape,dtype=float)forimg,idxinself.iter_by_image_data():img=img.astype(float)ifkernel=="gauss":data=ndi.gaussian_filter(img,width,**kwargs)elifkernel=="disk":disk=Tophat2DKernel(width)disk.normalize("integral")data=ndi.convolve(img,disk.array,**kwargs)elifkernel=="box":data=ndi.uniform_filter(img,width,**kwargs)else:raiseValueError(f"Invalid kernel: {kernel!r}")smoothed_data[idx]=datareturnself._init_copy(data=smoothed_data)
[docs]defcutout(self,position,width,mode="trim",odd_npix=False):""" Create a cutout around a given position. Parameters ---------- position : `~astropy.coordinates.SkyCoord` Center position of the cutout region. width : tuple of `~astropy.coordinates.Angle` Angular sizes of the region in (lon, lat) in that specific order. If only one value is passed, a square region is extracted. mode : {'trim', 'partial', 'strict'} Mode option for Cutout2D, for details see `~astropy.nddata.utils.Cutout2D`. odd_npix : bool Force width to odd number of pixels. Returns ------- cutout : `~gammapy.maps.WcsNDMap` Cutout map """geom_cutout=self.geom.cutout(position=position,width=width,mode=mode,odd_npix=odd_npix)cutout_info=geom_cutout.cutout_slices(self.geom,mode=mode)slices=cutout_info["parent-slices"]parent_slices=Ellipsis,slices[0],slices[1]slices=cutout_info["cutout-slices"]cutout_slices=Ellipsis,slices[0],slices[1]data=np.zeros(shape=geom_cutout.data_shape,dtype=self.data.dtype)data[cutout_slices]=self.data[parent_slices]returnself._init_copy(geom=geom_cutout,data=data)
[docs]defstack(self,other,weights=None,nan_to_num=True):"""Stack cutout into map. Parameters ---------- other : `WcsNDMap` Other map to stack weights : `WcsNDMap` Array to be used as weights. The spatial geometry must be equivalent to `other` and additional axes must be broadcastable. nan_to_num: bool Non-finite values are replaced by zero if True (default). """ifself.geom==other.geom:parent_slices,cutout_slices=None,Noneelifself.geom.is_aligned(other.geom):cutout_slices=other.geom.cutout_slices(self.geom)slices=cutout_slices["parent-slices"]parent_slices=Ellipsis,slices[0],slices[1]slices=cutout_slices["cutout-slices"]cutout_slices=Ellipsis,slices[0],slices[1]else:raiseValueError("Can only stack equivalent maps or cutout of the same map.")data=other.quantity[cutout_slices].to_value(self.unit)ifnan_to_num:not_finite=~np.isfinite(data)ifnp.any(not_finite):data=data.copy()data[not_finite]=0ifweightsisnotNone:ifnotother.geom.to_image()==weights.geom.to_image():raiseValueError("Incompatible spatial geoms between map and weights")data=data*weights.data[cutout_slices]self.data[parent_slices]+=data