import healpy as hp
import numpy as np
import matplotlib.pyplot as plt
# from scipy.special import gamma as gafun
# from scipy.stats import norm
# import scipy.integrate as integrate
import os
# import warnings
import subprocess
import pandas as pd
# import numpy.ma as ma
#def findmaxima
def _readmax_forhotspot(file,folder=''):
'''Read the output file generated by hotspot'''
data=np.loadtxt(f'{folder}{file}')
toreturn=pd.DataFrame(data=np.array(data),columns=['Pixel number','Intensity'])
toreturn['Pixel number']=toreturn['Pixel number'].astype('int')
return(toreturn.sort_values(by=['Intensity'],ascending=False))
[docs]def hotspot(mapp, minima=False, folder=None, name='Map', hotspot_command='hotspot'):
'''Find the maxima (or minima) in a map.
Notes
----
This function use the Healpix utility ``hotspot``, which is not yet implemented in Python distribution of Healpix.
The utility ``healpix`` is called on a terminal and it writes files on the hard drive. This function deletes these
files after the are read by Python unless a specific folder is given as an input. This function requires
``healpix`` to be installed on the computer for now, although this is expected to change soon.
The utility has several limitations regarding the intensity of the maxima. It only stores 4 digits before the
period and 6 digits after it. It is highly recommended to normalise the map before using this function. For extreme
cases, outliers have to be removed before applying the function, for example with::
>>> mapp[mapp>9999]=9999
Parameters
----------
mapp : np.ndarray
A Healpy map of the sky.
minima : bool, optional
If ``True``, this function return the minima instead of the maxima. We remind the user that the rest of the functions
are designed for set of maxima. If minima are to be used, we recommend changing the sign of the minima so the can be
regarded as maxima of the **opposite** map.
folder : str, optional
Use only if the files generated by ``hotspot`` are to be stored. This argument indicates where.
name : str, optional
The name of the map to use in the written files. Unnecesary unless ``folder`` is given.
hotspot_command : str, optional
The name of the command installed on the computer, usually 'hotspot' or 'hotspot_gcc'.
Returns
-------
pd.DataFrame
Pandas DataFrame with two columns: 'Pixel number' denotes the number of the pixel where the maximum is detected (depends
on the ``nside`` of ``mapp``); 'Intensity' denotes the value of this maxima in the input map
'''
rem=False
if folder is None:
rem=True
folder=os.getcwd()
if folder[-1] != '/':
folder=folder+'/'
hp.write_map(folder+name+'.fits',mapp);
hotspotfile= open(f'{folder}{name}.par','w')
hotspotfile.writelines([f'infile={folder}{name}.fits \n',
f'extrema_outfile = {folder}pixlminmax_{name}.fits \n',
f'maxima_outfile = {folder}maxima_{name}.dat \n',
f'minima_outfile = {folder}minima_{name}.dat'])
hotspotfile.close()
subprocess.run([f'{hotspot_command} {folder}{name}.par > /dev/null'],check=True,shell=True)
# !hotspot {name}.par > /dev/null
if minima != True:
toreturn=_readmax_forhotspot(f'{folder}maxima_{name}.dat') #maxima
else:
toreturn=_readmax_forhotspot(f'{folder}minima_{name}.dat') #minima
if rem: #This is True is there is not an input folder
os.remove(f'{name}.fits')
os.remove(f'{name}.par')
os.remove(f'maxima_{name}.dat')
os.remove(f'pixlminmax_{name}.fits')
os.remove(f'minima_{name}.dat')
return(toreturn)
[docs]def locate_maxima(maxima,nside,lonlat=True):
'''Find the location of the input maxima.
Parameters
----------
maxima : pd.DataFrame
Pandas DataFrame with the information of the maxima. Has to contain at least the column 'Pixel number'.
nside : int
An integer corresponding to the ``nside`` of the map where the maxima where detected. It needs to be a power of 2.
lonlat : bool, optional
If ``True``, the functions returns the longitude and latitude, in degrees. If ``False``, the function returns
colatitude and longitude, in radians.
Returns
-------
pd.DataFrame
The same Pandas DataFrame as the input ``maxima``, with two new columns denoting the position of each maxima on the sky.
'''
if lonlat==False:
[name1,name2]=['Colatitude [rad]','Longitude [rad]']
elif lonlat==True:
[name1,name2]=['Longitude [deg]', 'Latitude[deg]']
maxima[name1],maxima[name2]=hp.pix2ang(nside,maxima['Pixel number'],lonlat=lonlat)
return(maxima)
[docs]def max_inmask(maxima,mask,withvalue=1.):
'''Get the maxima located where the mask have a specific value (e.g., 1)
Parameters
----------
maxima : pd.DataFrame
Pandas DataFrame with the information of the maxima. Has to contain at least the column 'Pixel number'.
mask : np.ndarray
A Healpix map to be used as mask. **Its ``nside`` should be the same as the original map.**
withvalue : float, optional
Value to be matched. Only maxima that have this value on the ``mask`` will be returned.
Returns
-------
pd.DataFrame
The same Pandas DataFrame as the input ``maxima``, but containing only the maxima located at points where the
``mask`` has the value ``withvalue`` (1 by default).
'''
filtered_maxima=maxima.loc[(mask[maxima['Pixel number']]==withvalue)]
return(filtered_maxima)
[docs]def max_threshold(maxima,threshold):
'''Get the maxima with intensity above or equal to a threshold
Parameters
----------
maxima : pd.DataFrame
Pandas DataFrame with the information of the maxima. Has to contain at least the column 'Intensity'.
threshold : float
Minimum value of the maxima to be returned.
Returns
-------
pd.DataFrame
The same Pandas DataFrame as the input ``maxima``, but containing only the maxima with an intensity equal
or greater than ``threshold``.
'''
return(maxima[maxima['Intensity']>=threshold])
[docs]def plot_maxima(maxima,f=None,bins=50):
'''Visualise the maxima distribution and compare it with a theoretical one.
Parameters
----------
maxima : pd.DataFrame
Pandas DataFrame with the information of the maxima. Has to contain at least the column 'Intensity'.
f : function or None, optional
Theoretical distribution of the maxima, f.
bins : int, optional
Number of bins for the histogram of the maxima
Returns
-------
matplotlib.figure.Figure
Figure containing the histogram of the maxima intensity and, if included, the expected distribution from ``f``.
'''
fig=plt.figure()
num,xbins,_=plt.hist(maxima['Intensity'],bins=bins,density=1,label='Maxima distribution of the map')
if f is not None:
xvec=np.linspace(xbins[0],xbins[-1],num=bins*4)
plt.plot(xvec,f(xvec),label='Expected distribution')
plt.legend()
plt.tight_layout()
ax=plt.gca()
return(fig,ax)