Coverage for src/ufig/sampling_util.py: 88%

1# Copyright (C) 2018 ETH Zurich, Institute for Particle Physics and Astrophysics 

2 

3""" 

4Created on Mar 6, 2018 

5author: Joerg Herbel 

6""" 

7 

8import warnings 

9 

10import healpy as hp 

11import numpy as np 

12 

13from . import coordinate_util 

14 

15warnings.filterwarnings("once") 

16 

17 

18class Catalog: 

19 pass 

20 

21 

22def sample_position_uniform(numobj, w, pixel_index, nside): 

23 """ 

24 Sample a Healpix pixel uniformly 

25 

26 :param numobj: Number of uniform samples 

27 :param w: wcs-object containing all the relevant wcs-information 

28 :param pixel_index: Index of the Healpix pixels sampled 

29 :param nside: NSIDE of the Healpix map 

30 :return: Uniformly drawn x-coordinate (in pixels) 

31 :return: Uniformly drawn y-coordinate (in pixels) 

32 """ 

33 

34 corners = hp.boundaries(nside, pixel_index, 1) 

35 theta, phi = hp.vec2ang(np.transpose(corners)) 

36 

37 # removing this as it will cause sampling uniform to sample over too big range and 

38 # the while loop to take forever the subsequent code should handle this 

39 phi[phi > np.pi] -= 2 * np.pi # To deal with periodic boundaries 

40 

41 cos_theta_sample = np.random.uniform( 

42 low=np.min(np.cos(theta)), high=np.max(np.cos(theta)), size=numobj 

43 ) 

44 phi_sample = np.random.uniform(low=np.min(phi), high=np.max(phi), size=numobj) 

45 

46 n_reps = 0 

47 

48 while True: 

49 mask = hp.ang2pix(nside, np.arccos(cos_theta_sample), phi_sample) != pixel_index 

50 

51 if np.sum(mask) == 0: 

52 break 

53 

54 cos_theta_sample[mask] = np.random.uniform( 

55 low=np.min(np.cos(theta)), high=np.max(np.cos(theta)), size=np.sum(mask) 

56 ) 

57 

58 phi_sample[mask] = np.random.uniform( 

59 low=np.min(phi), high=np.max(phi), size=np.sum(mask) 

60 ) 

61 

62 n_reps += 1 

63 

64 if n_reps > 10_000: 64 ↛ 65line 64 didn't jump to line 65 because the condition on line 64 was never true

65 raise RuntimeError("more than 10_000 many iterations") 

66 if w is None: 66 ↛ 68line 66 didn't jump to line 68 because the condition on line 66 was never true

67 # If no wcs is provided, ra and dec are returned 

68 return coordinate_util.thetaphi2radec(np.arccos(cos_theta_sample), phi_sample) 

69 else: 

70 # If a wcs is provided, x and y are returned 

71 return coordinate_util.thetaphi2xy(w, np.arccos(cos_theta_sample), phi_sample)