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