aboutsummaryrefslogtreecommitdiffstats
path: root/bin/get-healpix-patch
blob: 0dc27f6bcc9349c5378f0e50e94271d8ca632fd2 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
#!/usr/bin/env python3
#
# Copyright (c) 2017 Weitian LI <weitian@aaronly.me>
# MIT license
#

"""
Extract a patch of sky from the all-sky HEALPix map.
"""


import os
import sys
import argparse
import logging

import scipy.ndimage
import healpy as hp
from reproject import reproject_from_healpix

from fg21sim.configs import configs
from fg21sim.sky import SkyPatch
from fg21sim.utils import setup_logging
from fg21sim.utils.fits import read_fits_healpix
from fg21sim.utils.units import UnitConversions as AUC


def main():
    parser = argparse.ArgumentParser(
        description="Extract a patch from the all-sky HEALPix map")
    parser.add_argument("-C", "--clobber", action="store_true",
                        help="overwrite the existing output files")
    parser.add_argument("-c", "--config", dest="config", required=False,
                        help="fg21sim configuration from which to " +
                        "obtain the sky patch properties")
    parser.add_argument("--center", dest="center",
                        help="center coordinate of the sky patch; " +
                        "format: ra,dec; unit: deg")
    parser.add_argument("--size", dest="size",
                        help="size of the sky patch; " +
                        "format: xsize,ysize; unit: pixel")
    parser.add_argument("--pixelsize", dest="pixelsize", type=float,
                        help="pixel size of the sky patch; unit: [arcsec]")
    parser.add_argument("-S", "--smooth", dest="smooth", action="store_true",
                        help="Smooth the output patch using a Gaussian " +
                        "filter whose sigma is 2x the pixel size " +
                        "of the input HEALPix map")
    parser.add_argument("infile", help="input all-sky HEALPix map")
    parser.add_argument("outfile", help="output extracted sky patch")
    args = parser.parse_args()

    setup_logging(dict_config=configs.logging)
    tool = os.path.basename(sys.argv[0])
    logger = logging.getLogger(tool)
    logger.info("COMMAND: {0}".format(" ".join(sys.argv)))

    if args.config:
        configs.read_userconfig(args.config)
        center = (configs.getn("sky/patch/xcenter"),
                  configs.getn("sky/patch/ycenter"))  # [ deg ]
        size = (configs.getn("sky/patch/xsize"),
                configs.getn("sky/patch/ysize"))
        pixelsize = configs.getn("sky/patch/pixelsize")  # [ arcsec ]
    elif not all([args.center, args.size, args.pixelsize]):
        raise ValueError("--center, --size, and --pixelsize are " +
                         "required when --config is missing!")

    if args.center:
        center = args.center.split(",")
        center = (float(center[0]), float(center[1]))
    if args.size:
        size = args.size.split(",")
        size = (int(size[0]), int(size[1]))
    if args.pixelsize:
        pixelsize = args.pixelsize

    logger.info("patch center: (%.3f, %.3f) [deg]" % center)
    logger.info("patch size: (%d, %d) pixels" % size)
    logger.info("patch pixel size: %.1f [arcsec]" % pixelsize)

    sky = SkyPatch(size=size, pixelsize=pixelsize, center=center)
    logger.info("Read HEALPix map from file: %s" % args.infile)
    hpdata, hpheader = read_fits_healpix(args.infile)
    nside = hpheader["NSIDE"]
    logger.info("Reprojecting HEALPix map to sky patch ...")
    image, __ = reproject_from_healpix(
        input_data=(hpdata, hpheader["COORDSYS"]),
        output_projection=sky.wcs, shape_out=size)
    if args.smooth:
        sigma = (2.0 * hp.nside2resol(nside, arcmin=True) *
                 AUC.arcmin2arcsec / pixelsize)
        image = scipy.ndimage.gaussian_filter(image, sigma=sigma)
        logger.info("Smoothed sky patch using Gaussian filter of " +
                    "sigma = %.2f [pixel]" % sigma)
    sky.header = hpheader.copy(strip=True)
    sky.header["OBJECT"] = "Sky Patch"
    sky.header["EXTNAME"] = "IMAGE"
    sky.header.add_history(" ".join(sys.argv))
    sky.data = image
    sky.write(args.outfile, clobber=args.clobber)
    logger.info("Write sky patch to file: %s" % args.outfile)


if __name__ == "__main__":
    main()