diff --git a/pyradiosky/data/cat_mock.dat b/pyradiosky/data/cat_mock.dat new file mode 100644 index 00000000..6193a674 --- /dev/null +++ b/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 diff --git a/pyradiosky/tests/test_utils.py b/pyradiosky/tests/test_utils.py index 89f23661..55d02f2c 100644 --- a/pyradiosky/tests/test_utils.py +++ b/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 diff --git a/pyradiosky/utils.py b/pyradiosky/utils.py index 1d235b72..e04496c3 100644 --- a/pyradiosky/utils.py +++ b/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