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
|
#!/usr/bin/env python3
#
# Copyright (c) 2017 Weitian LI <weitian@aaronly.me>
# MIT License
#
"""
Tile the given slices to the required FoV size, scale down to the
wanted size (for faster simulation later).
"""
import sys
import argparse
import numpy as np
from astropy.io import fits
from astropy.cosmology import FlatLambdaCDM
import astropy.units as au
import scipy.ndimage
from z2freq import z2freq, freq2z
# Adopted cosmology
cosmo = FlatLambdaCDM(H0=71.0, Om0=0.27, Ob0=0.046)
def main():
outfile_default = "{prefix}_f{freq:06.2f}_N{Nside}_fov{fov:.1f}.fits"
parser = argparse.ArgumentParser(
description="Tile slice to match FoV size and scale to required size")
parser.add_argument("-C", "--clobber", dest="clobber",
action="store_true",
help="overwrite existing files")
parser.add_argument("-F", "--fov", dest="fov", default=10.0, type=float,
help="required FoV [deg] of the output slice " +
"(default: 10.0 [deg])")
parser.add_argument("-N", "--n-side", dest="Nside", default=1800, type=int,
help="required image size of output slice " +
"(default: 1800)")
parser.add_argument("-i", "--infile", dest="infile", required=True,
help="input slice")
parser.add_argument("-o", "--outfile", dest="outfile",
default=outfile_default,
help="output tiled slice (default: %s)" %
outfile_default)
parser.add_argument("-p", "--prefix", dest="prefix", default="deltaTb",
help="prefix for the output tiled slice " +
"(default: 'deltaTb')")
exgrp = parser.add_mutually_exclusive_group(required=True)
exgrp.add_argument("-z", "--redshift-c", dest="zc", type=float,
help="central redshift of the selected cube")
exgrp.add_argument("-f", "--freq-c", dest="fc", type=float,
help="central frequency [MHz] of the selected cube")
args = parser.parse_args()
if args.zc:
zc = args.zc
fc = z2freq(zc, print_=False)
else:
fc = args.fc
zc = freq2z(fc, print_=False)
with fits.open(args.infile) as f:
img_in = f[0].data
header = f[0].header
freq = header["FREQ"] # [MHz]
Lside = header["LSIDE"] # [Mpc]
print("frequency = %.2f [MHz], Lside = %.1f [cMpc]" % (freq, Lside))
DM = cosmo.comoving_distance(zc).value # [Mpc]
print("Comoving distance (@z=%.3f/f=%.2fMHz) = %.2f [Mpc]" % (zc, fc, DM))
Nside = img_in.shape[0]
fov_in = (Lside / DM) * au.rad.to(au.deg) # [deg]
Nup = int(np.ceil(args.fov / fov_in))
print("Input FoV: %s [deg], Tiling repeats: %d" % (fov_in, Nup))
img_tiled = np.tile(img_in, reps=(Nup, Nup))
Nside2 = round(Nside * args.fov / fov_in)
img2 = img_tiled[:Nside2, :Nside2]
# Rescale to the output size
zoom = (args.Nside + 0.1) / Nside2 # +0.1 to workaround the scipy warning
img_out = scipy.ndimage.zoom(img2, zoom=zoom, order=1)
# Record information to header
header["Z_C"] = (zc, "Central redshift")
header["FREQ_C"] = (fc, "[MHz] Frequency w.r.t. to central redshift")
header["FoV"] = (args.fov, "[deg] FoV of this slice")
header["PixSize"] = (3600.0*args.fov/args.Nside, "[arcsec] Pixel size")
header.add_history(" ".join(sys.argv))
hdu = fits.PrimaryHDU(data=img_out, header=header)
outfile = args.outfile.format(prefix=args.prefix, freq=freq,
Nside=args.Nside, fov=args.fov)
try:
hdu.writeto(outfile, overwrite=args.clobber)
except TypeError:
hdu.writeto(outfile, clobber=args.clobber)
print("Tiled and scaled slice: %s" % outfile)
if __name__ == "__main__":
main()
|
