|
10 | 10 |
|
11 | 11 | import os |
12 | 12 | from datetime import datetime |
| 13 | + |
13 | 14 | from importlib.metadata import version |
| 15 | +import numpy as np |
| 16 | +import healpy as hp |
14 | 17 |
|
15 | 18 | from astropy.io import fits |
16 | 19 | from astropy.io import ascii |
@@ -251,19 +254,93 @@ def read_dndz(file_path): |
251 | 254 |
|
252 | 255 | Returns |
253 | 256 | ------- |
254 | | - list : |
| 257 | + np.array |
255 | 258 | redshift bin centers |
256 | | - list : |
| 259 | + np.array |
257 | 260 | number densities |
258 | | - list : |
259 | | - redshift bin edges |
| 261 | + np.array |
| 262 | + redshift bin edges; one less than centers and density arrays |
260 | 263 |
|
261 | 264 | """ |
262 | | - dat = ascii.read(file_path, format="commented_header") |
| 265 | + try: |
| 266 | + # Expecting header line "# z dn_dz" |
| 267 | + dat = ascii.read(file_path, format="commented_header") |
| 268 | + missing = [col for col in ("z", "dn_dz") if col not in dat.dtype.names] |
| 269 | + if missing: |
| 270 | + raise ValueError( |
| 271 | + f"Missing columns in dndz path {file_path}: {missing}" |
| 272 | + ) |
| 273 | + except: |
| 274 | + # No header line |
| 275 | + dat = ascii.read(file_path) |
| 276 | + dat.rename_column("col1", "z") |
| 277 | + dat.rename_column("col2", "dn_dz") |
| 278 | + |
| 279 | + # Remove last n(z) value which should be zero, to match bin centers |
| 280 | + tolerance = 1e-5 |
| 281 | + if dat["dn_dz"][-1] / sum(dat["dn_dz"]) > tolerance: |
| 282 | + raise ValueError("dn_dz at last z-edge = {dat['dn_dz'][-1]}, no zero") |
263 | 283 |
|
264 | | - # Remove last n(z) value which is zero, to match bin centers |
265 | 284 | nz = dat["dn_dz"][:-1] |
266 | 285 | z_edges = dat["z"] |
| 286 | + |
267 | 287 | z_centers = bin_edges2centers(z_edges) |
268 | 288 |
|
269 | 289 | return z_centers, nz, z_edges |
| 290 | + |
| 291 | + |
| 292 | +def read_hp_mask(input_path, verbose=False): |
| 293 | + """Read Hp Mask. |
| 294 | +
|
| 295 | + Read healpix mask FITS file. |
| 296 | +
|
| 297 | + Parameters |
| 298 | + ---------- |
| 299 | + input_path : str |
| 300 | + input file path |
| 301 | + verbose : bool, optional |
| 302 | + verbose output if ``True``; default is ``False`` |
| 303 | +
|
| 304 | + Returns |
| 305 | + ------- |
| 306 | + array |
| 307 | + mask information |
| 308 | + bool |
| 309 | + NEST (RING) ordering if ``True`` (``False``) |
| 310 | + int |
| 311 | + nside |
| 312 | +
|
| 313 | + """ |
| 314 | + if verbose: |
| 315 | + print(f"Reading mask {input_path}...") |
| 316 | + |
| 317 | + nest = False |
| 318 | + |
| 319 | + # Open input mask |
| 320 | + mask, header = hp.read_map( |
| 321 | + input_path, |
| 322 | + h=True, |
| 323 | + nest=nest, |
| 324 | + ) |
| 325 | + for key, value in header: |
| 326 | + if key == "ORDERING": |
| 327 | + if value == "RING": |
| 328 | + if nest: |
| 329 | + raise ValueError( |
| 330 | + "input mask has ORDENING=RING, set nest to False" |
| 331 | + ) |
| 332 | + elif value == "NEST": |
| 333 | + if not nest: |
| 334 | + raise ValueError( |
| 335 | + "input mask has ORDENING=NEST, set nest to True" |
| 336 | + ) |
| 337 | + |
| 338 | + # Get nside from header |
| 339 | + nside = None |
| 340 | + for key, value in header: |
| 341 | + if key == "NSIDE": |
| 342 | + nside = int(value) |
| 343 | + if not nside: |
| 344 | + raise KeyError("NSIDE not found in FITS mask header") |
| 345 | + |
| 346 | + return mask, nest, nside |
0 commit comments