WebSky bright source catalog workflow

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

/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"

import os

# for jupyter.nersc.gov otherwise the notebook only uses 2 cores
num_threads = 128
os.environ["OMP_NUM_THREADS"] = "1"

cutoff_flux = 1e-3

output_filename = "websky_high_flux_catalog_1mJy.h5"

plot = False

cd /global/cfs/cdirs/sobs/www/users/Radio_WebSky/matched_catalogs_2

/global/cfs/cdirs/sobs/www/users/Radio_WebSky/matched_catalogs_2

%ls

catalog_100.0.h5  catalog_232.0.h5  catalog_353.0.h5  catalog_643.0.h5
catalog_111.0.h5  catalog_24.5.h5   catalog_375.0.h5  catalog_67.8.h5
catalog_129.0.h5  catalog_256.0.h5  catalog_409.0.h5  catalog_70.0.h5
catalog_143.0.h5  catalog_27.3.h5   catalog_41.7.h5   catalog_729.0.h5
catalog_153.0.h5  catalog_275.0.h5  catalog_44.0.h5   catalog_73.7.h5
catalog_164.0.h5  catalog_294.0.h5  catalog_467.0.h5  catalog_79.6.h5
catalog_18.7.h5   catalog_30.0.h5   catalog_47.4.h5   catalog_817.0.h5
catalog_189.0.h5  catalog_306.0.h5  catalog_525.0.h5  catalog_857.0.h5
catalog_21.6.h5   catalog_314.0.h5  catalog_545.0.h5  catalog_90.2.h5
catalog_210.0.h5  catalog_340.0.h5  catalog_584.0.h5  catalog_906.0.h5
catalog_217.0.h5  catalog_35.9.h5   catalog_63.9.h5

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("catalog_100.0.h5", "r")

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

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

cat["polarized flux"][:3]

array([1.42910628e-09, 1.99535624e-08, 2.29563857e-09], dtype='>f8')

(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 = np.array(all_indices)

Generate 1 source only

argmax = np.array(cat["flux"]).argmax()
all_indices = np.array([argmax])
high_flux_sources_mask = high_flux_sources_mask * 0
high_flux_sources_mask[argmax] = 1

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 tqdm(freqs):
    cat = h5py.File(f"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]

100%|██████████| 13/13 [00:20<00:00,  1.54s/it]

# sources_xr = sources_xr.sortby(sources_xr.flux.loc[dict(freq=float(freqs[0]))])

sources_xr.coords["index"] = np.arange(len(sources_xr.coords["index"]))

if plot:
    for s in tqdm(range(len(all_indices))):
        sources_xr["fluxnorm"].loc[dict(index=s)] = sources_xr["flux"].loc[
            dict(index=s)
        ] / sources_xr["flux"].loc[dict(index=s)].sel(freq=100)

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

if plot:
    #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["logpolycoefpolflux"] = sources_xr["logpolycoefflux"].copy()

if plot:
    sources_xr["logpolycoefnorm"] = sources_xr["logpolycoefflux"].copy()

from numba import njit

@njit
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

from scipy.optimize import curve_fit

from dask.distributed import Client, LocalCluster

cluster = LocalCluster(n_workers=32, threads_per_worker=2, processes=True)
client = Client(cluster)

print(client)

<Client: 'tcp://127.0.0.1:41433' processes=32 threads=64, memory=503.14 GiB>

s=0
sources_xr["flux"].sel(index=s).data

array([0.02853885, 0.01962042, 0.0104193 , 0.00548068, 0.00337772,
       0.00218903, 0.0014772 , 0.00095432, 0.00064614, 0.000557  ,
       0.00049761, 0.00043035, 0.00040939])

def run_curve_fit_factory(field):
    def run_curve_fit(s):
        return curve_fit(
            model, sources_xr.coords["freq"].data, sources_xr[field].sel(index=s).data
        )[0]
    return run_curve_fit

from dask.diagnostics import ProgressBar

import dask.bag as db

from dask.distributed import progress

%%time

sources_xr["logpolycoefflux"].data = db.range(len(all_indices), npartitions=num_threads).map(run_curve_fit_factory("flux")).compute()

/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/distributed/client.py:3361: UserWarning: Sending large graph of size 142.02 MiB.
This may cause some slowdown.
Consider loading the data with Dask directly
 or using futures or delayed objects to embed the data into the graph without repetition.
See also https://docs.dask.org/en/stable/best-practices.html#load-data-with-dask for more information.
  warnings.warn(
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"
/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/scipy/__init__.py:155: UserWarning: A NumPy version >=1.18.5 and <1.26.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"

CPU times: user 18.7 s, sys: 11.5 s, total: 30.2 s
Wall time: 39 s

%%time

sources_xr["logpolycoefpolflux"].data = db.range(len(all_indices), npartitions=num_threads).map(run_curve_fit_factory("polarized_flux")).compute()

/global/common/software/cmb/zonca/conda/pycmb/lib/python3.10/site-packages/distributed/client.py:3361: UserWarning: Sending large graph of size 142.02 MiB.
This may cause some slowdown.
Consider loading the data with Dask directly
 or using futures or delayed objects to embed the data into the graph without repetition.
See also https://docs.dask.org/en/stable/best-practices.html#load-data-with-dask for more information.
  warnings.warn(

CPU times: user 16.2 s, sys: 10.1 s, total: 26.3 s
Wall time: 33.3 s

%%time

if plot:
    sources_xr["logpolycoefnorm"].data = db.range(len(all_indices), npartitions=num_threads).map(run_curve_fit_factory("fluxnorm")).compute()

CPU times: user 3 µs, sys: 3 µs, total: 6 µs
Wall time: 11 µs

# for s in range(len(all_indices)):
#     sources_xr["logpolycoefflux"].loc[dict(index=s)], cov = curve_fit(
#         model, sources_xr.coords["freq"], sources_xr.flux.sel(index=s)
#     )
#     sources_xr["logpolycoefpolflux"].loc[dict(index=s)], cov = curve_fit(
#         model, sources_xr.coords["freq"], sources_xr.polarized_flux.sel(index=s)
#     )
#     if plot:
#         sources_xr["logpolycoefnorm"].loc[dict(index=s)], cov = curve_fit(
#             model, sources_xr.coords["freq"], sources_xr.fluxnorm.sel(index=s)
#         )

if plot:
    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()

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

(<xarray.DataArray 'logpolycoefflux' ()>
 array(-17557.7742453),
 <xarray.DataArray 'logpolycoefflux' ()>
 array(23993.56673443))

if plot:
    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)

if plot:
    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"]]

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>
Dimensions:             (index: 372255, power: 5)
Coordinates:
  * index               (index) int64 11253 16428 24110 ... 281755430 281755795
  * power               (power) int64 4 3 2 1 0
    theta               (index) float64 2.655 2.659 2.596 ... 1.38 1.377 1.426
    phi                 (index) float64 4.177 3.977 4.17 ... 0.4638 0.5054
Data variables:
    logpolycoefflux     (index, power) float64 0.0002901 -0.006796 ... 0.01476
    logpolycoefpolflux  (index, power) float64 5.318e-05 -0.001129 ... -0.005808

xarray.Dataset

• Dimensions:
  • index: 372255
  • power: 5
• Coordinates: (4)
  • index
    (index)
    int64
    11253 16428 ... 281755430 281755795

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

  • 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.0002901 -0.006796 ... 0.01476

    units :

    Jy

    array([[ 2.90053416e-04, -6.79555871e-03,  6.00239441e-02,
            -2.37922456e-01,  3.59419360e-01],
           [ 1.31256483e-03, -2.96609404e-02,  2.50643946e-01,
            -9.40591152e-01,  1.32678333e+00],
           [ 4.98663676e-05, -1.20799507e-03,  1.12272497e-02,
            -4.79160078e-02,  8.03926627e-02],
           ...,
           [ 2.88200893e-05, -5.58959572e-04,  3.89880397e-03,
            -1.14879891e-02,  1.33987434e-02],
           [ 2.67619262e-05, -5.14763210e-04,  3.57754498e-03,
            -1.07709891e-02,  1.37275093e-02],
           [ 3.22674612e-05, -6.31555794e-04,  4.45254848e-03,
            -1.31392330e-02,  1.47624425e-02]])

  • logpolycoefpolflux
    (index, power)
    float64
    5.318e-05 -0.001129 ... -0.005808

    units :

    Jy

    array([[ 5.31765869e-05, -1.12884428e-03,  8.85423431e-03,
            -3.04050832e-02,  3.86182943e-02],
           [ 7.11030208e-04, -1.47609159e-02,  1.12776947e-01,
            -3.75177643e-01,  4.57867767e-01],
           [ 3.53213070e-06, -7.17390506e-05,  5.28384844e-04,
            -1.67717417e-03,  1.98338291e-03],
           ...,
           [-2.56380115e-05,  5.27290371e-04, -3.96098587e-03,
             1.28171900e-02, -1.49508938e-02],
           [ 6.19757587e-06, -1.10323374e-04,  6.88961558e-04,
            -1.75523928e-03,  1.54990839e-03],
           [-1.00676079e-05,  2.04517752e-04, -1.52694624e-03,
             4.95085915e-03, -5.80847642e-03]])

• 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))

  • power
    PandasIndex

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

• Attributes: (0)

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

flux_100 = np.polynomial.polynomial.polyval(np.log(100), np.array(output_catalog["logpolycoefflux"])[:,::-1].T, tensor=False)

output_catalog["flux_100"] = np.polynomial.polynomial.polyval(np.log(100), output_catalog["logpolycoefflux"][:,::-1].T, tensor=False)

output_catalog["flux_100"].max()

<xarray.DataArray 'flux_100' ()>
array(41.02290808)

xarray.DataArray

'flux_100'

• 41.02

  array(41.02290808)

• Coordinates: (0)
• Indexes: (0)
• Attributes: (0)

np.array(output_catalog["logpolycoefflux"]).shape

(372255, 5)

output_catalog = output_catalog.sortby("flux_100", ascending=False)

del output_catalog["flux_100"]

np.polyval(output_catalog.logpolycoefflux[0], np.log(100))

0.003478966281996676

output_catalog.attrs["notes"] = \
"""Catalog of sources where the flux in Jy at any frequency is calculated with a 5th order polynomial in the logarithm of the frequency in GHz, separately for temperature and polarization.
The catalog does not contain information about the polarization angle of a source.
The catalog sorted in descending order based on the source flux at 100 GHz"""

output_catalog

<xarray.Dataset>
Dimensions:             (index: 372255, power: 5)
Coordinates:
  * index               (index) int64 11253 16428 24110 ... 281755430 281755795
  * power               (power) int64 4 3 2 1 0
    theta               (index) float64 2.655 2.659 2.596 ... 1.38 1.377 1.426
    phi                 (index) float64 4.177 3.977 4.17 ... 0.4638 0.5054
Data variables:
    logpolycoefflux     (index, power) float64 0.0002901 -0.006796 ... 0.01476
    logpolycoefpolflux  (index, power) float64 5.318e-05 -0.001129 ... -0.005808
Attributes:
    notes:    Catalog of sources where the flux in Jy at any frequency is cal...

xarray.Dataset

• Dimensions:
  • index: 372255
  • power: 5
• Coordinates: (4)
  • index
    (index)
    int64
    11253 16428 ... 281755430 281755795

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

  • 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

    units :

    rad

    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

    units :

    rad

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

• Data variables: (2)
  • logpolycoefflux
    (index, power)
    float64
    0.0002901 -0.006796 ... 0.01476

    units :

    Jy

    array([[ 2.90053416e-04, -6.79555871e-03,  6.00239441e-02,
            -2.37922456e-01,  3.59419360e-01],
           [ 1.31256483e-03, -2.96609404e-02,  2.50643946e-01,
            -9.40591152e-01,  1.32678333e+00],
           [ 4.98663676e-05, -1.20799507e-03,  1.12272497e-02,
            -4.79160078e-02,  8.03926627e-02],
           ...,
           [ 2.88200893e-05, -5.58959572e-04,  3.89880397e-03,
            -1.14879891e-02,  1.33987434e-02],
           [ 2.67619262e-05, -5.14763210e-04,  3.57754498e-03,
            -1.07709891e-02,  1.37275093e-02],
           [ 3.22674612e-05, -6.31555794e-04,  4.45254848e-03,
            -1.31392330e-02,  1.47624425e-02]])

  • logpolycoefpolflux
    (index, power)
    float64
    5.318e-05 -0.001129 ... -0.005808

    units :

    Jy

    array([[ 5.31765869e-05, -1.12884428e-03,  8.85423431e-03,
            -3.04050832e-02,  3.86182943e-02],
           [ 7.11030208e-04, -1.47609159e-02,  1.12776947e-01,
            -3.75177643e-01,  4.57867767e-01],
           [ 3.53213070e-06, -7.17390506e-05,  5.28384844e-04,
            -1.67717417e-03,  1.98338291e-03],
           ...,
           [-2.56380115e-05,  5.27290371e-04, -3.96098587e-03,
             1.28171900e-02, -1.49508938e-02],
           [ 6.19757587e-06, -1.10323374e-04,  6.88961558e-04,
            -1.75523928e-03,  1.54990839e-03],
           [-1.00676079e-05,  2.04517752e-04, -1.52694624e-03,
             4.95085915e-03, -5.80847642e-03]])

• 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))

  • power
    PandasIndex

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

• Attributes: (1)

  notes :

  Catalog of sources where the flux in Jy at any frequency is calculated with a 5th order polynomial in the logarithm of the frequency in GHz, separately for temperature and polarization. The catalog does not contain information about the polarization angle of a source. The catalog sorted in descending order based on the source flux at 100 GHz

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

%ls -lah $output_filename

-rw-rw---- 1 zonca sobs 37M Sep 19 11:15 websky_high_flux_catalog_1mJy.h5

import xarray

xarray.open_dataset(output_filename)

<xarray.Dataset>
Dimensions:             (index: 372255, power: 5)
Coordinates:
  * index               (index) int64 11253 16428 24110 ... 281755430 281755795
  * power               (power) int64 4 3 2 1 0
    theta               (index) float64 ...
    phi                 (index) float64 ...
Data variables:
    logpolycoefflux     (index, power) float64 ...
    logpolycoefpolflux  (index, power) float64 ...
Attributes:
    notes:    Catalog of sources where the flux in Jy at any frequency is cal...

xarray.Dataset

• Dimensions:
  • index: 372255
  • power: 5
• Coordinates: (4)
  • index
    (index)
    int64
    11253 16428 ... 281755430 281755795

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

  • power
    (power)
    int64
    4 3 2 1 0

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

  • theta
    (index)
    float64
    ...

    units :

    rad

    [372255 values with dtype=float64]

  • phi
    (index)
    float64
    ...

    units :

    rad

    [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([    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))

  • power
    PandasIndex

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

• Attributes: (1)

  notes :

  Catalog of sources where the flux in Jy at any frequency is calculated with a 5th order polynomial in the logarithm of the frequency in GHz, separately for temperature and polarization. The catalog does not contain information about the polarization angle of a source. The catalog sorted in descending order based on the source flux at 100 GHz

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

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

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

b'Jy'

!mv $output_filename ~/p/issues/202405_pysm_catalog_pixell/