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
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
|
#!/usr/bin/env python3
#
# Copyright (c) 2017-2018 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
from fg21sim.share import CONFIGS
from fg21sim.utils import setup_logging
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", required=False,
help="fg21sim configuration from which to " +
"obtain the sky patch properties")
parser.add_argument("--center",
help="center coordinate of the sky patch; " +
"format: ra,dec; unit: deg")
parser.add_argument("--size",
help="size of the sky patch; " +
"format: xsize,ysize; unit: pixel")
parser.add_argument("--pixelsize", type=float,
help="pixel size of the sky patch; unit: [arcsec]")
parser.add_argument("-S", "--smooth", action="store_true",
help="Smooth the output patch with a Gaussian " +
"filter of sigma 'sigma-npix' (next argument) " +
"pixel size of the input HEALPix map")
parser.add_argument("--sigma-npix", dest="sigma_npix",
type=float, default=2.0,
help="number of pixels for the above Gaussian filter")
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 os.path.exists(args.outfile):
if args.clobber:
os.remove(args.outfile)
logger.warning("Removed existing output file: %s" % args.outfile)
else:
raise FileExistsError("Output file exists: %s" % args.outfile)
logger.info("Importing necessary modules, waiting ...")
import scipy.ndimage
import healpy as hp
from reproject import reproject_from_healpix
#
from fg21sim.sky import SkyPatch
from fg21sim.utils.io import read_fits_healpix
from fg21sim.utils.units import UnitConversions as AUC
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"]
try:
ordering = hpheader["ORDERING"]
except KeyError:
logger.warning("No 'ORDERING' keyword for file: %s" % args.infile)
logger.warning("Assume the 'NESTED' ordering")
ordering = "NESTED"
nested = ordering.upper() == "NESTED"
try:
coordsys = hpheader["COORDSYS"]
except KeyError:
logger.warning("No 'COORDSYS' keyword for file: %s" % args.infile)
logger.warning("Assume to use the 'Galactic' coordinate")
coordsys = "Galactic"
logger.info("Reprojecting HEALPix map to sky patch ...")
image, __ = reproject_from_healpix(input_data=(hpdata, coordsys),
output_projection=sky.wcs,
shape_out=size, nested=nested)
if args.smooth:
logger.info("Smoothing the sky patch with a Gaussian filter ...")
sigma = (args.sigma_npix * 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.merge_header(hpheader.copy(strip=True))
sky.add_history(" ".join(sys.argv))
sky.data = image
sky.write(args.outfile)
logger.info("Written extracted sky patch to file: %s" % args.outfile)
if __name__ == "__main__":
main()
|
