aboutsummaryrefslogtreecommitdiffstats
path: root/fg21sim
diff options
context:
space:
mode:
Diffstat (limited to 'fg21sim')
-rw-r--r--fg21sim/configs/checkers.py2
-rw-r--r--fg21sim/galactic/snr.py220
2 files changed, 119 insertions, 103 deletions
diff --git a/fg21sim/configs/checkers.py b/fg21sim/configs/checkers.py
index 2535f9e..1ba1f61 100644
--- a/fg21sim/configs/checkers.py
+++ b/fg21sim/configs/checkers.py
@@ -167,7 +167,7 @@ def check_galactic_snr(configs):
sky = get_sky(configs)
key = comp + "/resolution"
resolution = configs.getn(key) # [ arcmin ]
- if resolution < sky.pixelsize:
+ if resolution > sky.pixelsize:
results[key] = "resolution should be higher than output map"
return results
diff --git a/fg21sim/galactic/snr.py b/fg21sim/galactic/snr.py
index 088a2ae..f715aa0 100644
--- a/fg21sim/galactic/snr.py
+++ b/fg21sim/galactic/snr.py
@@ -1,4 +1,4 @@
-# Copyright (c) 2016 Weitian LI <liweitianux@live.com>
+# Copyright (c) 2016-2017 Weitian LI <weitian@aaronly.me>
# MIT license
"""
@@ -11,13 +11,14 @@ from datetime import datetime, timezone
import numpy as np
from astropy.io import fits
+from astropy.coordinates import SkyCoord
import astropy.units as au
-import healpy as hp
import pandas as pd
-from ..utils.fits import write_fits_healpix
+from ..sky import get_sky
+from ..utils.wcs import make_wcs
from ..utils.convert import Fnu_to_Tb_fast
-from ..utils.grid import make_grid_ellipse, map_grid_to_healpix
+from ..utils.grid import make_ellipse
logger = logging.getLogger(__name__)
@@ -71,6 +72,7 @@ class SuperNovaRemnants:
def __init__(self, configs):
self.configs = configs
+ self.sky = get_sky(configs)
self._set_configs()
def _set_configs(self):
@@ -78,7 +80,7 @@ class SuperNovaRemnants:
comp = "galactic/snr"
self.catalog_path = self.configs.get_path(comp+"/catalog")
self.catalog_outfile = self.configs.get_path(comp+"/catalog_outfile")
- self.resolution = self.configs.getn(comp+"/resolution") * au.arcmin
+ self.resolution = self.configs.getn(comp+"/resolution") # [ arcmin ]
self.prefix = self.configs.getn(comp+"/prefix")
self.save = self.configs.getn(comp+"/save")
self.output_dir = self.configs.get_path(comp+"/output_dir")
@@ -87,32 +89,30 @@ class SuperNovaRemnants:
self.use_float = self.configs.getn("output/use_float")
self.checksum = self.configs.getn("output/checksum")
self.clobber = self.configs.getn("output/clobber")
- self.nside = self.configs.getn("common/nside")
- self.pixsize = hp.nside2resol(self.nside, arcmin=True) / 60.0
- self.pixarea = self.pixsize ** 2 # [ deg^2 ]
self.freq_unit = au.Unit(self.configs.getn("frequency/unit"))
logger.info("Loaded and set up configurations")
def _load_catalog(self):
- """Load the Galactic SNRs catalog data."""
+ """
+ Load the Galactic SNRs catalog data.
+
+ Catalog column units:
+ * glon, glat : deg
+ * size_major, size_minor : arcmin
+ * flux : Jy
+ """
self.catalog = pd.read_csv(self.catalog_path)
nrow, ncol = self.catalog.shape
logger.info("Loaded SNRs catalog data from: {0}".format(
self.catalog_path))
logger.info("SNRs catalog data: {0} objects, {1} columns".format(
nrow, ncol))
- # Set the units for columns
- self.units = {
- "glon": au.deg,
- "glat": au.deg,
- "size": au.arcmin,
- "flux": au.Jy,
- }
# The flux densities are given at 1 GHz
self.catalog_flux_freq = (1.0*au.GHz).to(self.freq_unit).value
def _save_catalog_inuse(self):
- """Save the effective/inuse SNRs catalog data to a CSV file.
+ """
+ Save the effective/inuse SNRs catalog data to a CSV file.
NOTE
----
@@ -146,8 +146,9 @@ class SuperNovaRemnants:
logger.info("Saved SNRs catalog in use to: %s" % self.catalog_outfile)
def _filter_catalog(self):
- """Filter the catalog data to remove the objects with incomplete
- data.
+ """
+ Filter the catalog data to remove the objects with incomplete data,
+ as well as the SNRs lying outside the sky coverage.
The following cases are filtered out:
- Missing angular size
@@ -163,20 +164,31 @@ class SuperNovaRemnants:
cond3 = pd.isnull(self.catalog["flux"])
cond4 = pd.isnull(self.catalog["specindex"])
cond_keep = ~(cond1 | cond2 | cond3 | cond4)
- n_total = len(cond_keep)
n_remain = cond_keep.sum()
- n_delete = n_total - n_remain
- n_delete_p = n_delete / n_total * 100
+ n_delete = len(cond_keep) - n_remain
self.catalog = self.catalog[cond_keep]
+ logger.info("SNRs catalog: filtered out due to incomplete data: " +
+ "{0:d} objects".format(n_delete))
+ # Filter out the SNRs lying outside the sky region (e.g., a patch)
+ skycoords = SkyCoord(l=self.catalog["glon"], b=self.catalog["glat"],
+ frame="galactic", unit="deg")
+ inside = self.sky.contains(skycoords)
+ n_remain = inside.sum()
+ n_delete = len(inside) - n_remain
+ self.catalog = self.catalog[inside]
# Drop the index
self.catalog.reset_index(drop=True, inplace=True)
self.catalog_filtered = True
- logger.info("SNRs catalog: filtered out " +
- "{0:d} ({1:.1f}%) objects".format(n_delete, n_delete_p))
+ logger.info("SNRs catalog: filtered out due to sky coverage: " +
+ "{0:d} objects".format(n_delete))
logger.info("Filtered SNRs catalog: {0} objects".format(n_remain))
+ if n_remain == 0:
+ raise RuntimeError("NO remaining SNRs within simulation sky! " +
+ "Check the catalog or disable this component.")
def _add_random_rotation(self):
- """Add random rotation angles for each SNR as column "rotation"
+ """
+ Add random rotation angles for each SNR as column "rotation"
within the catalog data frame.
The rotation angles are uniformly distributed within [0, 360).
@@ -188,11 +200,11 @@ class SuperNovaRemnants:
angles = np.random.uniform(low=0.0, high=360.0, size=num)
rotation = pd.Series(data=angles, name="rotation")
self.catalog["rotation"] = rotation
- self.units["rotation"] = au.deg
logger.info("Added random rotation angles as the 'rotation' column")
def _calc_Tb(self, flux, specindex, frequency, size):
- """Calculate the brightness temperature at requested frequency
+ """
+ Calculate the brightness temperature at requested frequency
by assuming a power-law spectral shape.
Parameters
@@ -225,58 +237,63 @@ class SuperNovaRemnants:
freq_ref = self.catalog_flux_freq # [ MHz ]
Fnu = flux * (frequency / freq_ref) ** (-specindex) # [ Jy ]
omega = size[0] * size[1] # [ deg^2 ]
- if omega < self.pixarea:
+ pixelarea = (self.sky.pixelsize/60.0) ** 2 # [ deg^2 ]
+ if omega < pixelarea:
# The object is smaller than a pixel, so round up to a pixel area
- omega = self.pixarea
+ omega = pixelarea
Tb = Fnu_to_Tb_fast(Fnu, omega, frequency)
return Tb
def _simulate_templates(self):
- """Simulate the template (HEALPix) images for each SNR, and cache
- these templates within the class.
+ """
+ Simulate the template images/maps for each SNR, and cache these
+ templates within this object.
The template images/maps have values of (or approximate) ones for
these effective pixels, excluding the pixels corresponding to the
edges of original rotated ellipse, which may have values of
significantly less than 1 due to the rotation.
- Therefore, simulating the HEALPix map of one SNR at a requested
- frequency is simply multiplying the cached template image by the
+ Therefore, simulating the sky map of one SNR at a requested
+ frequency is simply multiplying these cached templates by the
calculated brightness temperature (Tb) at that frequency.
- Furthermore, the total HEALPix map of all SNRs are straightforward
- additions of all the maps of each SNR.
+ Furthermore, the final output sky map of all SNRs are just additions
+ of all the maps of each SNR.
Attributes
----------
templates : dict
A dictionary containing the simulated templates for each SNR.
The dictionary keys are the names (`self.catalog["name"]`)
- of the SNRs, and the values are `(hpidx, hpval)` tuples with
- `hpidx` the indexes of effective HEALPix pixels (RING ordering)
- and `hpval` the values of the corresponding pixels.
+ of the SNRs, and the values are `(idx, val)` tuples with
+ `idx` the indexes of effective image pixels and `hpval` the
+ values of the corresponding pixels.
e.g.,
- ``{ name1: (hpidx1, hpval1), name2: (hpidx2, hpval2), ... }``
+ ``{ name1: (idx1, val1), name2: (idx2, val2), ... }``
"""
templates = {}
- resolution = self.resolution.to(au.deg).value
- logger.info("Simulate HEALPix template for each SNR")
- coef_size2deg = self.units["size"].to(au.deg)
+ logger.info("Simulate sky template for each SNR")
for row in self.catalog.itertuples():
name = row.name
- logger.debug("Simulate HEALPix template for SNR: {0}".format(name))
- center = (row.glon, row.glat)
- size = (row.size_major * coef_size2deg,
- row.size_minor * coef_size2deg)
- rotation = row.rotation
- grid = make_grid_ellipse(center, size, resolution, rotation)
- hpidx, hpval = map_grid_to_healpix(grid, self.nside)
- templates[name] = (hpidx, hpval)
+ logger.debug("Simulate sky template for SNR: {0}".format(name))
+ gcenter = (row.glon, row.glat) # [ deg ]
+ radii = (int(np.ceil(row.size_major * 0.5 / self.resolution)),
+ int(np.ceil(row.size_minor * 0.5 / self.resolution)))
+ rmax = max(radii)
+ pcenter = (rmax, rmax)
+ image = make_ellipse(pcenter, radii, row.rotation)
+ wcs = make_wcs(center=gcenter, size=image.shape,
+ pixelsize=self.resolution,
+ frame="Galactic", projection="CAR")
+ idx, val = self.sky.reproject_from(image, wcs, squeeze=True)
+ templates[name] = (idx, val)
logger.info("Done simulate {0} SNR templates".format(len(templates)))
self.templates = templates
def _simulate_single(self, data, frequency):
- """Simulate one single SNR at the specified frequency.
+ """
+ Simulate one single SNR at the specified frequency.
Parameters
----------
@@ -290,31 +307,29 @@ class SuperNovaRemnants:
Returns
-------
- hpidx : 1D `~numpy.ndarray`
- The indexes (in RING ordering) of the effective HEALPix
- pixels for the SNR.
- hpval : 1D `~numpy.ndarray`
- The values (i.e., brightness temperature) of each HEALPix
- pixel with respect the above indexes.
+ idx : 1D `~numpy.ndarray`
+ The indexes of the effective map pixels for the SNR.
+ val : 1D `~numpy.ndarray`
+ The values (i.e., brightness temperature) of each map
+ pixel with respect to the above indexes.
See Also
--------
`self._simulate_template()` for more detailed description.
"""
name = data.name
- hpidx, hpval = self.templates[name]
+ idx, val = self.templates[name]
# Calculate the brightness temperature
flux = data.flux
specindex = data.specindex
- coef_size2deg = self.units["size"].to(au.deg)
- size = (data.size_major * coef_size2deg,
- data.size_minor * coef_size2deg) # [ deg ]
+ size = (data.size_major/60.0, data.size_minor/60.0) # [ deg ]
Tb = self._calc_Tb(flux, specindex, frequency, size)
- hpval = hpval * Tb
- return (hpidx, hpval)
+ val = val * Tb
+ return (idx, val)
def _make_filepath(self, **kwargs):
- """Make the path of output file according to the filename pattern
+ """
+ Make the path of output file according to the filename pattern
and output directory loaded from configurations.
"""
data = {
@@ -326,12 +341,12 @@ class SuperNovaRemnants:
return filepath
def _make_header(self):
- """Make the header with detail information (e.g., parameters and
+ """
+ Make the header with detail information (e.g., parameters and
history) for the simulated products.
"""
header = fits.Header()
- header["COMP"] = ("Galactic supernova remnants (SNRs)",
- "Emission component")
+ header["COMP"] = (self.name, "Emission component")
header["UNIT"] = ("Kelvin", "Map unit")
header["CREATOR"] = (__name__, "File creator")
# TODO:
@@ -344,20 +359,17 @@ class SuperNovaRemnants:
self.header = header
logger.info("Created FITS header")
- def output(self, hpmap, frequency):
- """Write the simulated free-free map to disk with proper header
+ def output(self, skymap, frequency):
+ """
+ Write the simulated free-free map to disk with proper header
keywords and history.
Returns
-------
- filepath : str
- The (absolute) path to the output HEALPix map file.
+ outfile : str
+ The (absolute) path to the output sky map file.
"""
- if not os.path.exists(self.output_dir):
- os.mkdir(self.output_dir)
- logger.info("Created output dir: {0}".format(self.output_dir))
- #
- filepath = self._make_filepath(frequency=frequency)
+ outfile = self._make_filepath(frequency=frequency)
if not hasattr(self, "header"):
self._make_header()
header = self.header.copy()
@@ -367,14 +379,16 @@ class SuperNovaRemnants:
"File creation date"
)
if self.use_float:
- hpmap = hpmap.astype(np.float32)
- write_fits_healpix(filepath, hpmap, header=header,
- clobber=self.clobber, checksum=self.checksum)
- logger.info("Write simulated map to file: {0}".format(filepath))
- return filepath
+ skymap = skymap.astype(np.float32)
+ sky = self.sky.copy()
+ sky.data = skymap
+ sky.header = header
+ sky.write(outfile, clobber=self.clobber, checksum=self.checksum)
+ return outfile
def preprocess(self):
- """Perform the preparation procedures for the final simulations.
+ """
+ Perform the preparation procedures for the final simulations.
Attributes
----------
@@ -396,8 +410,9 @@ class SuperNovaRemnants:
self._preprocessed = True
def simulate_frequency(self, frequency):
- """Simulate the emission (HEALPix) map of all Galactic SNRs at
- the specified frequency.
+ """
+ Simulate the sky map of all Galactic SNRs emission at the specified
+ frequency.
Parameters
----------
@@ -406,10 +421,10 @@ class SuperNovaRemnants:
Returns
-------
- hpmap_f : 1D `~numpy.ndarray`
- The HEALPix map (RING ordering) at the input frequency.
+ skymap_f : 1D/2D `~numpy.ndarray`
+ The sky map at the input frequency.
filepath : str
- The (absolute) path to the output HEALPix file if saved,
+ The (absolute) path to the output sky map if saved,
otherwise ``None``.
See Also
@@ -420,20 +435,21 @@ class SuperNovaRemnants:
#
logger.info("Simulating {name} map at {freq} ({unit}) ...".format(
name=self.name, freq=frequency, unit=self.freq_unit))
- hpmap_f = np.zeros(hp.nside2npix(self.nside))
+ skymap_f = np.zeros(self.sky.shape)
for row in self.catalog.itertuples():
- hpidx, hpval = self._simulate_single(row, frequency)
- hpmap_f[hpidx] += hpval
+ index, value = self._simulate_single(row, frequency)
+ skymap_f[index] += value
#
if self.save:
- filepath = self.output(hpmap_f, frequency)
+ filepath = self.output(skymap_f, frequency)
else:
filepath = None
- return (hpmap_f, filepath)
+ return (skymap_f, filepath)
def simulate(self, frequencies):
- """Simulate the emission (HEALPix) maps of all Galactic SNRs for
- every specified frequency.
+ """
+ Simulate the sky maps of all Galactic SNRs emission at every
+ specified frequency.
Parameters
----------
@@ -443,18 +459,18 @@ class SuperNovaRemnants:
Returns
-------
- hpmaps : list[1D `~numpy.ndarray`]
- List of HEALPix maps (in RING ordering) at each frequency.
+ skymaps : list[1D `~numpy.ndarray`]
+ List of sky maps at each frequency.
paths : list[str]
- List of (absolute) path to the output HEALPix maps.
+ List of (absolute) path to the output sky maps.
"""
- hpmaps = []
+ skymaps = []
paths = []
for f in np.array(frequencies, ndmin=1):
- hpmap_f, filepath = self.simulate_frequency(f)
- hpmaps.append(hpmap_f)
- paths.append(filepath)
- return (hpmaps, paths)
+ skymap_f, outfile = self.simulate_frequency(f)
+ skymaps.append(skymap_f)
+ paths.append(outfile)
+ return (skymaps, paths)
def postprocess(self):
"""Perform the post-simulation operations before the end."""