Inspect WebSky high-flux catalog output

import h5py
import numpy as np
import healpy as hp
import matplotlib.pyplot as plt

freqs = [
    "18.7",
    "24.5",
    "44.0",
    "70.0",
    "100.0",
    "143.0",
    "217.0",
    "353.0",
    "545.0",
    "643.0",
    "729.0",
    "857.0",
    "906.0",
]

cat = h5py.File("data/matched_catalogs_2/catalog_100.0.h5", "r")

There are no metadata in the file, I guess fluxes are in Jy

cutoff_flux = 1e-3

high_flux_sources_mask = cat["flux"][:] > cutoff_flux

(high_flux_sources_mask).sum()

372255

high_flux_sources_mask.mean() * 100

0.13211945911740433

for k, v in cat.items():
    print(k, v[:3])

flux [3.24291534e-07 3.16862867e-07 3.17171157e-07]
phi [3.22861886 3.22861886 3.22861886]
polarized flux [1.42910628e-09 1.99535624e-08 2.29563857e-09]
theta [1.64009452 1.64009452 1.64009452]

(all_indices,) = np.nonzero(high_flux_sources_mask)

len(all_indices)

372255

all_indices = np.array(sorted(all_indices))

all_indices contains the index of the sources in the original Websky catalog

all_indices

array([    11253,     16428,     24110, ..., 281755117, 281755430,
       281755795])

import pandas as pd
import xarray as xr

columns = ["theta", "phi", "flux", "polarized flux"]

flux = xr.DataArray(
    data=np.zeros((len(all_indices), len(freqs)), dtype=np.float64),
    coords={"index": all_indices, "freq": list(map(float, freqs))},
    name="flux",
)
fluxnorm = flux.copy()

polarized_flux = flux.copy()

sources_xr = xr.Dataset(
    {"flux": flux, "polarized_flux": polarized_flux, "fluxnorm": fluxnorm}
)
for freq in freqs:
    print(freq)
    cat = h5py.File(f"data/matched_catalogs_2/catalog_{freq}.h5", "r")
    for column in ["flux", "polarized_flux"]:
        sources_xr[column].loc[dict(index=all_indices, freq=float(freq))] = cat[
            column.replace("_", " ")
        ][high_flux_sources_mask]

18.7
24.5
44.0
70.0
100.0
143.0
217.0
353.0
545.0
643.0
729.0
857.0
906.0

sources_xr

<xarray.Dataset> Size: 119MB
Dimensions:         (index: 372255, freq: 13)
Coordinates:
  * index           (index) int64 3MB 11253 16428 24110 ... 281755430 281755795
  * freq            (freq) float64 104B 18.7 24.5 44.0 ... 729.0 857.0 906.0
Data variables:
    flux            (index, freq) float64 39MB 0.02854 0.01962 ... 0.001688
    polarized_flux  (index, freq) float64 39MB 0.001121 0.0004986 ... 3.246e-05
    fluxnorm        (index, freq) float64 39MB 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0

xarray.Dataset

• Dimensions:
  • index: 372255
  • freq: 13
• Coordinates: (2)
  • index
    (index)
    int64
    11253 16428 ... 281755430 281755795

    array([    11253,     16428,     24110, ..., 281755117, 281755430, 281755795])

  • freq
    (freq)
    float64
    18.7 24.5 44.0 ... 857.0 906.0

    array([ 18.7,  24.5,  44. ,  70. , 100. , 143. , 217. , 353. , 545. , 643. ,
           729. , 857. , 906. ])

• Data variables: (3)
  • flux
    (index, freq)
    float64
    0.02854 0.01962 ... 0.001688

    array([[0.02853885, 0.01962042, 0.0104193 , ..., 0.00049761, 0.00043035,
            0.00040939],
           [0.07462281, 0.04701015, 0.01973497, ..., 0.00036078, 0.00029885,
            0.00028011],
           [0.00986449, 0.00733716, 0.00481509, ..., 0.0005461 , 0.00049098,
            0.00047334],
           ...,
           [0.00131828, 0.00119997, 0.00135572, ..., 0.0014146 , 0.00140503,
            0.00140175],
           [0.00196843, 0.00173202, 0.00178621, ..., 0.00128392, 0.00125407,
            0.00124397],
           [0.00101562, 0.00095178, 0.00116291, ..., 0.00167089, 0.00168358,
            0.00168796]])

  • polarized_flux
    (index, freq)
    float64
    0.001121 0.0004986 ... 3.246e-05

    array([[1.12068269e-03, 4.98555717e-04, 7.24046113e-05, ...,
            3.11687218e-05, 2.45648989e-05, 3.17132094e-06],
           [9.41744819e-03, 2.96976970e-04, 1.30501488e-04, ...,
            3.22987471e-05, 5.65086248e-06, 1.69376062e-06],
           [7.76644811e-05, 6.61808008e-06, 9.53326016e-05, ...,
            1.76184112e-05, 6.33791444e-07, 1.04503160e-05],
           ...,
           [4.60562251e-06, 1.46073426e-05, 2.50197045e-04, ...,
            3.62357605e-05, 5.36425478e-06, 6.18042104e-05],
           [5.87950126e-06, 1.00558582e-05, 1.04939129e-04, ...,
            1.93261112e-05, 1.08091372e-05, 5.17453154e-05],
           [5.70268412e-06, 6.20402307e-06, 1.40497898e-04, ...,
            1.50881197e-05, 3.55578559e-05, 3.24648208e-05]])

  • fluxnorm
    (index, freq)
    float64
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    array([[0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           ...,
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.],
           [0., 0., 0., ..., 0., 0., 0.]])

• Indexes: (2)
  • index
    PandasIndex

    PandasIndex(Index([    11253,     16428,     24110,     25342,     26268,     27427,
               27855,     40649,     43434,     44634,
           ...
           281751967, 281752822, 281753210, 281753640, 281753799, 281754105,
           281754532, 281755117, 281755430, 281755795],
          dtype='int64', name='index', length=372255))

  • freq
    PandasIndex

    PandasIndex(Index([ 18.7,  24.5,  44.0,  70.0, 100.0, 143.0, 217.0, 353.0, 545.0, 643.0,
           729.0, 857.0, 906.0],
          dtype='float64', name='freq'))

• Attributes: (0)

min_idx = int(sources_xr.flux.sel(freq=100).argmin().item())
max_idx = int(sources_xr.flux.sel(freq=100).argmax().item())
print("Index with minimum flux at 100 GHz:", min_idx)
print("Index with maximum flux at 100 GHz:", max_idx)

Index with minimum flux at 100 GHz: 64456
Index with maximum flux at 100 GHz: 9474

sources_xr = sources_xr.sortby(sources_xr.flux.sel(freq=100.0), ascending=False)

len(sources_xr)

3

max_loc = int(sources_xr.flux.sel(freq=100).argmax().item())
min_loc = int(sources_xr.flux.sel(freq=100).argmin().item())
print("Location (position) of maximum flux at 100 GHz:", max_loc)
print("Location (position) of minimum flux at 100 GHz:", min_loc)

Location (position) of maximum flux at 100 GHz: 0
Location (position) of minimum flux at 100 GHz: 372254

sources_xr.coords["index"]

<xarray.DataArray 'index' (index: 372255)> Size: 3MB
array([ 27897657, 242175600, 270516228, ..., 246732122, 263347056, 166742332])
Coordinates:
  * index    (index) int64 3MB 27897657 242175600 ... 263347056 166742332

xarray.DataArray

'index'

• index: 372255

• 27897657 242175600 270516228 ... 246732122 263347056 166742332

  array([ 27897657, 242175600, 270516228, ..., 246732122, 263347056, 166742332])

• Coordinates: (1)
  • index
    (index)
    int64
    27897657 242175600 ... 166742332

    array([ 27897657, 242175600, 270516228, ..., 246732122, 263347056, 166742332])

• Indexes: (1)
  • index
    PandasIndex

    PandasIndex(Index([ 27897657, 242175600, 270516228, 276407774, 256213631, 281544643,
           263719520, 259533227, 251679008, 252717495,
           ...
           238190295, 145826702, 108046020, 232265270, 174768449, 256408684,
           236382968, 246732122, 263347056, 166742332],
          dtype='int64', name='index', length=372255))

• Attributes: (0)

# this was needed for some missing functionality in an older version of xarray
#sources_xr.coords["index"] = np.arange(len(sources_xr.coords["index"]))

sources_xr["fluxnorm"] = sources_xr["flux"] / sources_xr["flux"].sel(freq=100)

#print(sources_xr["fluxnorm"].loc[dict(index=s)], sources_xr["flux"].loc[dict(index=s)])

#sources_xr.fluxnorm.plot(vmin=0, vmax=100)
#plt.figure()
#sources_xr.flux.plot(vmin=0, vmax=100)

sources_xr["logpolycoefflux"] = xr.DataArray(
    np.zeros((len(all_indices), 5), dtype=np.float64),
    dims=["index", "power"],
    coords={"power": np.arange(5)[::-1]},
)
sources_xr["logpolycoefnorm"] = sources_xr["logpolycoefflux"].copy()
sources_xr["logpolycoefpolflux"] = sources_xr["logpolycoefflux"].copy()

from scipy.optimize import curve_fit
import time

def is_interactive():
    try:
        from IPython import get_ipython
        return get_ipython() is not None
    except ImportError:
        return False

def model(freq, a, b, c, d, e):
    log_freq = np.log(freq)
    return a * log_freq**4 + b * log_freq**3 + c * log_freq**2 + d * log_freq + e

start_time = time.time()
total = len(sources_xr.coords["index"])
indices = sources_xr.coords["index"]
if is_interactive():
    from tqdm.notebook import tqdm
    iterator = tqdm(indices, desc="Fitting sources")
else:
    iterator = indices

for i, s in enumerate(iterator):
    sources_xr["logpolycoefflux"].loc[dict(index=s)] = curve_fit(
        model, sources_xr.coords["freq"], sources_xr.flux.sel(index=s)
    )[0]
    sources_xr["logpolycoefpolflux"].loc[dict(index=s)] = curve_fit(
        model, sources_xr.coords["freq"], sources_xr.polarized_flux.sel(index=s)
    )[0]
    sources_xr["logpolycoefnorm"].loc[dict(index=s)] = curve_fit(
        model, sources_xr.coords["freq"], sources_xr.fluxnorm.sel(index=s)
    )[0]
    if (i + 1) % 100 == 0 or (i + 1) == total:
        elapsed = time.time() - start_time
        rate = (i + 1) / elapsed
        remaining = (total - (i + 1)) / rate if rate > 0 else float("inf")
        msg = f"ETA: {int(remaining // 3600)}h {int((remaining % 3600) // 60)}m"
        if is_interactive():
            iterator.set_postfix_str(msg)
        else:
            print(f"Processed {i+1}/{total} | {msg}")

# for s in range(len(all_indices)):
#     plt.figure()
#     sources_xr.flux.sel(index=s).plot(marker="o", linestyle="none")  # , xscale="log")
#     sources_xr.fluxnorm.sel(index=s).plot(
#         marker="o", linestyle="none"
#     )  # , xscale="log")

#     plt.loglog(
#         sources_xr.coords["freq"],
#         model(sources_xr.coords["freq"], *sources_xr.logpolycoefflux.sel(index=s)),
#     )
#     plt.loglog(
#         sources_xr.coords["freq"],
#         model(sources_xr.coords["freq"], *sources_xr.logpolycoefnorm.sel(index=s)),
#     )
#     plt.grid()
#     break

sources_xr.logpolycoefflux.min(), sources_xr.logpolycoefflux.max()

(<xarray.DataArray 'logpolycoefflux' ()> Size: 8B
 array(-17557.80083613),
 <xarray.DataArray 'logpolycoefflux' ()> Size: 8B
 array(23993.59700202))

# plt.figure(figsize=(12, 5))
# plt.subplot(121)
# sources_xr.logpolycoefflux.plot(vmax=50, vmin=-50)
# plt.subplot(122)
# sources_xr.logpolycoefnorm.plot(vmax=50, vmin=-50)pwd

# plt.figure(figsize=(15, 8))

# for power in range(5):
#     plt.subplot(231 + power)

#     np.fabs(sources_xr.logpolycoefflux.loc[dict(power=power)]).plot.hist(
#         bins=np.logspace(-0, 4, 20), density=False, lw=3, label="fluxes"
#     )

#     np.fabs(sources_xr.logpolycoefnorm.loc[dict(power=power)]).plot.hist(
#         bins=np.logspace(-0, 4, 20),
#         density=False,
#         histtype="step",
#         lw=2,
#         label="normalized fluxes",
#         linestyle="--",
#     )
#     plt.grid()
#     plt.title(f"Power {power}")
#     plt.legend()
#     plt.xscale("log")
#     plt.xlabel(None)

output_catalog = sources_xr[["logpolycoefflux","logpolycoefpolflux"]]

#no need to override this again
#output_catalog["index"] = all_indices

output_catalog.logpolycoefflux.attrs["units"] = "Jy"
output_catalog.logpolycoefpolflux.attrs["units"] = "Jy"

for coord in ["theta", "phi"]:
    output_catalog = output_catalog.assign_coords(**{coord:(("index"), cat[coord][high_flux_sources_mask].astype(np.float64))})

output_catalog

<xarray.Dataset> Size: 39MB
Dimensions:             (index: 372255, power: 5)
Coordinates:
  * index               (index) int64 3MB 27897657 242175600 ... 166742332
  * power               (power) int64 40B 4 3 2 1 0
    theta               (index) float64 3MB 2.655 2.659 2.596 ... 1.377 1.426
    phi                 (index) float64 3MB 4.177 3.977 4.17 ... 0.4638 0.5054
Data variables:
    logpolycoefflux     (index, power) float64 15MB 0.5214 -12.34 ... 0.08684
    logpolycoefpolflux  (index, power) float64 15MB -1.76 36.18 ... 0.02163

xarray.Dataset

• Dimensions:
  • index: 372255
  • power: 5
• Coordinates: (4)
  • index
    (index)
    int64
    27897657 242175600 ... 166742332

    array([ 27897657, 242175600, 270516228, ..., 246732122, 263347056, 166742332])

  • power
    (power)
    int64
    4 3 2 1 0

    array([4, 3, 2, 1, 0])

  • theta
    (index)
    float64
    2.655 2.659 2.596 ... 1.377 1.426

    array([2.65485716, 2.65869761, 2.59632301, ..., 1.38046789, 1.37715018,
           1.42648816])

  • phi
    (index)
    float64
    4.177 3.977 4.17 ... 0.4638 0.5054

    array([4.17657375, 3.97703004, 4.17026234, ..., 0.45735464, 0.46375346,
           0.50541824])

• Data variables: (2)
  • logpolycoefflux
    (index, power)
    float64
    0.5214 -12.34 ... -0.05652 0.08684

    units :

    Jy

    array([[ 5.21352400e-01, -1.23354024e+01,  1.13774189e+02,
            -4.93836646e+02,  8.72595399e+02],
           [ 2.61441524e+01, -5.80425233e+02,  4.80132667e+03,
            -1.75578008e+04,  2.39935970e+04],
           [ 2.66880868e-01, -5.23683404e+00,  3.70532139e+01,
            -1.09576844e+02,  1.22323272e+02],
           ...,
           [ 4.05660108e-05, -9.72696923e-04,  8.86015615e-03,
            -3.65758529e-02,  5.83135637e-02],
           [ 5.34233369e-05, -1.27064302e-03,  1.14441393e-02,
            -4.65193710e-02,  7.26237416e-02],
           [ 6.67844595e-05, -1.57728282e-03,  1.40750621e-02,
            -5.65244951e-02,  8.68435068e-02]])

  • logpolycoefpolflux
    (index, power)
    float64
    -1.76 36.18 ... -0.01785 0.02163

    units :

    Jy

    array([[-1.75996758e+00,  3.61777496e+01, -2.71679505e+02,
             8.77739069e+02, -1.01911362e+03],
           [-1.99578096e+00,  4.08465942e+01, -3.05044085e+02,
             9.78760643e+02, -1.12740632e+03],
           [ 2.11514845e-01, -4.23964325e+00,  3.08037399e+01,
            -9.56311446e+01,  1.07191630e+02],
           ...,
           [-1.54668160e-05,  3.13226756e-04, -2.31025331e-03,
             7.29200519e-03, -8.18540545e-03],
           [ 8.16796007e-06, -1.65100357e-04,  1.22179409e-03,
            -3.92277016e-03,  4.63592341e-03],
           [ 3.48534352e-05, -7.16857842e-04,  5.42127642e-03,
            -1.78541146e-02,  2.16335531e-02]])

• Indexes: (2)
  • index
    PandasIndex

    PandasIndex(Index([ 27897657, 242175600, 270516228, 276407774, 256213631, 281544643,
           263719520, 259533227, 251679008, 252717495,
           ...
           238190295, 145826702, 108046020, 232265270, 174768449, 256408684,
           236382968, 246732122, 263347056, 166742332],
          dtype='int64', name='index', length=372255))

  • power
    PandasIndex

    PandasIndex(Index([4, 3, 2, 1, 0], dtype='int64', name='power'))

• Attributes: (0)

output_filename = "data/websky_high_flux_catalog_1mJy.h5"

!ls -l $output_filename

ls: cannot access 'data/websky_high_flux_catalog_1mJy.h5': No such file or directory

output_catalog.coords["theta"].attrs["units"] = "rad"
output_catalog.coords["phi"].attrs["units"] = "rad"
output_catalog.coords["theta"].attrs["reference_frame"] = "Galactic"
output_catalog.coords["phi"].attrs["reference_frame"] = "Galactic"

output_catalog.attrs["description"] = (
    "Websky catalog of sources with flux > 1 mJy at 100 GHz, fitted with a 4th order polynomial in log frequency. "
    "Galactic reference frame. Sorted by flux at 100 GHz (descending). "
    "The 'index' coordinate gives the original index in the Websky catalog."
)

output_catalog.to_netcdf(output_filename, format="NETCDF4", mode="w")  # requires netcdf4 package

%ls -lah $output_filename

-rw-rw-r-- 1 azonca azonca 37M May 12 22:41 data/websky_high_flux_catalog_1mJy.h5

import xarray

xarray.open_dataset(output_filename)

<xarray.Dataset> Size: 39MB
Dimensions:             (index: 372255, power: 5)
Coordinates:
  * index               (index) int64 3MB 27897657 242175600 ... 166742332
  * power               (power) int64 40B 4 3 2 1 0
    theta               (index) float64 3MB ...
    phi                 (index) float64 3MB ...
Data variables:
    logpolycoefflux     (index, power) float64 15MB ...
    logpolycoefpolflux  (index, power) float64 15MB ...
Attributes:
    description:  Websky catalog of sources with flux > 1 mJy at 100 GHz, fit...

xarray.Dataset

• Dimensions:
  • index: 372255
  • power: 5
• Coordinates: (4)
  • index
    (index)
    int64
    27897657 242175600 ... 166742332

    array([ 27897657, 242175600, 270516228, ..., 246732122, 263347056, 166742332])

  • power
    (power)
    int64
    4 3 2 1 0

    array([4, 3, 2, 1, 0])

  • theta
    (index)
    float64
    ...

    units :

    rad

    reference_frame :

    Galactic

    [372255 values with dtype=float64]

  • phi
    (index)
    float64
    ...

    units :

    rad

    reference_frame :

    Galactic

    [372255 values with dtype=float64]

• Data variables: (2)
  • logpolycoefflux
    (index, power)
    float64
    ...

    units :

    Jy

    [1861275 values with dtype=float64]

  • logpolycoefpolflux
    (index, power)
    float64
    ...

    units :

    Jy

    [1861275 values with dtype=float64]

• Indexes: (2)
  • index
    PandasIndex

    PandasIndex(Index([ 27897657, 242175600, 270516228, 276407774, 256213631, 281544643,
           263719520, 259533227, 251679008, 252717495,
           ...
           238190295, 145826702, 108046020, 232265270, 174768449, 256408684,
           236382968, 246732122, 263347056, 166742332],
          dtype='int64', name='index', length=372255))

  • power
    PandasIndex

    PandasIndex(Index([4, 3, 2, 1, 0], dtype='int64', name='power'))

• Attributes: (1)

  description :

  Websky catalog of sources with flux > 1 mJy at 100 GHz, fitted with a 4th order polynomial in log frequency. Galactic reference frame. Sorted by flux at 100 GHz (descending). The 'index' coordinate gives the original index in the Websky catalog.

import h5py
f = h5py.File(output_filename, 'r')
f["logpolycoefflux"]

<HDF5 dataset "logpolycoefflux": shape (372255, 5), type "<f8">

f["logpolycoefflux"].attrs["units"]

b'Jy'