aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorAaron LI <aly@aaronly.me>2017-06-12 14:04:19 +0800
committerAaron LI <aly@aaronly.me>2017-06-12 14:04:19 +0800
commit086f5a5ab0c27162b3d8a725b4de19bc35679a2f (patch)
tree4061904f4fe8b6d0bc34e9961ac488925715c33e
parentf79809f04ed784e802b94e23529ab86c124c1f1d (diff)
downloadfg21sim-086f5a5ab0c27162b3d8a725b4de19bc35679a2f.tar.bz2
bin/get-healpix-patch: Add argument --smooth and write history
-rwxr-xr-xbin/get-healpix-patch13
1 files changed, 13 insertions, 0 deletions
diff --git a/bin/get-healpix-patch b/bin/get-healpix-patch
index 8295d0b..7f01f78 100755
--- a/bin/get-healpix-patch
+++ b/bin/get-healpix-patch
@@ -14,6 +14,8 @@ import sys
import argparse
import logging
+import scipy.ndimage
+import healpy as hp
from reproject import reproject_from_healpix
from fg21sim.configs import configs
@@ -38,6 +40,10 @@ def main():
"format: xsize,ysize; unit: pixel")
parser.add_argument("--pixelsize", dest="pixelsize", type=float,
help="pixel size of the sky patch; unit: arcmin")
+ 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()
@@ -74,13 +80,20 @@ def main():
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 * hp.nside2resol(nside, arcmin=True) / 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)