Merge pull request #63 from CosmoStat/develop

v0.1.9

Author: martinkilbinger

cs_util/args.py

@@ -9,6 +9,7 @@
 
 import sys
 import os
+import shlex
 
 from optparse import OptionParser
 
@@ -103,6 +104,9 @@ def my_string_split(string, num=-1, verbose=False, stop=False, sep=None):
     the first in the list [space, underscore] that occurs in the string.
     (Thus, if both occur, use space.)
 
+    Handles quoted strings: entries containing spaces can be enclosed
+    in double quotes, e.g., 'value1 "entry with spaces" value3'.
+
     Parameters
     ----------
     string : str
@@ -131,6 +135,19 @@ def my_string_split(string, num=-1, verbose=False, stop=False, sep=None):
     if string is None:
         return None
 
+    # Handle quoted strings with shlex
+    if '"' in string or "'" in string:
+        try:
+            res = shlex.split(string)
+            if num != -1 and num != len(res) and stop:
+                raise ValueError(
+                    f"String has {len(res)} elements, required is {num}"
+                )
+            return res
+        except ValueError:
+            # Fall through to regular splitting if shlex fails
+            pass
+
     if sep is None:
         has_space = string.find(" ")
         has_underscore = string.find("_")
@@ -162,7 +179,7 @@ def my_string_split(string, num=-1, verbose=False, stop=False, sep=None):
 
     if num != -1 and num != len(res) and stop:
         raise ValueError(
-            f"String '{len(res)}' has length {num}, required is {num}"
+            f"String has {len(res)} elements, required is {num}"
         )
 
     return res

cs_util/cat.py

@@ -10,6 +10,7 @@
 
 import os
 from datetime import datetime
+
 from importlib.metadata import version
 import numpy as np
 import healpy as hp
@@ -253,118 +254,93 @@ def read_dndz(file_path):
 
     Returns
     -------
-    list :
+    np.array
         redshift bin centers
-    list :
+    np.array
         number densities
-    list :
-        redshift bin edges
+    np.array
+        redshift bin edges; one less than centers and density arrays
 
     """
-    dat = ascii.read(file_path, format="commented_header")
+    try:
+        # Expecting header line "# z dn_dz"
+        dat = ascii.read(file_path, format="commented_header")
+        missing = [col for col in ("z", "dn_dz") if col not in dat.dtype.names]
+        if missing:
+            raise ValueError(
+                f"Missing columns in dndz path {file_path}: {missing}"
+            )
+    except:
+        # No header line
+        dat = ascii.read(file_path)
+        dat.rename_column("col1", "z")
+        dat.rename_column("col2", "dn_dz")
+
+    # Remove last n(z) value which should be zero, to match bin centers
+    tolerance = 1e-5
+    if dat["dn_dz"][-1] / sum(dat["dn_dz"]) > tolerance:
+        raise ValueError("dn_dz at last z-edge = {dat['dn_dz'][-1]}, no zero")
 
-    # Remove last n(z) value which is zero, to match bin centers
     nz = dat["dn_dz"][:-1]
     z_edges = dat["z"]
+
     z_centers = bin_edges2centers(z_edges)
 
     return z_centers, nz, z_edges
 
 
-def read_hp_mask(input_name, verbose=False):
+def read_hp_mask(input_path, verbose=False):
     """Read Hp Mask.
 
-    Read healpy mask.
+    Read healpix mask FITS file.
 
-    Parameters
-    ----------
-    input_name : str
-        input file name
-    verbose : bool, optional
-        verbose output if ``True``;; default is ``False``
+        Parameters
+        ----------
+        input_path : str
+                input file path
+        verbose : bool, optional
+                verbose output if ``True``; default is ``False``
 
-    Returns
-    -------
-    array
-        mask values
-    bool
-        nest value (always ``False``)
-    int
-        nside
+        Returns
+        -------
+        array
+                mask information
+        bool
+                NEST (RING) ordering if ``True`` (``False``)
+        int
+                nside
 
     """
     if verbose:
-        print(f'Reading mask {input_name}...')
+        print(f"Reading mask {input_path}...")
 
     nest = False
 
-    # Open input mask                                                           
-    mask, header= hp.read_map(
-        input_name,
+    # Open input mask
+    mask, header = hp.read_map(
+        input_path,
         h=True,
         nest=nest,
     )
-    for (key, value) in header:
-        if key == 'ORDERING':
-            if value == 'RING':
+    for key, value in header:
+        if key == "ORDERING":
+            if value == "RING":
                 if nest:
                     raise ValueError(
-                        'input mask has ORDENING=RING, set nest to False'
+                        "input mask has ORDENING=RING, set nest to False"
                     )
-            elif value == 'NEST':
+            elif value == "NEST":
                 if not nest:
                     raise ValueError(
-                        'input mask has ORDENING=NEST, set nest to True'
+                        "input mask has ORDENING=NEST, set nest to True"
                     )
 
-    # Get nside from header                                                     
+    # Get nside from header
     nside = None
     for key, value in header:
-        if key == 'NSIDE':
+        if key == "NSIDE":
             nside = int(value)
     if not nside:
-        raise KeyError('NSIDE not found in FITS mask header')
+        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