aboutsummaryrefslogtreecommitdiffstats
path: root/fg21sim/sky.py
diff options
context:
space:
mode:
Diffstat (limited to 'fg21sim/sky.py')
-rw-r--r--fg21sim/sky.py110
1 files changed, 102 insertions, 8 deletions
diff --git a/fg21sim/sky.py b/fg21sim/sky.py
index dbf445e..08ca0df 100644
--- a/fg21sim/sky.py
+++ b/fg21sim/sky.py
@@ -14,11 +14,12 @@ import numpy as np
from scipy import ndimage
from astropy.io import fits
import astropy.units as au
-from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
from regions import PixCoord, RectanglePixelRegion
+from reproject import reproject_interp, reproject_to_healpix
import healpy as hp
+from .utils.wcs import make_wcs
from .utils.fits import read_fits_healpix, write_fits_healpix
@@ -173,13 +174,10 @@ class SkyPatch:
--------
Currently only support the `TAN` (tangential) projection.
"""
- w = WCS(naxis=2)
- w.wcs.ctype = ["RA---TAN", "DEC--TAN"]
- w.wcs.crval = np.array([self.xcenter, self.ycenter])
- w.wcs.crpix = np.array([self.xsize/2.0-0.5, self.ysize/2.0-0.5])
- w.wcs.cdelt = np.array([-self.pixelsize/60.0, self.pixelsize/60.0])
- w.wcs.cunit = ["deg", "deg"]
- w.wcs.equinox = 2000.0
+ w = make_wcs(center=(self.xcenter, self.ycenter),
+ size=(self.xsize, self.ysize),
+ pixelsize=self.pixelsize,
+ frame="ICRS", projection="TAN")
return w
def contains(self, skycoord):
@@ -209,6 +207,50 @@ class SkyPatch:
width=self.xsize, height=self.ysize)
return region.contains(pixcoord)
+ def reproject_from(self, data, wcs, squeeze=False):
+ """
+ Reproject the given image/data together with WCS information
+ onto the grid of this sky.
+
+ Parameters
+ ----------
+ data : 2D float `~numpy.ndarray`
+ The input data/image to be reprojected
+ wcs : `~astropy.wcs.WCS`
+ The WCS information of the input data/image (naxis=2)
+ squeeze : bool, optional
+ Whether to squeeze the reprojected data to only keep
+ the positive pixels.
+
+ Returns
+ -------
+ If ``squeeze=True``, then returns tuple of ``(indexes, values)``,
+ otherwise, returns the reprojected image/data array.
+
+ indexes : 1D int `~numpy.ndarray`
+ The indexes of the pixels with positive values.
+ values : 1D float `~numpy.ndarray`
+ The values of the above pixels.
+
+ reprojected : 1D `~numpy.ndarray`
+ The reprojected data/image with same shape of this sky,
+ i.e., ``self.data``.
+ """
+ eps = 1e-5
+ wcs_out = self.wcs
+ shape_out = (self.ysize, self.xsize)
+ reprojected, __ = reproject_interp(
+ input_data=(data, wcs), output_projection=wcs_out,
+ shape_out=shape_out)
+ reprojected = reprojected.flatten()
+ if squeeze:
+ with np.errstate(invalid="ignore"):
+ indexes = reprojected > eps
+ values = reprojected[indexes]
+ return (indexes, values)
+ else:
+ return reprojected
+
class SkyHealpix:
"""
@@ -305,6 +347,58 @@ class SkyHealpix:
clobber=clobber, checksum=checksum)
logger.info("Write sky map to file: %s" % outfile)
+ def contains(self, skycoord):
+ """
+ Shim method to be consistent with ``SkyPatch``.
+
+ Always returns ``True``, since the HEALPix map covers all sky.
+ """
+ if skycoord.isscalar:
+ return True
+ else:
+ return np.ones(shape=len(skycoord), dtype=np.bool)
+
+ def reproject_from(self, data, wcs, squeeze=False):
+ """
+ Reproject the given image/data together with WCS information
+ onto the grid of this sky.
+
+ Parameters
+ ----------
+ data : 2D float `~numpy.ndarray`
+ The input data/image to be reprojected
+ wcs : `~astropy.wcs.WCS`
+ The WCS information of the input data/image (naxis=2)
+ squeeze : bool, optional
+ Whether to squeeze the reprojected data to only keep
+ the positive pixels.
+
+ Returns
+ -------
+ If ``squeeze=True``, then returns tuple of ``(indexes, values)``,
+ otherwise, returns the reprojected image/data array.
+
+ indexes : 1D int `~numpy.ndarray`
+ The indexes of the pixels with positive values.
+ values : 1D float `~numpy.ndarray`
+ The values of the above pixels.
+
+ reprojected : 1D `~numpy.ndarray`
+ The reprojected data/image with same shape of this sky,
+ i.e., ``self.data.shape``.
+ """
+ eps = 1e-5
+ reprojected, __ = reproject_to_healpix(
+ input_data=(data, wcs), coord_system_out="galactic",
+ nested=False, nside=self.nside)
+ if squeeze:
+ with np.errstate(invalid="ignore"):
+ indexes = reprojected > eps
+ values = reprojected[indexes]
+ return (indexes, values)
+ else:
+ return reprojected
+
def get_sky(configs):
"""