Coverage for src/ufig/coordinate_util.py: 88%
102 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-12 19:08 +0000
1# Copyright (c) 2016 ETH Zurich, Institute of Astronomy, Claudio Bruderer
2# <claudio.bruderer@phys.ethz.ch>
4"""
5Created on Jul 5, 2016
6@author: Claudio Bruderer
7"""
9import contextlib
11import healpy as hp
12import numpy as np
13from astropy import wcs
14from cosmic_toolbox import logger
16LOGGER = logger.get_logger(__file__)
19def thetaphi2pix(theta, phi, nside, **kwargs):
20 """
21 Convert (RA, DEC) to HEALPix pixel number.
23 :param theta: Theta angle on the sphere following the HEALPix convention (in
24 radians)
25 :param phi: Phi angle on the sphere following the HEALPix convention (in radians)
26 :param nside: nside of HEALPix map
27 :param kwargs: additional keyword arguments for healpy.ang2pix-function; e.g. nest
28 :return: pixel number position lies in
29 """
31 pix = hp.ang2pix(nside=nside, theta=theta, phi=phi, **kwargs)
33 return pix
36def radec2pix(ra, dec, nside, **kwargs):
37 """
38 Convert (RA, DEC) to HEALPix pixel number.
40 :param ra: RA coordinate (in degrees)
41 :param dec: DEC coordinate (in degrees)
42 :param nside: nside of HEALPix map
43 :param kwargs: additional keyword arguments for healpy.ang2pix-function; e.g. nest
44 :return: pixel number position lies in
45 """
47 theta, phi = radec2thetaphi(ra, dec)
48 pix = thetaphi2pix(theta, phi, nside, **kwargs)
50 return pix
53def xy2pix(w, x, y, nside, **kwargs):
54 """
55 Convert (RA, DEC) to HEALPix pixel number.
57 :param w: wcs-object containing all the relevant wcs-information
58 :param x: x coordinate (in pixels)
59 :param y: y coordinate (in pixels)
60 :param nside: nside of HEALPix map
61 :param kwargs: additional keyword arguments for healpy.ang2pix-function; e.g. nest
62 :return: pixel number position lies in
63 """
65 theta, phi = xy2thetaphi(w, x, y)
66 pix = thetaphi2pix(theta, phi, nside, **kwargs)
68 return pix
71def radec2xy(w, ra, dec):
72 """
73 Convert (RA, DEC) to (x, y).
75 :param w: wcs-object containing all the relevant wcs-information
76 :param ra: RA coordinate (in degrees)
77 :param dec: DEC coordinate (in degrees)
78 :return: x coordinate (in pixels)
79 :return: y coordinate (in pixels)
80 """
82 # Check if ra and dec are empty arrays and treat this case special, since astropy
83 # cannot handle this case. The try-statement is necessary to handle single numbers
84 # as input
85 with contextlib.suppress(TypeError):
86 if len(ra) == 0:
87 x = np.array([], dtype=ra.dtype)
88 y = np.array([], dtype=dec.dtype)
89 return x, y
91 skycoords = np.array([ra, dec])
92 if np.shape(skycoords)[0] == skycoords.size:
93 skycoords = skycoords.reshape((1, skycoords.size))
94 else:
95 skycoords = skycoords.transpose()
97 x, y = np.transpose(w.wcs_world2pix(skycoords, 1))
99 # Bottom left corner is (0, 0) in UFig instead of (0.5, 0.5) as required by the
100 # WCS-transformation
101 x -= 0.5
102 y -= 0.5
104 return x, y
107def radec2thetaphi(ra, dec):
108 """
109 Convert (RA, DEC) to (theta, phi).
111 :param ra: RA coordinate (in degrees)
112 :param dec: DEC coordinate (in degrees)
113 :return: Theta angle on the sphere following the HEALPix convention (in radians)
114 :return: Phi angle on the sphere following the HEALPix convention (in radians)
115 """
117 theta = np.pi / 2 - np.pi / 180 * dec
118 phi = np.pi / 180 * ra
120 return theta, phi
123def xy2radec(w, x, y):
124 """
125 Convert (x, y) to (RA, DEC).
127 :param w: wcs-object containing all the relevant wcs-information
128 :param x: x coordinate (in pixels)
129 :param y: y coordinate (in pixels)
130 :return: RA coordinate (in degrees)
131 :return: DEC coordinate (in degrees)
132 """
134 # Check if x and y are empty arrays and treat this case special, since astropy
135 # cannot handle this case The try-statement is necessary to handle single numbers as
136 # input
137 with contextlib.suppress(TypeError):
138 if len(x) == 0:
139 ra = np.array([], dtype=x.dtype)
140 dec = np.array([], dtype=y.dtype)
141 return ra, dec
143 pixelcoords = (
144 np.array([x, y]) + 0.5
145 ) # +0.5 is to convert it into origin-1 convention
146 if np.shape(pixelcoords)[0] == pixelcoords.size:
147 pixelcoords = pixelcoords.reshape((1, pixelcoords.size))
148 else:
149 pixelcoords = pixelcoords.transpose()
151 ra, dec = np.transpose(w.wcs_pix2world(pixelcoords, 1))
153 return ra, dec
156def xy2thetaphi(w, x, y):
157 """
158 Convert (x, y) to (theta, phi).
160 :param w: wcs-object containing all the relevant wcs-information
161 :param x: x coordinate (in pixels)
162 :param y: y coordinate (in pixels)
163 :return: Theta angle on the sphere following the HEALPix convention (in radians)
164 :return: Phi angle on the sphere following the HEALPix convention (in radians)
165 """
167 ra, dec = xy2radec(w, x, y)
168 return radec2thetaphi(ra, dec)
171def thetaphi2radec(theta, phi):
172 """
173 Convert (theta, phi) to (RA, DEC).
175 :param theta: Theta angle on the sphere following the HEALPix convention (in
176 radians)
177 :param phi: Phi angle on the sphere following the HEALPix convention (in radians)
178 :return: RA coordinate (in degrees)
179 :return: DEC coordinate (in degrees)
180 """
182 ra = phi * 180.0 / np.pi
183 dec = 90.0 - theta * 180.0 / np.pi
185 return ra, dec
188def thetaphi2xy(w, theta, phi):
189 """
190 Convert (theta, phi) to (x, y).
192 :param w: wcs-object containing all the relevant wcs-information
193 :param theta: Theta angle on the sphere following the HEALPix convention (in
194 radians)
195 :param phi: Phi angle on the sphere following the HEALPix convention (in radians)
196 :return: x coordinate (in pixels)
197 :return: y coordinate (in pixels)
198 """
200 ra, dec = thetaphi2radec(theta, phi)
201 return radec2xy(w, ra, dec)
204def tile_in_skycoords(pixscale, ra0, dec0, crpix_ra, crpix_dec):
205 """
206 Maps a pixel tile onto the sky. The tile, which has pixels with width pixscale, is
207 defined around a center point (ra0, dec0). The projection employed is a tangent
208 plane gnonomic projection.
210 :param pixscale: Pixelscale of the image
211 :param size_x: size of image in x-direction (pixel)
212 :param size_y: size of image in x-direction (pixel)
213 :param ra0: Right ascension of the center of the tile
214 :param dec0: Declination of the center of the tile
215 :param crpix_ra: RA axis reference pixel (x-axis)
216 :param crpix_dec: DEC axis reference pixel (y-axis)
217 :return wcs_trans: wcs-object containing all the relevant wcs-transformation
218 information
219 """
221 c1_1 = -1 * pixscale / 60.0 / 60.0 # pixelscale in degrees
222 c1_2 = 0.0
223 c2_1 = 0.0
224 c2_2 = pixscale / 60.0 / 60.0 # pixelscale in degrees
226 # DES uses Gnonomic tangent projection as standard
227 w = wcs.WCS(naxis=2)
228 # Note, throughout UFig use convention that bottom left corner
229 # is (0, 0) --> Center is Size_x / 2
230 w.wcs.crpix = [crpix_ra, crpix_dec]
231 if c1_2 == c2_1 == 0.0: 231 ↛ 234line 231 didn't jump to line 234 because the condition on line 231 was always true
232 w.wcs.cdelt = [c1_1, c2_2]
233 else:
234 raise AssertionError(
235 "WCS transformation only implemented for a vanishing rotation at the moment"
236 )
237 w.wcs.crval = [ra0, dec0]
238 w.wcs.ctype = [b"RA---TAN", b"DEC--TAN"]
240 return w
243def get_healpix_pixels(nside, w, size_x, size_y, margin=50, npoints=300, nest=False):
244 """
245 For a given wcs-information of a tile find overlapping HEALPix pixels.
246 These pixels are found, by setting up a grid of npoints x npoints points, finding
247 the corresponding HEALPix pixels of these gridpoints and computing the unique pixel
248 IDs.
250 Note: Pixels are in RING-ordering
252 :param nside: NSIDE of the HEALPix map
253 :param w: wcs-object containing all the relevant wcs-information
254 :param size_x: size of image in x-direction (pixel)
255 :param size_y: size of image in x-direction (pixel)
256 :param margin: Number of pixels margin to ensure that the whole tile is covered
257 :param npoints: Number of points in one dimension sampled
258 :param nest: whether the HEALPix map is ordered in the NESTED-format or not;
259 default = False
260 :return pixels: HEALPix pixels overlapping a given tile
261 """
263 # The image-corners are in the origin-1 convention
264 img_corners_x = np.array([-margin, size_x + margin, -margin, size_x + margin])
265 img_corners_y = np.array([-margin, -margin, size_y + margin, size_y + margin])
266 x_vec = np.linspace(np.min(img_corners_x), np.max(img_corners_x), npoints)
267 y_vec = np.linspace(np.min(img_corners_y), np.max(img_corners_y), npoints)
269 X, Y = np.meshgrid(x_vec, y_vec)
270 x, y = X.flatten(), Y.flatten()
272 pixels = xy2pix(w, x, y, nside, nest=nest)
273 unique_pixels = np.unique(pixels)
275 return unique_pixels
278def wcs_from_parameters(par):
279 """
280 Creates an astropy WCS objects from the parameters in stored in the context.
281 This is a wrapper for tile_in_skycoords
282 to avoid code duplicates when getting crpix_ra and crpix_dec.
283 :param par: parameters, type: ivy.utils.struct.Struct
284 :return: wcs object created by tile_in_skycoords
285 """
287 try:
288 crpix_ra = par.crpix_ra
289 except AttributeError:
290 crpix_ra = par.size_x / 2 + 0.5
291 LOGGER.debug("Setting RA reference pixel to the center of the field")
293 try:
294 crpix_dec = par.crpix_dec
295 except AttributeError:
296 crpix_dec = par.size_y / 2 + 0.5
297 LOGGER.debug("Setting DEC reference pixel to the center of the field")
299 w = tile_in_skycoords(
300 pixscale=par.pixscale,
301 ra0=par.ra0,
302 dec0=par.dec0,
303 crpix_ra=crpix_ra,
304 crpix_dec=crpix_dec,
305 )
307 return w
310def get_healpix_pixels_from_map(par):
311 """
312 Returns the pixels of a boolean HEALPix map that are True.
313 Also adapts the nside_sampling parameter if it is not set to the same value as
314 the nside of the map.
315 """
316 if par.healpix_map is None:
317 raise ValueError(
318 "The healpix_map parameter is not set."
319 " Use sampling_mode=wcs or define a healpix_map."
320 )
322 pixels = np.where(par.healpix_map)[0]
324 nside_sampling = hp.get_nside(par.healpix_map)
325 if nside_sampling != par.nside_sampling:
326 LOGGER.warning(
327 "The nside_sampling parameter is different from the nside of the map. "
328 "Setting nside_sampling to the nside of the map."
329 )
330 par.nside_sampling = nside_sampling
332 return pixels