Modified gleam

pyradiosky/data/cat_mock.dat

@@ -0,0 +1,2 @@
+2.013e+02 -4.301e+01 1.937e+03 -5.000e-01
+1.395e+02 -1.209e+01 5.446e+02 -9.599e-01

pyradiosky/tests/test_utils.py

@@ -134,3 +134,38 @@ def jy2ksr_nonastropy(freq_arr):
     conv1 = jy2ksr_nonastropy(freqs) * units.K * units.sr / units.Jy
 
     assert np.allclose(conv0, conv1)
+
+
+def test_modified_gleam():
+
+    pytest.importorskip("astropy_healpix")
+
+    # if file is not present in that directory
+    with pytest.raises(ValueError, match="File not found"):
+        skyutils.modified_gleam("testgleam.dat")
+
+    # read the dummy cat file present in the directory
+    cat1 = skyutils.modified_gleam(
+        "pyradiosky/data/cat_mock.dat",
+        usecols=(0, 1, 2, 3),
+        modified_gleam_filename="pyradiosky/data/cat_moc_out.dat",
+    )
+
+    # read the dummy cat file and fill the blank regions
+    cat2 = skyutils.modified_gleam(
+        "pyradiosky/data/cat_mock.dat", usecols=(0, 1, 2, 3), fill_blank=True, nside=8
+    )
+
+    # read the dummy cat file and fill the blank regions,and also add the peeled sources
+    cat3 = skyutils.modified_gleam(
+        "pyradiosky/data/cat_mock.dat",
+        usecols=(0, 1, 2, 3),
+        fill_blank=True,
+        nside=8,
+        add_peeled_sources=True,
+    )
+
+    assert len(cat2) > len(cat1)
+
+    # 9 is added peeled source
+    assert len(cat3) - len(cat2) == 9

pyradiosky/utils.py

@@ -4,6 +4,7 @@
 """Utility methods."""
 import os
 import warnings
+import copy
 
 import numpy as np
 
@@ -19,12 +20,13 @@
 import astropy.units as units
 from astropy.units import Quantity
 
-
 # The frame radio astronomers call the apparent or current epoch is the
 # "true equator & equinox" frame, notated E_upsilon in the USNO circular
 # astropy doesn't have this frame but it's pretty easy to adapt the CIRS frame
 # by modifying the ra to reflect the difference between
 # GAST (Grenwich Apparent Sidereal Time) and the earth rotation angle (theta)
+
+
 def _tee_to_cirs_ra(tee_ra, time):
     """
     Convert from the true equator & equinox frame to the CIRS frame.
@@ -278,3 +280,113 @@ def download_gleam(path=".", filename="gleam.vot", overwrite=False, row_limit=No
     table.write(opath, format="votable", overwrite=overwrite)
 
     print("GLEAM catalog downloaded and saved to " + opath)
+
+
+def modified_gleam(
+    gleam_filename="gleamegc.dat",
+    usecols=(10, 12, 77, -5),
+    fill_blank=False,
+    nside=32,
+    add_peeled_sources=False,
+    modified_gleam_filename="",
+):
+    """
+    Write modified gleam catalogue to fill the blank regions, and also to add the peeled sources.
+
+    Parameters
+    ----------
+    gleam_filename : str
+        name of the gleam catalogue file.
+    usecols : tuple
+        default = (10,12,77,-5)
+        4 columns, 1st-RA (deg), 2nd-Dec (deg), 3rd - Flux at 100 MHz (Jy), 4th - sp. index
+    fill_blank : bool
+        If True, it will fill the blank regions
+    nside : int
+        Use for healpix pixalization
+    add_peeled_sources : bool
+        If True, it will add the peeled sources from Table 2 Hurley-Walker, N et al. 2017
+    modified_gleam_filename : str
+        name of the modified gleam catalogue file.
+
+    """
+    import astropy_healpix as astrohp
+    import astropy_healpix.healpy as hp
+
+    # check if file exists in the given path
+    if os.path.exists(gleam_filename):
+        print("true gleam catalogue file - ", gleam_filename)
+    else:
+        raise ValueError("File not found")
+
+    # read array of RA (deg), Dec (deg), Flux at 100 MHz (Jy), sp. index
+    aa = np.genfromtxt(gleam_filename, usecols=usecols)
+
+    # use sources with finite flux and sp. index values
+    cat_gleam = aa[(aa[:, 2] >= 0.0) & np.isfinite(aa[:, 2]) & np.isfinite(aa[:, 3])]
+
+    # fill blank regions
+    if fill_blank:
+
+        npix = hp.nside2npix(nside)
+        numhpxmap = np.zeros(npix, dtype=np.float)
+        indices = hp.ang2pix(nside, cat_gleam[:, 0], cat_gleam[:, 1], lonlat=True)
+
+        # fill hpmap with number of source
+        for i in range(npix):
+            wr = np.where(indices == i)
+            numhpxmap[i] = wr[0].size
+
+        # non zero indices in the hpmap
+        non_zero_ind = np.nonzero(numhpxmap)[0]
+
+        # find zero indices with neighbour >2
+        nn = astrohp.neighbours(np.arange(npix), nside)
+        zero_ind = []
+
+        for i in range(npix):
+            if np.where(numhpxmap[nn[:, i]] == 0)[0].size > 2:
+                zero_ind.append(i)
+        zero_ind = np.array(zero_ind)
+
+        # fill zero indices with randomly chosen non zero indices
+        np.random.seed(10)
+        zero_fill_ind = np.random.choice(non_zero_ind, size=len(zero_ind))
+
+        numhpxmap_fill = copy.deepcopy(numhpxmap)
+        numhpxmap_fill[zero_ind] = numhpxmap[zero_fill_ind]
+
+        ra_zero, dec_zero = hp.pix2ang(nside, zero_ind, lonlat=True)
+
+        cat_random = []
+        for i in range(len(ra_zero)):
+            wr = np.where(indices == zero_fill_ind[i])
+            for j in wr[0]:
+                cat_random.append(
+                    (ra_zero[i], dec_zero[i], cat_gleam[j, 2], cat_gleam[j, 3])
+                )
+        cat_random = np.array(cat_random)
+        cat_gleam = np.vstack((cat_gleam, cat_random))
+
+    # add peeled sources
+    if add_peeled_sources:
+        gleam_peeled = np.array(
+            [
+                [201.3667, -43.0192, 1937.472, -0.50],  # Centaurus A
+                [139.5250, -12.0956, 544.686, -0.96],  # Hydra A
+                [79.9583, -45.7789, 774.612, -0.99],  # Pictor A
+                [252.7833, 4.9925, 791.486, -1.07],  # Hercules A
+                [187.7042, 12.3911, 1562.747, -0.86],  # Virgo A
+                [83.6333, 22.0144, 1560.743, -0.22],  # Crab
+                [299.8667, 40.7339, 13599.676, -0.78],  # Cygnus A
+                [350.8667, 58.8117, 15811.361, -0.41],  # Cassiopeia A
+                [50.6708, -37.2083, 1256.016, -0.81],
+            ]
+        )  # Fornax A
+        cat_gleam = np.vstack((gleam_peeled, cat_gleam))
+
+    # save in a file
+    if modified_gleam_filename != "":
+        np.savetxt(modified_gleam_filename, cat_gleam)
+
+    return cat_gleam