black

Author: martinkilbinger

cs_util/cat.py

@@ -11,6 +11,7 @@
 import os
 from datetime import datetime
 from importlib.metadata import version
+import numpy as np
 import healpy as hp
 
 from astropy.io import fits
@@ -271,7 +272,27 @@ def read_dndz(file_path):
 
 
 def read_hp_mask(input_name, verbose=False):
+    """Read Hp Mask.
 
+    Read healpy mask.
+
+    Parameters
+    ----------
+    input_name : str
+        input file name
+    verbose : bool, optional
+        verbose output if ``True``;; default is ``False``
+
+    Returns
+    -------
+    array
+        mask values
+    bool
+        nest value (always ``False``)
+    int
+        nside
+
+    """
     if verbose:
         print(f'Reading mask {input_name}...')
 
@@ -305,3 +326,45 @@ def read_hp_mask(input_name, verbose=False):
         raise KeyError('NSIDE not found in FITS mask header')
 
     return mask, nest, nside
+
+
+def get_binned_area(ra, dec, nside=512, return_pix=False):                                         
+    """Get Binned Area.                                                          
+                                                                                 
+    Return sky area corresponding to occupied pixels of binned catalogue.        
+
+    Parameters
+    ----------
+    ra : np.ndarray
+        Right Ascension coordinates
+    dec : np.ndarray
+        Declination coordinates
+    nside : int, optional
+        nside for binning; default is 512
+    return_pix: bool, optional
+        return valid pixel list if ``True``; default is ``False``
+                                                                                 
+    Returns                                                                      
+    -------                                                                      
+    float                                                                        
+        area in square degrees                                                   
+    list
+        pixels of input data positions
+                                                                                 
+    """                                                                          
+    # Pixel list of input data                                                   
+    pix = hp.ang2pix(nside, ra, dec, lonlat=True)                               
+                                                                                 
+    # Number of occupied pixels                                                  
+    Nocc = np.unique(pix).size                                                  
+                                                                                 
+    # Pixel area in square degrees                                               
+    pix_area_deg2 = hp.nside2pixarea(nside, degrees=True)                        
+                                                                                 
+    # Footprint area in square degrees                                           
+    area_deg2 = Nocc * pix_area_deg2                                             
+
+    if return_pix:
+        return area_deg2, pix
+    else:
+        return area_deg2

cs_util/plots.py

@@ -61,7 +61,20 @@ def show():
     plt.close()
 
 
-def dx(idx, nx=3, fx=1.015, log=True):
+def get_x_dx(x_arr, shift_x, idx, log=True):
+
+    if shift_x:
+        if log:
+            this_x = x_arr[idx] * dx(idx, len(x_arr), log=log)
+        else:
+            raise ValueError("shift_x without log not implemented yet")
+    else:
+        this_x = x_arr[idx]
+
+    return this_x
+
+
+def dx(idx, nx=3, fx=1.03, log=True):
     """Dx.
 
     Return small shift useful to diplace points along the the x-axis
@@ -288,12 +301,8 @@ def plot_data_1d(
         fig = ax.figure
 
     for idx in range(len(x)):
-        this_x = x[idx]
-        if shift_x:
-            if xlog:
-                this_x *= dx(idx, len(x), log=xlog)
-            else:
-                raise ValueError("shift_x without log not implemented yet")
+        this_x = get_x_dx(x, shift_x, idx, log=xlog)
+
         if np.isnan(yerr[idx]).all():
             eb = ax.plot(
                 this_x,