# Licensed under a 3-clause BSD style license - see LICENSE.rst"""Utilities for dealing with HEALPix projections and mappings."""importcopyimportnumpyasnpfromastropyimportunitsasufromastropy.coordinatesimportSkyCoordfromastropy.ioimportfitsfromastropy.unitsimportQuantityfromgammapy.utils.arrayimportis_power2from..axesimportMapAxesfrom..coordimportMapCoord,skycoord_to_lonlatfrom..geomimportGeom,pix_tuple_to_idxfrom..utilsimportINVALID_INDEX,coordsys_to_frame,frame_to_coordsysfrom.ioimportHPX_FITS_CONVENTIONS,HpxConvfrom.utilsimport(coords_to_vec,get_nside_from_pix_size,get_pix_size_from_nside,get_subpixels,get_superpixels,match_hpx_pix,nside_to_order,parse_hpxregion,ravel_hpx_index,unravel_hpx_index,)# Not sure if we should expose this in the docs or not:# HPX_FITS_CONVENTIONS, HpxConv__all__=["HpxGeom"]
[docs]classHpxGeom(Geom):"""Geometry class for HEALPix maps. This class performs mapping between partial-sky indices (pixel number within a HEALPix region) and all-sky indices (pixel number within an all-sky HEALPix map). Multi-band HEALPix geometries use a global indexing scheme that assigns a unique pixel number based on the all-sky index and band index. In the single-band case the global index is the same as the HEALPix index. By default, the constructor will return an all-sky map. Partial-sky maps can be defined with the ``region`` argument. Parameters ---------- nside : `~numpy.ndarray` HEALPix NSIDE parameter, the total number of pixels is 12*nside*nside. For multi-dimensional maps one can pass either a single ``nside`` value or a vector of ``nside`` values defining the pixel size for each image plane. If ``nside`` is not a scalar then its dimensionality should match that of the non-spatial axes. If nest is True, ``nside`` must be a power of 2, less than 2**30. nest : bool Indexing scheme. If True, "NESTED" scheme. If False, "RING" scheme. frame : {"icrs", "galactic"} Coordinate system. Default is "icrs". region : str or tuple Spatial geometry for partial-sky maps. If None, the map will encompass the whole sky. String input will be parsed according to HPX_REG header keyword conventions. Tuple input can be used to define an explicit list of pixels encompassed by the geometry. axes : list Axes for non-spatial dimensions. """is_hpx=Trueis_region=Falsedef__init__(self,nside,nest=True,frame="icrs",region=None,axes=None):fromhealpy.pixelfuncimportcheck_nsidecheck_nside(nside,nest=nest)self._nside=np.array(nside,ndmin=1)self._axes=MapAxes.from_default(axes,n_spatial_axes=1)ifself.nside.size>1andself.nside.shape!=self.shape_axes:raiseValueError("Wrong dimensionality for nside. nside must ""be a scalar or have a dimensionality consistent ""with the axes argument.")self._frame=frameself._nest=nestself._ipix=Noneself._region=regionself._create_lookup(region)self._npix=self._npix*np.ones(self.shape_axes,dtype=int)def_create_lookup(self,region):"""Create local-to-global pixel lookup table."""ifisinstance(region,str):ipix=[self.get_index_list(nside,self._nest,region)fornsideinself._nside.flat]self._ipix=[ravel_hpx_index((p,i*np.ones_like(p)),np.ravel(self.npix_max))fori,pinenumerate(ipix)]self._region=regionself._indxschm="EXPLICIT"self._npix=np.array([len(t)fortinself._ipix])ifself.nside.ndim>1:self._npix=self._npix.reshape(self.nside.shape)self._ipix=np.concatenate(self._ipix)elifisinstance(region,tuple):region=[np.asarray(t)fortinregion]m=np.any(np.stack([t>=0fortinregion]),axis=0)region=[t[m]fortinregion]self._ipix=ravel_hpx_index(region,self.npix_max)self._ipix=np.unique(self._ipix)region=unravel_hpx_index(self._ipix,self.npix_max)self._region="explicit"self._indxschm="EXPLICIT"iflen(region)==1:self._npix=np.array([len(region[0])])else:self._npix=np.zeros(self.shape_axes,dtype=int)idx=np.ravel_multi_index(region[1:],self.shape_axes)cnt=np.unique(idx,return_counts=True)self._npix.flat[cnt[0]]=cnt[1]elifregionisNone:self._region=Noneself._indxschm="IMPLICIT"self._npix=self.npix_maxelse:raiseValueError(f"Invalid region string: {region!r}")
[docs]deflocal_to_global(self,idx_local):"""Compute a global index (all-sky) from a local (partial-sky) index. Parameters ---------- idx_local : tuple A tuple of pixel indices with local HEALPix pixel indices. Returns ------- idx_global : tuple A tuple of pixel index vectors with global HEALPix pixel indices. """ifself._ipixisNone:returnidx_localifself.nside.size>1:idx=ravel_hpx_index(idx_local,self._npix)else:idx_tmp=tuple([idx_local[0]]+[np.zeros(t.shape,dtype=int)fortinidx_local[1:]])idx=ravel_hpx_index(idx_tmp,self._npix)idx_global=unravel_hpx_index(self._ipix[idx],self.npix_max)returnidx_global[:1]+tuple(idx_local[1:])
[docs]defglobal_to_local(self,idx_global,ravel=False):"""Compute global (all-sky) index from a local (partial-sky) index. Parameters ---------- idx_global : tuple A tuple of pixel indices with global HEALPix pixel indices. ravel : bool, optional Return a raveled index. Default is False. Returns ------- idx_local : tuple A tuple of pixel indices with local HEALPix pixel indices. """if(isinstance(idx_global,int)or(isinstance(idx_global,tuple)andisinstance(idx_global[0],int))orisinstance(idx_global,np.ndarray)):idx_global=unravel_hpx_index(np.array(idx_global,ndmin=1),self.npix_max)ifself.nside.size==1:idx=np.array(idx_global[0],ndmin=1)else:idx=ravel_hpx_index(idx_global,self.npix_max)ifself._ipixisnotNone:retval=np.full(idx.size,-1,"i")m=np.isin(idx.flat,self._ipix)retval[m]=np.searchsorted(self._ipix,idx.flat[m])retval=retval.reshape(idx.shape)else:retval=idxifself.nside.size==1:idx_local=tuple([retval]+list(idx_global[1:]))else:idx_local=unravel_hpx_index(retval,self._npix)m=np.any(np.stack([t==INVALID_INDEX.intfortinidx_local]),axis=0)fori,tinenumerate(idx_local):idx_local[i][m]=INVALID_INDEX.intifnotravel:returnidx_localelse:returnravel_hpx_index(idx_local,self.npix)
[docs]defcutout(self,position,width,**kwargs):"""Create a cutout around a given position. Parameters ---------- position : `~astropy.coordinates.SkyCoord` Center position of the cutout region. width : `~astropy.coordinates.Angle` or `~astropy.units.Quantity` Diameter of the circular cutout region. Returns ------- cutout : `~gammapy.maps.WcsNDMap` Cutout map. """ifnotself.is_regular:raiseValueError("Can only do a cutout from a regular map.")width=u.Quantity(width,"deg").valuereturnself.create(nside=self.nside,nest=self.nest,width=width,skydir=position,frame=self.frame,axes=self.axes,)
[docs]defcoord_to_pix(self,coords):importhealpyashpcoords=MapCoord.create(coords,frame=self.frame,axis_names=self.axes.names).broadcastedtheta,phi=coords.theta,coords.phiifself.axes:idxs=self.axes.coord_to_idx(coords,clip=True)bins=self.axes.coord_to_pix(coords)# FIXME: Figure out how to handle coordinates out of# bounds of non-spatial dimensionsifself.nside.size>1:nside=self.nside[tuple(idxs)]else:nside=self.nsidem=~np.isfinite(theta)theta[m]=0.0phi[m]=0.0pix=hp.ang2pix(nside,theta,phi,nest=self.nest)pix=tuple([pix])+binsifnp.any(m):forpinpix:p[m]=INVALID_INDEX.intelse:pix=(hp.ang2pix(self.nside,theta,phi,nest=self.nest),)returnpix
[docs]defpix_to_idx(self,pix,clip=False):# FIXME: Look for better method to clip HPX indices# TODO: copy idx to avoid modifying input pix?# pix_tuple_to_idx seems to always make a copy!?idx=pix_tuple_to_idx(pix)idx_local=self.global_to_local(idx)fori,_inenumerate(idx):ifclip:ifi>0:np.clip(idx[i],0,self.axes[i-1].nbin-1,out=idx[i])else:np.clip(idx[i],0,None,out=idx[i])else:ifi>0:mask=(idx[i]<0)|(idx[i]>=self.axes[i-1].nbin)np.putmask(idx[i],mask,-1)else:mask=(idx_local[i]<0)|(idx[i]<0)np.putmask(idx[i],mask,-1)returntuple(idx)
@propertydefaxes(self):"""List of non-spatial axes."""returnself._axes@propertydefaxes_names(self):"""All axes names."""return["skycoord"]+self.axes.names@propertydefshape_axes(self):"""Shape of non-spatial axes."""returnself.axes.shape@propertydefdata_shape(self):"""Shape of the `~numpy.ndarray` matching this geometry."""npix_shape=tuple([np.max(self.npix)])return(npix_shape+self.axes.shape)[::-1]@propertydefdata_shape_axes(self):"""Shape of data of the non-spatial axes and unit spatial axes."""returnself.axes.shape[::-1]+(1,)@propertydefndim(self):"""Number of dimensions as an integer."""returnlen(self._axes)+2@propertydefordering(self):"""HEALPix ordering ('NESTED' or 'RING')."""return"NESTED"ifself.nestelse"RING"@propertydefnside(self):"""NSIDE in each band."""returnself._nside@propertydeforder(self):"""The order in each band (``NSIDE = 2 ** ORDER``). Set to -1 for bands with NSIDE that is not a power of 2. """returnnside_to_order(self.nside)@propertydefnest(self):"""Whether HEALPix order is nested as a boolean."""returnself._nest@propertydefnpix(self):"""Number of pixels in each band. For partial-sky geometries this can be less than the number of pixels for the band NSIDE. """returnself._npix@propertydefnpix_max(self):"""Maximum number of pixels."""maxpix=12*self.nside**2returnmaxpix*np.ones(self.shape_axes,dtype=int)@propertydefframe(self):returnself._frame@propertydefprojection(self):"""Map projection."""return"HPX"@propertydefregion(self):"""Region string."""returnself._region@propertydefis_allsky(self):"""Flag for all-sky maps."""returnself._regionisNone@propertydefis_regular(self):"""Flag identifying whether this geometry is regular in non-spatial dimensions. False for multi-resolution or irregular geometries. If True, all image planes have the same pixel geometry. """ifself.nside.size>1orself.region=="explicit":returnFalseelse:returnTrue@propertydefcenter_coord(self):"""Map coordinates of the center of the geometry as a tuple."""lon,lat,frame=skycoord_to_lonlat(self.center_skydir)returntuple([lon,lat])+self.axes.center_coord@propertydefcenter_pix(self):"""Pixel coordinates of the center of the geometry as a tuple."""returnself.coord_to_pix(self.center_coord)@propertydefcenter_skydir(self):"""Sky coordinate of the center of the geometry. Returns ------- center : `~astropy.coordinates.SkyCoord` Center position. """# TODO: simplifyimporthealpyashpifself.is_allsky:lon,lat=0.0,0.0elifself.region=="explicit":idx=unravel_hpx_index(self._ipix,self.npix_max)nside=self._get_nside(idx)vec=hp.pix2vec(nside,idx[0],nest=self.nest)vec=np.array([np.mean(t)fortinvec])lonlat=hp.vec2ang(vec,lonlat=True)lon,lat=lonlat[0],lonlat[1]else:tokens=parse_hpxregion(self.region)iftokens[0]in["DISK","DISK_INC"]:lon,lat=float(tokens[1]),float(tokens[2])eliftokens[0]=="HPX_PIXEL":nside_pix=int(tokens[2])ipix_pix=int(tokens[3])iftokens[1]=="NESTED":nest_pix=Trueeliftokens[1]=="RING":nest_pix=Falseelse:raiseValueError(f"Invalid ordering scheme: {tokens[1]!r}")theta,phi=hp.pix2ang(nside_pix,ipix_pix,nest_pix)lat=np.degrees((np.pi/2)-theta)lon=np.degrees(phi)returnSkyCoord(lon,lat,frame=self.frame,unit="deg")@propertydefpixel_scales(self):self.angle_="""Pixel scale. Returns ------- angle: `~astropy.coordinates.Angle` """returnget_pix_size_from_nside(self.nside)*u.deg
[docs]definterp_weights(self,coords,idxs=None):"""Get interpolation weights for given coordinates. Parameters ---------- coords : `MapCoord` or dict Input coordinates. idxs : `~numpy.ndarray`, optional Indices for non-spatial axes. Default is None. Returns ------- weights : `~numpy.ndarray` Interpolation weights. """importhealpyashpcoords=MapCoord.create(coords,frame=self.frame).broadcastedifidxsisNone:idxs=self.coord_to_idx(coords,clip=True)[1:]theta,phi=coords.theta,coords.phim=~np.isfinite(theta)theta[m]=0phi[m]=0ifnotself.is_regular:nside=self.nside[tuple(idxs)]else:nside=self.nsidepix,wts=hp.get_interp_weights(nside,theta,phi,nest=self.nest)wts[:,m]=0pix[:,m]=INVALID_INDEX.intifnotself.is_regular:pix_local=[self.global_to_local([pix]+list(idxs))[0]]else:pix_local=[self.global_to_local(pix,ravel=True)]# If a pixel lies outside of the geometry set its index to the center pixelm=pix_local[0]==INVALID_INDEX.intifm.any():coords_ctr=[coords.lon,coords.lat]coords_ctr+=[ax.pix_to_coord(t)forax,tinzip(self.axes,idxs)]idx_ctr=self.coord_to_idx(coords_ctr)idx_ctr=self.global_to_local(idx_ctr)pix_local[0][m]=(idx_ctr[0]*np.ones(pix.shape,dtype=int))[m]pix_local+=[np.broadcast_to(t,pix_local[0].shape)fortinidxs]returnpix_local,wts
@propertydefipix(self):"""HEALPix pixel and band indices for every pixel in the map."""returnself.get_idx()
[docs]defis_aligned(self,other):"""Check if HEALPix geoms and extra axes are aligned. Parameters ---------- other : `HpxGeom` Other geometry. Returns ------- aligned : bool Whether geometries are aligned. """foraxis,otheraxisinzip(self.axes,other.axes):ifaxis!=otheraxis:returnFalseifnotself.nside==other.nside:returnFalseelifnotself.frame==other.frame:returnFalseelifnotself.nest==other.nest:returnFalseelse:returnTrue
[docs]defto_nside(self,nside):"""Upgrade or downgrade the resolution to a given NSIDE. Parameters ---------- nside : int HEALPix NSIDE parameter. Returns ------- geom : `~HpxGeom` A HEALPix geometry object. """ifnotself.is_regular:raiseValueError("Upgrade and degrade only implemented for standard maps")axes=copy.deepcopy(self.axes)returnself.__class__(nside=nside,nest=self.nest,frame=self.frame,region=self.region,axes=axes)
[docs]defto_binsz(self,binsz):"""Change pixel size of the geometry. Parameters ---------- binsz : float or `~astropy.units.Quantity` New pixel size. A float is assumed to be in degree. Returns ------- geom : `WcsGeom` Geometry with new pixel size. """binsz=u.Quantity(binsz,"deg").valueifself.is_allsky:returnself.create(binsz=binsz,frame=self.frame,axes=copy.deepcopy(self.axes),)else:returnself.create(skydir=self.center_skydir,binsz=binsz,width=self.width.to_value("deg"),frame=self.frame,axes=copy.deepcopy(self.axes),)
[docs]defseparation(self,center):"""Compute sky separation with respect to a given center. Parameters ---------- center : `~astropy.coordinates.SkyCoord` Center position. Returns ------- separation : `~astropy.coordinates.Angle` Separation angle array (1D). """coord=self.to_image().get_coord()returncenter.separation(coord.skycoord)
[docs]defto_swapped(self):"""Geometry copy with swapped ORDERING (NEST->RING or vice versa). Returns ------- geom : `~HpxGeom` A HEALPix geometry object. """axes=copy.deepcopy(self.axes)returnself.__class__(self.nside,notself.nest,frame=self.frame,region=self.region,axes=axes,)
def_get_neighbors(self,idx):importhealpyashpnside=self._get_nside(idx)idx_nb=(hp.get_all_neighbours(nside,idx[0],nest=self.nest),)idx_nb+=tuple([t[None,...]*np.ones_like(idx_nb[0])fortinidx[1:]])returnidx_nbdef_pad_spatial(self,pad_width):ifself.is_allsky:raiseValueError("Cannot pad an all-sky map.")idx=self.get_idx(flat=True)idx_r=ravel_hpx_index(idx,self.npix_max)# TODO: Pre-filter indices to find those close to the edgeidx_nb=self._get_neighbors(idx)idx_nb=ravel_hpx_index(idx_nb,self.npix_max)for_inrange(pad_width):mask_edge=np.isin(idx_nb,idx_r,invert=True)idx_edge=idx_nb[mask_edge]idx_edge=np.unique(idx_edge)idx_r=np.sort(np.concatenate((idx_r,idx_edge)))idx_nb=unravel_hpx_index(idx_edge,self.npix_max)idx_nb=self._get_neighbors(idx_nb)idx_nb=ravel_hpx_index(idx_nb,self.npix_max)idx=unravel_hpx_index(idx_r,self.npix_max)returnself.__class__(self.nside.copy(),self.nest,frame=self.frame,region=idx,axes=copy.deepcopy(self.axes),)
[docs]defcrop(self,crop_width):ifself.is_allsky:raiseValueError("Cannot crop an all-sky map.")idx=self.get_idx(flat=True)idx_r=ravel_hpx_index(idx,self.npix_max)# TODO: Pre-filter indices to find those close to the edgeidx_nb=self._get_neighbors(idx)idx_nb=ravel_hpx_index(idx_nb,self.npix_max)for_inrange(crop_width):# Mask of pixels that have at least one neighbor not# contained in the geometrymask_edge=np.any(np.isin(idx_nb,idx_r,invert=True),axis=0)idx_r=idx_r[~mask_edge]idx_nb=idx_nb[:,~mask_edge]idx=unravel_hpx_index(idx_r,self.npix_max)returnself.__class__(self.nside.copy(),self.nest,frame=self.frame,region=idx,axes=copy.deepcopy(self.axes),)
[docs]defupsample(self,factor):ifnotis_power2(factor):raiseValueError("Upsample factor must be a power of 2.")ifself.is_allsky:returnself.__class__(self.nside*factor,self.nest,frame=self.frame,region=self.region,axes=copy.deepcopy(self.axes),)idx=list(self.get_idx(flat=True))nside=self._get_nside(idx)idx_new=get_subpixels(idx[0],nside,nside*factor,nest=self.nest)foriinrange(1,len(idx)):idx[i]=idx[i][...,None]*np.ones(idx_new.shape,dtype=int)idx[0]=idx_newreturnself.__class__(self.nside*factor,self.nest,frame=self.frame,region=tuple(idx),axes=copy.deepcopy(self.axes),)
[docs]defdownsample(self,factor,axis_name=None):ifnotis_power2(factor):raiseValueError("Downsample factor must be a power of 2.")ifaxis_nameisnotNone:raiseValueError("Currently the only valid axis name is None.")ifself.is_allsky:returnself.__class__(self.nside//factor,self.nest,frame=self.frame,region=self.region,axes=copy.deepcopy(self.axes),)idx=list(self.get_idx(flat=True))nside=self._get_nside(idx)idx_new=get_superpixels(idx[0],nside,nside//factor,nest=self.nest)idx[0]=idx_newreturnself.__class__(self.nside//factor,self.nest,frame=self.frame,region=tuple(idx),axes=copy.deepcopy(self.axes),)
[docs]@classmethoddefcreate(cls,nside=None,binsz=None,nest=True,frame="icrs",region=None,axes=None,skydir=None,width=None,):"""Create an HpxGeom object. Parameters ---------- nside : int or `~numpy.ndarray`, optional HEALPix NSIDE parameter. This parameter sets the size of the spatial pixels in the map. If nest is True, ``nside`` must be a power of 2, less than 2**30. Default is None. binsz : float or `~numpy.ndarray`, optional Approximate pixel size in degrees. An ``nside`` will be chosen that corresponds to a pixel size closest to this value. This option is superseded by ``nside``. Default is None. nest : bool, optional Indexing scheme. If True, "NESTED" scheme. If False, "RING" scheme. Default is True. frame : {"icrs", "galactic"} Coordinate system, either Galactic ("galactic") or Equatorial ("icrs"). Default is "icrs". region : str, optional HEALPix region string. Allows for partial-sky maps. Default is None. axes : list, optional List of axes for non-spatial dimensions. Default is None. skydir : tuple or `~astropy.coordinates.SkyCoord`, optional Sky position of map center. Can be either a SkyCoord object or a tuple of longitude and latitude in deg in the coordinate system of the map. Default is None. width : float, optional Diameter of the map in degrees. If set the map will encompass all pixels within a circular region centered on ``skydir``. Default is None. Returns ------- geom : `~HpxGeom` A HEALPix geometry object. Examples -------- >>> from gammapy.maps import HpxGeom, MapAxis >>> axis = MapAxis.from_bounds(0,1,2) >>> geom = HpxGeom.create(nside=16) # doctest: +SKIP >>> geom = HpxGeom.create(binsz=0.1, width=10.0) # doctest: +SKIP >>> geom = HpxGeom.create(nside=64, width=10.0, axes=[axis]) # doctest: +SKIP >>> geom = HpxGeom.create(nside=[32,64], width=10.0, axes=[axis]) # doctest: +SKIP """ifnsideisNoneandbinszisNone:raiseValueError("Either nside or binsz must be defined.")ifnsideisNoneandbinszisnotNone:nside=get_nside_from_pix_size(binsz)ifskydirisNone:lon,lat=(0.0,0.0)elifisinstance(skydir,tuple):lon,lat=skydirelifisinstance(skydir,SkyCoord):lon,lat,frame=skycoord_to_lonlat(skydir,frame=frame)else:raiseValueError(f"Invalid type for skydir: {type(skydir)!r}")ifregionisNoneandwidthisnotNone:region=f"DISK({lon},{lat},{width/2})"returncls(nside,nest=nest,frame=frame,region=region,axes=axes)
[docs]@classmethoddeffrom_header(cls,header,hdu_bands=None,format=None):"""Create an HPX object from a FITS header. Parameters ---------- header : `~astropy.io.fits.Header` The FITS header. hdu_bands : `~astropy.io.fits.BinTableHDU`, optional The BANDS table HDU. Default is None. format : str, optional FITS convention. Default is None. If None the format is guessed. The following formats are supported: - "gadf" - "fgst-ccube" - "fgst-ltcube" - "fgst-bexpcube" - "fgst-srcmap" - "fgst-template" - "fgst-srcmap-sparse" - "galprop" - "galprop2" - "healpy" Returns ------- hpx : `~HpxGeom` HEALPix geometry. """ifformatisNone:format=HpxConv.identify_hpx_format(header)conv=HPX_FITS_CONVENTIONS[format]axes=MapAxes.from_table_hdu(hdu_bands,format=format)ifheader["PIXTYPE"]!="HEALPIX":raiseValueError(f"Invalid header PIXTYPE: {header['PIXTYPE']} (must be HEALPIX)")ifheader["ORDERING"]=="RING":nest=Falseelifheader["ORDERING"]=="NESTED":nest=Trueelse:raiseValueError(f"Invalid header ORDERING: {header['ORDERING']} (must be RING or NESTED)")ifhdu_bandsisnotNoneand"NSIDE"inhdu_bands.columns.names:nside=hdu_bands.data.field("NSIDE").reshape(axes.shape).astype(int)elif"NSIDE"inheader:nside=header["NSIDE"]elif"ORDER"inheader:nside=2**header["ORDER"]else:raiseValueError("Failed to extract NSIDE or ORDER.")try:frame=coordsys_to_frame(header[conv.frame])exceptKeyError:frame=header.get("COORDSYS","icrs")try:region=header["HPX_REG"]exceptKeyError:try:region=header["HPXREGION"]exceptKeyError:region=Nonereturncls(nside,nest,frame=frame,region=region,axes=axes)
[docs]@classmethoddeffrom_hdu(cls,hdu,hdu_bands=None):"""Create an HPX object from a BinTable HDU. Parameters ---------- hdu : `~astropy.io.fits.BinTableHDU` The FITS HDU. hdu_bands : `~astropy.io.fits.BinTableHDU`, optional The BANDS table HDU. Default is None. Returns ------- hpx : `~HpxGeom` HEALPix geometry. """# FIXME: Need correct handling of IMPLICIT and EXPLICIT maps# if HPX region is not defined then geometry is defined by# the set of all pixels in the tableif"HPX_REG"notinhdu.header:pix=(hdu.data.field("PIX"),hdu.data.field("CHANNEL"))else:pix=Nonereturncls.from_header(hdu.header,hdu_bands=hdu_bands,pix=pix)
[docs]defto_header(self,format="gadf",**kwargs):"""Build and return FITS header for this HEALPix map."""header=fits.Header()format=kwargs.get("format",HPX_FITS_CONVENTIONS[format])# FIXME: For some sparse maps we may want to allow EXPLICIT# with an empty region stringindxschm=kwargs.get("indxschm",None)ifindxschmisNone:ifself._regionisNone:indxschm="IMPLICIT"elifself.is_regular==1:indxschm="EXPLICIT"else:indxschm="LOCAL"if"FGST"informat.convname.upper():header["TELESCOP"]="GLAST"header["INSTRUME"]="LAT"header[format.frame]=frame_to_coordsys(self.frame)header["PIXTYPE"]="HEALPIX"header["ORDERING"]=self.orderingheader["INDXSCHM"]=indxschmheader["ORDER"]=np.max(self.order)header["NSIDE"]=np.max(self.nside)header["FIRSTPIX"]=0header["LASTPIX"]=np.max(self.npix_max)-1header["HPX_CONV"]=format.convname.upper()ifself.frame=="icrs":header["EQUINOX"]=(2000.0,"Equinox of RA & DEC specifications")ifself.region:header["HPX_REG"]=self._regionreturnheader
[docs]@staticmethoddefget_index_list(nside,nest,region):"""Get list of pixels indices for all the pixels in a region. Parameters ---------- nside : int HEALPix NSIDE parameter. nest : bool Indexing scheme. If True, "NESTED" scheme. If False, "RING" scheme. region : str HEALPix region string. Returns ------- ilist : `~numpy.ndarray` List of pixel indices. """importhealpyashp# TODO: this should return something more friendly than a tuple# e.g. a namedtuple or a dicttokens=parse_hpxregion(region)reg_type=tokens[0]ifreg_type=="DISK":lon,lat=float(tokens[1]),float(tokens[2])radius=np.radians(float(tokens[3]))vec=coords_to_vec(lon,lat)[0]ilist=hp.query_disc(nside,vec,radius,inclusive=False,nest=nest)elifreg_type=="DISK_INC":lon,lat=float(tokens[1]),float(tokens[2])radius=np.radians(float(tokens[3]))vec=coords_to_vec(lon,lat)[0]fact=int(tokens[4])ilist=hp.query_disc(nside,vec,radius,inclusive=True,nest=nest,fact=fact)elifreg_type=="HPX_PIXEL":nside_pix=int(tokens[2])iftokens[1]=="NESTED":ipix_ring=hp.nest2ring(nside_pix,int(tokens[3]))eliftokens[1]=="RING":ipix_ring=int(tokens[3])else:raiseValueError(f"Invalid ordering scheme: {tokens[1]!r}")ilist=match_hpx_pix(nside,nest,nside_pix,ipix_ring)else:raiseValueError(f"Invalid region type: {reg_type!r}")returnilist
@propertydefwidth(self):"""Width of the map."""# TODO: simplifyimporthealpyashpifself.is_allsky:width=180.0elifself.region=="explicit":idx=unravel_hpx_index(self._ipix,self.npix_max)nside=self._get_nside(idx)ang=hp.pix2ang(nside,idx[0],nest=self.nest,lonlat=True)dirs=SkyCoord(ang[0],ang[1],unit="deg",frame=self.frame)width=np.max(dirs.separation(self.center_skydir))else:tokens=parse_hpxregion(self.region)iftokens[0]in{"DISK","DISK_INC"}:width=float(tokens[3])eliftokens[0]=="HPX_PIXEL":pix_size=get_pix_size_from_nside(int(tokens[2]))width=2.0*pix_sizereturnu.Quantity(width,"deg")def_get_nside(self,idx):ifself.nside.size>1:returnself.nside[tuple(idx[1:])]else:returnself.nside
[docs]defto_wcs_geom(self,proj="AIT",oversample=2,width_pix=None):"""Make a WCS projection appropriate for this HEALPix pixelization. Parameters ---------- proj : str, optional Projection type of WCS geometry. Default is "AIT". oversample : float, optional Oversampling factor for WCS map. This will be the approximate ratio of the width of a HEALPix pixel to a WCS pixel. If this parameter is None then the width will be set from ``width_pix``. Default is 2. width_pix : int, optional Width of the WCS geometry in pixels. The pixel size will be set to the number of pixels satisfying ``oversample`` or ``width_pix`` whichever is smaller. If this parameter is None then the width will be set from ``oversample``. Default is None. Returns ------- wcs : `~gammapy.maps.WcsGeom` WCS geometry. """fromgammapy.mapsimportWcsGeompix_size=get_pix_size_from_nside(self.nside)binsz=np.min(pix_size)/oversamplewidth=2.0*self.width.to_value("deg")+np.max(pix_size)ifwidth_pixisnotNoneandint(width/binsz)>width_pix:binsz=width/width_pixifwidth>90.0:width=min(360.0,width),min(180.0,width)axes=copy.deepcopy(self.axes)returnWcsGeom.create(width=width,binsz=binsz,frame=self.frame,axes=axes,skydir=self.center_skydir,proj=proj,)
[docs]defto_wcs_tiles(self,nside_tiles=4,margin="0 deg"):"""Create WCS tiles geometries from HPX geometry with given nside. The HEALPix geom is divide into superpixels defined by ``nside_tiles``, which are then represented by a WCS geometry using a tangential projection. The number of WCS tiles is given by the number of pixels for the given ``nside_tiles``. Parameters ---------- nside_tiles : int, optional HEALPix NSIDE parameter for super pixel tiles. Default is 4. margin : `~astropy.units.Quantity`, optional Width margin of the WCS tile. Default is "0 deg". Returns ------- wcs_tiles : list List of WCS tile geometries. """importhealpyashpfromgammapy.mapsimportWcsGeommargin=u.Quantity(margin)ifnside_tiles>=self.nside:raiseValueError(f"nside_tiles must be < {self.nside}")ifnotself.is_allsky:raiseValueError("to_wcs_tiles() is only supported for all sky geoms")binsz=np.degrees(hp.nside2resol(self.nside))*u.deghpx=self.to_image().to_nside(nside=nside_tiles)wcs_tiles=[]forpixinrange(int(hpx.npix[0])):skydir=hpx.pix_to_coord([pix])vtx=hp.boundaries(nside=hpx.nside.item(),pix=pix,nest=hpx.nest,step=1)lon,lat=hp.vec2ang(vtx.T,lonlat=True)boundaries=SkyCoord(lon*u.deg,lat*u.deg,frame=hpx.frame)# Compute maximum separation between all pairs of boundaries and take it# as widthwidth=boundaries.separation(boundaries[:,np.newaxis]).max()wcs_tile_geom=WcsGeom.create(skydir=(float(skydir[0].item()),float(skydir[1].item())),width=width+margin,binsz=binsz,frame=hpx.frame,proj="TAN",axes=self.axes,)wcs_tiles.append(wcs_tile_geom)returnwcs_tiles
[docs]defget_idx(self,idx=None,local=False,flat=False,sparse=False,mode="center",axis_name=None,):# TODO: simplify this!!!ifidxisnotNoneandnp.any(np.array(idx)>=np.array(self.shape_axes)):raiseValueError(f"Image index out of range: {idx!r}")# Regular all- and partial-sky mapsifself.is_regular:pix=[np.arange(np.max(self._npix))]ifidxisNone:foraxinself.axes:ifmode=="edges"andax.name==axis_name:pix+=[np.arange(-0.5,ax.nbin,dtype=float)]else:pix+=[np.arange(ax.nbin,dtype=int)]else:pix+=[tfortinidx]pix=np.meshgrid(*pix[::-1],indexing="ij",sparse=sparse)[::-1]pix=self.local_to_global(pix)# Non-regular all-skyelifself.is_allskyandnotself.is_regular:shape=(np.max(self.npix),)ifidxisNone:shape=shape+self.shape_axeselse:shape=shape+(1,)*len(self.axes)pix=[np.full(shape,-1,dtype=int)foriinrange(1+len(self.axes))]foridx_imginnp.ndindex(self.shape_axes):ifidxisnotNoneandidx_img!=idx:continuenpix=self._npix[idx_img]ifidxisNone:s_img=(slice(0,npix),)+idx_imgelse:s_img=(slice(0,npix),)+(0,)*len(self.axes)pix[0][s_img]=np.arange(self._npix[idx_img])forjinrange(len(self.axes)):pix[j+1][s_img]=idx_img[j]pix=[p.Tforpinpix]# Explicit pixel indiceselse:ifidxisnotNone:npix_sum=np.concatenate(([0],np.cumsum(self._npix)))idx_ravel=np.ravel_multi_index(idx,self.shape_axes)s=slice(npix_sum[idx_ravel],npix_sum[idx_ravel+1])else:s=slice(None)pix_flat=unravel_hpx_index(self._ipix[s],self.npix_max)shape=(np.max(self.npix),)ifidxisNone:shape=shape+self.shape_axeselse:shape=shape+(1,)*len(self.axes)pix=[np.full(shape,-1,dtype=int)for_inrange(1+len(self.axes))]foridx_imginnp.ndindex(self.shape_axes):ifidxisnotNoneandidx_img!=idx:continuenpix=int(self._npix[idx_img].item())ifidxisNone:s_img=(slice(0,npix),)+idx_imgelse:s_img=(slice(0,npix),)+(0,)*len(self.axes)ifself.axes:m=np.all(np.stack([pix_flat[i+1]==tfori,tinenumerate(idx_img)]),axis=0,)pix[0][s_img]=pix_flat[0][m]else:pix[0][s_img]=pix_flat[0]forjinrange(len(self.axes)):pix[j+1][s_img]=idx_img[j]pix=[p.Tforpinpix]iflocal:pix=self.global_to_local(pix)ifflat:pix=tuple([p[p!=INVALID_INDEX.int]forpinpix])returnpix
[docs]defregion_mask(self,regions):"""Create a mask from a given list of regions. The mask is filled such that a pixel inside the region is filled with "True". To invert the mask, e.g. to create a mask with exclusion regions the tilde (~) operator can be used (see example below). Parameters ---------- regions : str, `~regions.Region` or list of `~regions.Region` Region or list of regions (pixel or sky regions accepted). A region can be defined as a string ind DS9 format as well. See http://ds9.si.edu/doc/ref/region.html for details. Returns ------- mask_map : `~gammapy.maps.WcsNDMap` of boolean type Boolean region mask. """fromgammapy.mapsimportMap,RegionGeomifnotself.is_regular:raiseValueError("Multi-resolution maps not supported yet")# TODO: use spatial coordinates only...geom=RegionGeom.from_regions(regions)coords=self.get_coord()mask=geom.contains(coords)returnMap.from_geom(self,data=mask)
[docs]defget_coord(self,idx=None,flat=False,sparse=False,mode="center",axis_name=None):ifmode=="edges"andaxis_nameisNone:raiseValueError("Mode 'edges' requires axis name to be defined")pix=self.get_idx(idx=idx,flat=flat,sparse=sparse,mode=mode,axis_name=axis_name)data=self.pix_to_coord(pix)coords=MapCoord.create(data=data,frame=self.frame,axis_names=self.axes.names)returncoords
[docs]defsolid_angle(self):"""Solid angle array as a `~astropy.units.Quantity` in ``sr``. The array has the same dimensionality as ``map.nside`` as all pixels have the same solid angle. """importhealpyashpreturnQuantity(hp.nside2pixarea(self.nside),"sr")
[docs]defis_allclose(self,other,rtol_axes=1e-6,atol_axes=1e-6):"""Compare two data IRFs for equivalency. Parameters ---------- other : `HpxGeom` Geometry to compare against. rtol_axes : float, optional Relative tolerance for axes comparison. Default is 1e-6. atol_axes : float, optional Relative tolerance for 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.is_allskyandnotother.is_allsky:returnFalseifself.data_shape!=other.data_shape:returnFalseaxes_eq=self.axes.is_allclose(other.axes,rtol=rtol_axes,atol=atol_axes)hpx_eq=(self.nside==other.nsideandself.frame==other.frameandself.order==other.orderandself.nest==other.nest)returnaxes_eqandhpx_eq