Add mask folder with point source finding and point source mask scripts

Author: merryduparc

project/SPT_dev/mask/make_source_mask.py

@@ -0,0 +1,39 @@
+import numpy as np
+from matplotlib import pyplot as plt
+from pixell import enmap, enplot, reproject
+from pspy import so_map
+import healpy as hp
+from os.path import join as opj
+
+catalog = np.loadtxt("maps/catalogs/cat.txt").T
+
+save_dir = ""
+spt_mask_dir = opj(save_dir, "masks")
+
+header = ['ra', 'dec', 'SNR', 'Tamp', 'dTamp', 'Qamp', 'dQamp', 'Uamp', 'dUamp', 'Tflux', 'dTflux', 'Qflux', 'dQflux', 'Uflux', 'dUflux', 'npix', 'status']
+catalog_dict = {k:v for k, v in zip(header, catalog)}
+
+car_template = so_map.from_enmap(enmap.read_map(opj(spt_mask_dir, 'pixel_mask_binary_borders_only_CAR.fits')))
+
+SNR_cut = 5
+
+for flux_cut in [100, 50, 15]:
+    flux_mask = (catalog_dict["SNR"] >  5) & (catalog_dict["Tflux"] > flux_cut)
+
+    ps_mask = so_map.generate_source_mask(
+        car_template, 
+        coordinates=np.deg2rad([np.array(catalog_dict['dec'][flux_mask]), np.array(catalog_dict['ra'][flux_mask])]), 
+        point_source_radius_arcmin=12
+    ).data
+    enmap.write_map(opj(spt_mask_dir, f"point_source_test_{flux_cut}_CAR.fits"), ps_mask)
+    plot = enplot.get_plots(
+        ps_mask, ticks=10, mask=0, colorbar=True, downgrade=2, range=(1)
+    )
+    enplot.write(opj(spt_mask_dir, f"point_source_test_{flux_cut}_CAR"), plot)
+    
+    ps_mask_healpix = reproject.map2healpix(ps_mask, nside=8096, method="spline")
+    hp.mollview(ps_mask_healpix)
+    plt.savefig(opj(spt_mask_dir, f"point_source_test_{flux_cut}.png"))
+    hp.write_map(opj(spt_mask_dir, f"point_source_test_{flux_cut}.fits"), ps_mask_healpix)
+
+

project/SPT_dev/mask/point_sources.rst

@@ -0,0 +1,37 @@
+**************************
+Finding sources in SPT maps
+**************************
+
+We use ``dory`` from tenki (https://github.com/amaurea/tenki/tree/master) to find sources in SPT maps. 
+We first need to project SPT masks and maps to CAR so dory can read these. 
+``project_maps_SPT_patch.py`` defines a CAR template around SPT survey and projects 90GHz full maps, ivar (normalized ivar * 1e-2 to mimic a 5uK.arcmin depth, but we should look for raw ivar maps) and masks.
+Go in ``python/mask/project_maps_SPT_patch.py`` and modify ``save_dir`` if you want to save maps somewhere in particular, then run:
+
+.. code:: shell
+    salloc --nodes 1 --qos interactive --time 00:20:00 --constraint cpu
+    OMP_NUM_THREADS=256 srun -n 1 -c 256 --cpu-bind=cores python mask/project_maps_SPT_patch.py
+
+You can now get dory by cloning ``tenki``.
+Before running ``dory``, I had some problems with floating point precision (?) in ``tenki/point_sources/dory.py``, so go to line 119 and change ``> 0`` to ``> 1e-10``.
+Otherwise, masks will look very weird and the code will not find any sources.
+We can now run ``dory`` (change your path to dory if needed), you can remove "debug" in --output if you trust what I did (I wouldn't):
+.. code:: shell
+    salloc --nodes 1 --qos interactive --time 01:00:00 --constraint cpu
+    OMP_NUM_THREADS=64 srun -n 4 -c 64 --cpu-bind=cores python tenki/point_sources/dory.py find \
+        maps/spt/full_095ghz_CAR.fits \
+        maps/spt/ivar_095ghz_CAR.fits \
+        maps/catalogs \
+        --freq 90 \
+        --beam /global/cfs/cdirs/sobs/users/tlouis/spt_data/ancillary_products/generally_applicable/beam_c26_cmb_bl_095ghz.txt \
+        -R "tile:3000" \
+        --pad 60 \
+        -m maps/full/pixel_mask_binary_borders_only_CAR_REV.fits \
+        --nsigma 3 \
+        --output "full,reg,debug" \
+        --verbose
+
+You now have a catalog of sources at ``maps/catalogs/cat.{txt,fits}``.
+A first try to make masks with it is in ``make_source_mask.py``:
+.. code:: shell
+    salloc --nodes 1 --qos interactive --time 00:05:00 --constraint cpu
+    OMP_NUM_THREADS=256 srun -n 1 -c 256 --cpu-bind=cores python mask/make_source_mask.py

project/SPT_dev/mask/project_maps_SPT_patch.py

@@ -0,0 +1,93 @@
+"""
+This script projects SPT and ACT maps and masks into a CAR template defined around SPT survey
+"""
+
+from pixell import enplot, enmap, reproject, utils
+import numpy as np
+from pspy import so_map
+from os.path import join as opj
+import os
+
+save_dir = ""
+spt_maps_dir = opj(save_dir, "maps/spt")
+spt_mask_dir = opj(save_dir, "masks")
+
+os.makedirs(spt_maps_dir, exist_ok=True)
+os.makedirs(spt_mask_dir, exist_ok=True)
+
+### Define project geometry (arbitrary box around SPT patch + few degrees)
+
+print("create SPT CAR geometry")
+box = [[-77, -62],[-35,60]] * utils.adeg
+shape,wcs = enmap.geometry2(pos=box,res=0.5 * utils.arcmin,proj='car', variant='fejer1')
+print(f"Project on geometry: {shape=}, {wcs=}")
+template = enmap.ones(shape, wcs)
+enmap.write_map(opj(save_dir, "SPT_CAR_template.fits"), template)
+
+
+### Project SPT full maps
+
+spt_maps_read_dir = "/global/cfs/cdirs/sobs/users/tlouis/spt_data/real_data_maps"
+spt_full_maps_fn_list = [opj(spt_maps_read_dir, f"full/full_{freq}ghz.fits") for freq in ["095"]]
+
+for spt_maps_fn in spt_full_maps_fn_list: 
+    # Project maps
+    maps_so = so_map.read_map(spt_maps_fn, fields_healpix=(0, 1, 2))
+    maps = maps_so.data
+
+    maps_proj = reproject.healpix2map(maps, shape, wcs, method='harm')    # For maps we use harm
+    
+    enmap.write_map(opj(spt_maps_dir, os.path.basename(spt_maps_fn)[:-5] + "_CAR.fits"), maps_proj)
+
+    plot = enplot.get_plots(
+        maps_proj, ticks=10, mask=0, downgrade=2, colorbar=True, range=(100, 30, 30)
+    )
+    enplot.write(opj(spt_maps_dir, os.path.basename(spt_maps_fn)[:-5] + "_CAR"), plot)
+    
+    # We also need ivar, we use a fudge factor to get close to the real ivar: here I use 1e-2 which corresponds to 5uk.arcmin depth
+    fudge_ivar = 1e-2
+    maps_so = so_map.read_map(spt_maps_fn, fields_healpix=(3))
+    maps = maps_so.data * fudge_ivar
+
+    maps_proj = reproject.healpix2map(maps, shape, wcs, method='harm')    # For maps we use harm
+    
+    enmap.write_map(opj(spt_maps_dir, os.path.basename(spt_maps_fn)[:-5].replace("full", "ivar") + "_CAR.fits"), maps_proj)
+
+    plot = enplot.get_plots(
+        maps_proj, ticks=10, mask=0, downgrade=2, colorbar=True, range=(1)
+    )
+    enplot.write(opj(spt_maps_dir, os.path.basename(spt_maps_fn)[:-5].replace("full", "ivar") + "_CAR"), plot)    
+
+
+spt_masks_read_dir = "/global/cfs/cdirs/sobs/users/tlouis/spt_data/ancillary_products/generally_applicable"
+spt_masks_fn_list = [
+    opj(spt_masks_read_dir, "pixel_mask_apodized_borders_objects.fits"),
+    opj(spt_masks_read_dir, "pixel_mask_apodized_borders_only.fits"),
+    opj(spt_masks_read_dir, "pixel_mask_binary_borders_objects.fits"),
+    opj(spt_masks_read_dir, "pixel_mask_binary_borders_only.fits"),
+]
+
+for spt_mask_fn in spt_masks_fn_list: 
+    # Project masks
+    maps_so = so_map.read_map(spt_mask_fn)
+    maps = maps_so.data
+
+    maps_proj = reproject.healpix2map(maps, shape, wcs, method='spline')    # For masks we use spline
+    
+    enmap.write_map(opj(spt_mask_dir, os.path.basename(spt_mask_fn)[:-5] + "_CAR.fits"), maps_proj)
+
+    plot = enplot.get_plots(
+        maps_proj, ticks=10, mask=0, downgrade=2, colorbar=True, range=1
+    )
+    enplot.write(opj(spt_mask_dir, os.path.basename(spt_mask_fn)[:-5] + "_CAR"), plot)
+    
+    # We also need reverted masks for dory
+    enmap.write_map(opj(spt_mask_dir, os.path.basename(spt_mask_fn)[:-5] + "_CAR_rev.fits"), 1 - maps_proj)
+
+    plot = enplot.get_plots(
+        1 - maps_proj, ticks=10, mask=0, downgrade=2, colorbar=True, range=1
+    )
+    enplot.write(opj(spt_mask_dir, os.path.basename(spt_mask_fn)[:-5] + "_CAR_rev"), plot)
+
+
+