1# 2# This file is part of Healpy. 3# 4# Healpy is free software; you can redistribute it and/or modify 5# it under the terms of the GNU General Public License as published by 6# the Free Software Foundation; either version 2 of the License, or 7# (at your option) any later version. 8# 9# Healpy is distributed in the hope that it will be useful, 10# but WITHOUT ANY WARRANTY; without even the implied warranty of 11# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12# GNU General Public License for more details. 13# 14# You should have received a copy of the GNU General Public License 15# along with Healpy; if not, write to the Free Software 16# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 17# 18# For more information about Healpy, see http://code.google.com/p/healpy 19# 20""" 21===================================================== 22pixelfunc.py : Healpix pixelization related functions 23===================================================== 24 25This module provides functions related to Healpix pixelization scheme. 26 27conversion from/to sky coordinates 28---------------------------------- 29 30- :func:`pix2ang` converts pixel number to angular coordinates 31- :func:`pix2vec` converts pixel number to unit 3-vector direction 32- :func:`ang2pix` converts angular coordinates to pixel number 33- :func:`vec2pix` converts 3-vector to pixel number 34- :func:`vec2ang` converts 3-vector to angular coordinates 35- :func:`ang2vec` converts angular coordinates to unit 3-vector 36- :func:`pix2xyf` converts pixel number to coordinates within face 37- :func:`xyf2pix` converts coordinates within face to pixel number 38- :func:`get_interp_weights` returns the 4 nearest pixels for given 39 angular coordinates and the relative weights for interpolation 40- :func:`get_all_neighbours` return the 8 nearest pixels for given 41 angular coordinates 42 43conversion between NESTED and RING schemes 44------------------------------------------ 45 46- :func:`nest2ring` converts NESTED scheme pixel numbers to RING 47 scheme pixel number 48- :func:`ring2nest` converts RING scheme pixel number to NESTED 49 scheme pixel number 50- :func:`reorder` reorders a healpix map pixels from one scheme to another 51 52nside/npix/resolution 53--------------------- 54 55- :func:`nside2npix` converts healpix nside parameter to number of pixel 56- :func:`npix2nside` converts number of pixel to healpix nside parameter 57- :func:`nside2order` converts nside to order 58- :func:`order2nside` converts order to nside 59- :func:`nside2resol` converts nside to mean angular resolution 60- :func:`nside2pixarea` converts nside to pixel area 61- :func:`isnsideok` checks the validity of nside 62- :func:`isnpixok` checks the validity of npix 63- :func:`get_map_size` gives the number of pixel of a map 64- :func:`get_min_valid_nside` gives the minimum nside possible for a given 65 number of pixel 66- :func:`get_nside` returns the nside of a map 67- :func:`maptype` checks the type of a map (one map or sequence of maps) 68- :func:`ud_grade` upgrades or degrades the resolution (nside) of a map 69 70Masking pixels 71-------------- 72 73- :const:`UNSEEN` is a constant value interpreted as a masked pixel 74- :func:`mask_bad` returns a map with ``True`` where map is :const:`UNSEEN` 75- :func:`mask_good` returns a map with ``False`` where map is :const:`UNSEEN` 76- :func:`ma` returns a masked array as map, with mask given by :func:`mask_bad` 77 78Map data manipulation 79--------------------- 80 81- :func:`fit_dipole` fits a monopole+dipole on the map 82- :func:`fit_monopole` fits a monopole on the map 83- :func:`remove_dipole` fits and removes a monopole+dipole from the map 84- :func:`remove_monopole` fits and remove a monopole from the map 85- :func:`get_interp_val` computes a bilinear interpolation of the map 86 at given angular coordinates, using 4 nearest neighbours 87""" 88 89try: 90 from exceptions import NameError 91except: 92 pass 93 94import numpy as np 95from functools import wraps 96 97UNSEEN = None 98 99try: 100 from . import _healpy_pixel_lib as pixlib 101 102 #: Special value used for masked pixels 103 UNSEEN = pixlib.UNSEEN 104except: 105 import warnings 106 107 warnings.warn("Warning: cannot import _healpy_pixel_lib module") 108 109# We are using 64-bit integer types. 110# nside > 2**29 requires extended integer types. 111max_nside = 1 << 29 112 113__all__ = [ 114 "pix2ang", 115 "pix2vec", 116 "ang2pix", 117 "vec2pix", 118 "ang2vec", 119 "vec2ang", 120 "get_interp_weights", 121 "get_interp_val", 122 "get_all_neighbours", 123 "max_pixrad", 124 "nest2ring", 125 "ring2nest", 126 "reorder", 127 "ud_grade", 128 "UNSEEN", 129 "mask_good", 130 "mask_bad", 131 "ma", 132 "fit_dipole", 133 "remove_dipole", 134 "fit_monopole", 135 "remove_monopole", 136 "nside2npix", 137 "npix2nside", 138 "nside2order", 139 "order2nside", 140 "nside2resol", 141 "nside2pixarea", 142 "isnsideok", 143 "isnpixok", 144 "get_map_size", 145 "get_min_valid_nside", 146 "get_nside", 147 "maptype", 148 "ma_to_array", 149] 150 151 152def check_theta_valid(theta): 153 """Raises exception if theta is not within 0 and pi""" 154 theta = np.asarray(theta) 155 if not ((theta >= 0).all() and (theta <= np.pi + 1e-5).all()): 156 raise ValueError("THETA is out of range [0,pi]") 157 158 159def lonlat2thetaphi(lon, lat): 160 """ Transform longitude and latitude (deg) into co-latitude and longitude (rad) 161 162 Parameters 163 ---------- 164 lon : int or array-like 165 Longitude in degrees 166 lat : int or array-like 167 Latitude in degrees 168 169 Returns 170 ------- 171 theta, phi : float, scalar or array-like 172 The co-latitude and longitude in radians 173 """ 174 return np.pi / 2. - np.radians(lat), np.radians(lon) 175 176 177def thetaphi2lonlat(theta, phi): 178 """ Transform co-latitude and longitude (rad) into longitude and latitude (deg) 179 180 Parameters 181 ---------- 182 theta : int or array-like 183 Co-latitude in radians 184 phi : int or array-like 185 Longitude in radians 186 187 Returns 188 ------- 189 lon, lat : float, scalar or array-like 190 The longitude and latitude in degrees 191 """ 192 return np.degrees(phi), 90. - np.degrees(theta) 193 194 195def maptype(m): 196 """Describe the type of the map (valid, single, sequence of maps). 197 Checks : the number of maps, that all maps have same length and that this 198 length is a valid map size (using :func:`isnpixok`). 199 200 Parameters 201 ---------- 202 m : sequence 203 the map to get info from 204 205 Returns 206 ------- 207 info : int 208 -1 if the given object is not a valid map, 0 if it is a single map, 209 *info* > 0 if it is a sequence of maps (*info* is then the number of 210 maps) 211 212 Examples 213 -------- 214 >>> import healpy as hp 215 >>> hp.pixelfunc.maptype(np.arange(12)) 216 0 217 >>> hp.pixelfunc.maptype([np.arange(12), np.arange(12)]) 218 2 219 """ 220 if not hasattr(m, "__len__"): 221 raise TypeError("input map is a scalar") 222 if len(m) == 0: 223 raise TypeError("input map has length zero") 224 if hasattr(m[0], "__len__"): 225 npix = len(m[0]) 226 for mm in m[1:]: 227 if len(mm) != npix: 228 raise TypeError("input maps have different npix") 229 if isnpixok(len(m[0])): 230 return len(m) 231 else: 232 raise TypeError("bad number of pixels") 233 else: 234 if isnpixok(len(m)): 235 return 0 236 else: 237 raise TypeError("bad number of pixels") 238 239 240def ma_to_array(m): 241 """Converts a masked array or a list of masked arrays to filled numpy arrays 242 243 Parameters 244 ---------- 245 m : a map (may be a sequence of maps) 246 247 Returns 248 ------- 249 m : filled map or tuple of filled maps 250 251 Examples 252 -------- 253 >>> import healpy as hp 254 >>> m = hp.ma(np.array([2., 2., 3, 4, 5, 0, 0, 0, 0, 0, 0, 0])) 255 >>> m.mask = np.array([0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], dtype=np.bool) 256 >>> print(m.data[1]) # data is not affected by mask 257 2.0 258 >>> print(m[1]) # shows that the value is masked 259 -- 260 >>> print(ma_to_array(m)[1]) # filled array, masked values replace by UNSEEN 261 -1.6375e+30 262 """ 263 264 try: 265 return m.filled() 266 except AttributeError: 267 try: 268 return np.array([mm.filled() for mm in m]) 269 except AttributeError: 270 pass 271 return m 272 273 274def is_ma(m): 275 """Converts a masked array or a list of masked arrays to filled numpy arrays 276 277 Parameters 278 ---------- 279 m : a map (may be a sequence of maps) 280 281 Returns 282 ------- 283 is_ma : bool 284 whether the input map was a ma or not 285 """ 286 return hasattr(m, "filled") or hasattr(m[0], "filled") 287 288 289def accept_ma(f): 290 """Wraps a function in order to convert the input map from 291 a masked to a regular numpy array, and convert back the 292 output from a regular array to a masked array""" 293 294 @wraps(f) 295 def wrapper(map_in, *args, **kwds): 296 return_ma = is_ma(map_in) 297 m = ma_to_array(map_in) 298 out = f(m, *args, **kwds) 299 return ma(out) if return_ma else out 300 301 return wrapper 302 303 304def mask_bad(m, badval=UNSEEN, rtol=1.e-5, atol=1.e-8): 305 """Returns a bool array with ``True`` where m is close to badval. 306 307 Parameters 308 ---------- 309 m : a map (may be a sequence of maps) 310 badval : float, optional 311 The value of the pixel considered as bad (:const:`UNSEEN` by default) 312 rtol : float, optional 313 The relative tolerance 314 atol : float, optional 315 The absolute tolerance 316 317 Returns 318 ------- 319 mask 320 a bool array with the same shape as the input map, ``True`` where input map is 321 close to badval, and ``False`` elsewhere. 322 323 See Also 324 -------- 325 mask_good, ma 326 327 Examples 328 -------- 329 >>> import healpy as hp 330 >>> import numpy as np 331 >>> m = np.arange(12.) 332 >>> m[3] = hp.UNSEEN 333 >>> hp.mask_bad(m) 334 array([False, False, False, True, False, False, False, False, False, 335 False, False, False], dtype=bool) 336 """ 337 m = np.asarray(m) 338 atol = np.absolute(atol) 339 rtol = np.absolute(rtol) 340 return np.absolute(m - badval) <= atol + rtol * np.absolute(badval) 341 342 343def mask_good(m, badval=UNSEEN, rtol=1.e-5, atol=1.e-8): 344 """Returns a bool array with ``False`` where m is close to badval. 345 346 Parameters 347 ---------- 348 m : a map (may be a sequence of maps) 349 badval : float, optional 350 The value of the pixel considered as bad (:const:`UNSEEN` by default) 351 rtol : float, optional 352 The relative tolerance 353 atol : float, optional 354 The absolute tolerance 355 356 Returns 357 ------- 358 a bool array with the same shape as the input map, ``False`` where input map is 359 close to badval, and ``True`` elsewhere. 360 361 See Also 362 -------- 363 mask_bad, ma 364 365 Examples 366 -------- 367 >>> import healpy as hp 368 >>> m = np.arange(12.) 369 >>> m[3] = hp.UNSEEN 370 >>> hp.mask_good(m) 371 array([ True, True, True, False, True, True, True, True, True, 372 True, True, True], dtype=bool) 373 """ 374 m = np.asarray(m) 375 atol = np.absolute(atol) 376 rtol = np.absolute(rtol) 377 return np.absolute(m - badval) > atol + rtol * np.absolute(badval) 378 379 380def ma(m, badval=UNSEEN, rtol=1e-5, atol=1e-8, copy=True): 381 """Return map as a masked array, with ``badval`` pixels masked. 382 383 Parameters 384 ---------- 385 m : a map (may be a sequence of maps) 386 badval : float, optional 387 The value of the pixel considered as bad (:const:`UNSEEN` by default) 388 rtol : float, optional 389 The relative tolerance 390 atol : float, optional 391 The absolute tolerance 392 copy : bool, optional 393 If ``True``, a copy of the input map is made. 394 395 Returns 396 ------- 397 a masked array with the same shape as the input map, 398 masked where input map is close to badval. 399 400 See Also 401 -------- 402 mask_good, mask_bad, numpy.ma.masked_values 403 404 Examples 405 -------- 406 >>> import healpy as hp 407 >>> m = np.arange(12.) 408 >>> m[3] = hp.UNSEEN 409 >>> hp.ma(m) 410 masked_array(data = [0.0 1.0 2.0 -- 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0], 411 mask = [False False False True False False False False False False False False], 412 fill_value = -1.6375e+30) 413 <BLANKLINE> 414 """ 415 return np.ma.masked_values(np.array(m), badval, rtol=rtol, atol=atol, copy=copy) 416 417 418def ang2pix(nside, theta, phi, nest=False, lonlat=False): 419 """ang2pix : nside,theta[rad],phi[rad],nest=False,lonlat=False -> ipix (default:RING) 420 421 Parameters 422 ---------- 423 nside : int, scalar or array-like 424 The healpix nside parameter, must be a power of 2, less than 2**30 425 theta, phi : float, scalars or array-like 426 Angular coordinates of a point on the sphere 427 nest : bool, optional 428 if True, assume NESTED pixel ordering, otherwise, RING pixel ordering 429 lonlat : bool 430 If True, input angles are assumed to be longitude and latitude in degree, 431 otherwise, they are co-latitude and longitude in radians. 432 433 Returns 434 ------- 435 pix : int or array of int 436 The healpix pixel numbers. Scalar if all input are scalar, array otherwise. 437 Usual numpy broadcasting rules apply. 438 439 See Also 440 -------- 441 pix2ang, pix2vec, vec2pix 442 443 Examples 444 -------- 445 >>> import healpy as hp 446 >>> hp.ang2pix(16, np.pi/2, 0) 447 1440 448 449 >>> hp.ang2pix(16, [np.pi/2, np.pi/4, np.pi/2, 0, np.pi], [0., np.pi/4, np.pi/2, 0, 0]) 450 array([1440, 427, 1520, 0, 3068]) 451 452 >>> hp.ang2pix(16, np.pi/2, [0, np.pi/2]) 453 array([1440, 1520]) 454 455 >>> hp.ang2pix([1, 2, 4, 8, 16], np.pi/2, 0) 456 array([ 4, 12, 72, 336, 1440]) 457 458 >>> hp.ang2pix([1, 2, 4, 8, 16], 0, 0, lonlat=True) 459 array([ 4, 12, 72, 336, 1440]) 460 """ 461 462 check_nside(nside, nest=nest) 463 464 if lonlat: 465 theta, phi = lonlat2thetaphi(theta, phi) 466 check_theta_valid(theta) 467 check_nside(nside, nest=nest) 468 if nest: 469 return pixlib._ang2pix_nest(nside, theta, phi) 470 else: 471 return pixlib._ang2pix_ring(nside, theta, phi) 472 473 474def pix2ang(nside, ipix, nest=False, lonlat=False): 475 """pix2ang : nside,ipix,nest=False,lonlat=False -> theta[rad],phi[rad] (default RING) 476 477 Parameters 478 ---------- 479 nside : int or array-like 480 The healpix nside parameter, must be a power of 2, less than 2**30 481 ipix : int or array-like 482 Pixel indices 483 nest : bool, optional 484 if True, assume NESTED pixel ordering, otherwise, RING pixel ordering 485 lonlat : bool, optional 486 If True, return angles will be longitude and latitude in degree, 487 otherwise, angles will be longitude and co-latitude in radians (default) 488 489 Returns 490 ------- 491 theta, phi : float, scalar or array-like 492 The angular coordinates corresponding to ipix. Scalar if all input 493 are scalar, array otherwise. Usual numpy broadcasting rules apply. 494 495 See Also 496 -------- 497 ang2pix, vec2pix, pix2vec 498 499 Examples 500 -------- 501 >>> import healpy as hp 502 >>> hp.pix2ang(16, 1440) 503 (1.5291175943723188, 0.0) 504 505 >>> hp.pix2ang(16, [1440, 427, 1520, 0, 3068]) 506 (array([ 1.52911759, 0.78550497, 1.57079633, 0.05103658, 3.09055608]), array([ 0. , 0.78539816, 1.61988371, 0.78539816, 0.78539816])) 507 508 >>> hp.pix2ang([1, 2, 4, 8], 11) 509 (array([ 2.30052398, 0.84106867, 0.41113786, 0.2044802 ]), array([ 5.49778714, 5.89048623, 5.89048623, 5.89048623])) 510 511 >>> hp.pix2ang([1, 2, 4, 8], 11, lonlat=True) 512 (array([ 315. , 337.5, 337.5, 337.5]), array([-41.8103149 , 41.8103149 , 66.44353569, 78.28414761])) 513 """ 514 check_nside(nside, nest=nest) 515 if nest: 516 theta, phi = pixlib._pix2ang_nest(nside, ipix) 517 else: 518 theta, phi = pixlib._pix2ang_ring(nside, ipix) 519 520 if lonlat: 521 return thetaphi2lonlat(theta, phi) 522 else: 523 return theta, phi 524 525 526def xyf2pix(nside, x, y, face, nest=False): 527 """xyf2pix : nside,x,y,face,nest=False -> ipix (default:RING) 528 529 Parameters 530 ---------- 531 nside : int, scalar or array-like 532 The healpix nside parameter, must be a power of 2 533 x, y : int, scalars or array-like 534 Pixel indices within face 535 face : int, scalars or array-like 536 Face number 537 nest : bool, optional 538 if True, assume NESTED pixel ordering, otherwise, RING pixel ordering 539 540 Returns 541 ------- 542 pix : int or array of int 543 The healpix pixel numbers. Scalar if all input are scalar, array otherwise. 544 Usual numpy broadcasting rules apply. 545 546 See Also 547 -------- 548 pix2xyf 549 550 Examples 551 -------- 552 >>> import healpy as hp 553 >>> hp.xyf2pix(16, 8, 8, 4) 554 1440 555 556 >>> hp.xyf2pix(16, [8, 8, 8, 15, 0], [8, 8, 7, 15, 0], [4, 0, 5, 0, 8]) 557 array([1440, 427, 1520, 0, 3068]) 558 """ 559 check_nside(nside, nest=nest) 560 if nest: 561 return pixlib._xyf2pix_nest(nside, x, y, face) 562 else: 563 return pixlib._xyf2pix_ring(nside, x, y, face) 564 565 566def pix2xyf(nside, ipix, nest=False): 567 """pix2xyf : nside,ipix,nest=False -> x,y,face (default RING) 568 569 Parameters 570 ---------- 571 nside : int or array-like 572 The healpix nside parameter, must be a power of 2 573 ipix : int or array-like 574 Pixel indices 575 nest : bool, optional 576 if True, assume NESTED pixel ordering, otherwise, RING pixel ordering 577 578 Returns 579 ------- 580 x, y : int, scalars or array-like 581 Pixel indices within face 582 face : int, scalars or array-like 583 Face number 584 585 See Also 586 -------- 587 xyf2pix 588 589 Examples 590 -------- 591 >>> import healpy as hp 592 >>> hp.pix2xyf(16, 1440) 593 (8, 8, 4) 594 595 >>> hp.pix2xyf(16, [1440, 427, 1520, 0, 3068]) 596 (array([ 8, 8, 8, 15, 0]), array([ 8, 8, 7, 15, 0]), array([4, 0, 5, 0, 8])) 597 598 >>> hp.pix2xyf([1, 2, 4, 8], 11) 599 (array([0, 1, 3, 7]), array([0, 0, 2, 6]), array([11, 3, 3, 3])) 600 """ 601 check_nside(nside, nest=nest) 602 if nest: 603 return pixlib._pix2xyf_nest(nside, ipix) 604 else: 605 return pixlib._pix2xyf_ring(nside, ipix) 606 607 608def vec2pix(nside, x, y, z, nest=False): 609 """vec2pix : nside,x,y,z,nest=False -> ipix (default:RING) 610 611 Parameters 612 ---------- 613 nside : int or array-like 614 The healpix nside parameter, must be a power of 2, less than 2**30 615 x,y,z : floats or array-like 616 vector coordinates defining point on the sphere 617 nest : bool, optional 618 if True, assume NESTED pixel ordering, otherwise, RING pixel ordering 619 620 Returns 621 ------- 622 ipix : int, scalar or array-like 623 The healpix pixel number corresponding to input vector. Scalar if all input 624 are scalar, array otherwise. Usual numpy broadcasting rules apply. 625 626 See Also 627 -------- 628 ang2pix, pix2ang, pix2vec 629 630 Examples 631 -------- 632 >>> import healpy as hp 633 >>> hp.vec2pix(16, 1, 0, 0) 634 1504 635 636 >>> hp.vec2pix(16, [1, 0], [0, 1], [0, 0]) 637 array([1504, 1520]) 638 639 >>> hp.vec2pix([1, 2, 4, 8], 1, 0, 0) 640 array([ 4, 20, 88, 368]) 641 """ 642 if nest: 643 return pixlib._vec2pix_nest(nside, x, y, z) 644 else: 645 return pixlib._vec2pix_ring(nside, x, y, z) 646 647 648def pix2vec(nside, ipix, nest=False): 649 """pix2vec : nside,ipix,nest=False -> x,y,z (default RING) 650 651 Parameters 652 ---------- 653 nside : int, scalar or array-like 654 The healpix nside parameter, must be a power of 2, less than 2**30 655 ipix : int, scalar or array-like 656 Healpix pixel number 657 nest : bool, optional 658 if True, assume NESTED pixel ordering, otherwise, RING pixel ordering 659 660 Returns 661 ------- 662 x, y, z : floats, scalar or array-like 663 The coordinates of vector corresponding to input pixels. Scalar if all input 664 are scalar, array otherwise. Usual numpy broadcasting rules apply. 665 666 See Also 667 -------- 668 ang2pix, pix2ang, vec2pix 669 670 Examples 671 -------- 672 >>> import healpy as hp 673 >>> hp.pix2vec(16, 1504) 674 (0.99879545620517241, 0.049067674327418015, 0.0) 675 676 >>> hp.pix2vec(16, [1440, 427]) 677 (array([ 0.99913157, 0.5000534 ]), array([ 0. , 0.5000534]), array([ 0.04166667, 0.70703125])) 678 679 >>> hp.pix2vec([1, 2], 11) 680 (array([ 0.52704628, 0.68861915]), array([-0.52704628, -0.28523539]), array([-0.66666667, 0.66666667])) 681 """ 682 check_nside(nside, nest=nest) 683 if nest: 684 return pixlib._pix2vec_nest(nside, ipix) 685 else: 686 return pixlib._pix2vec_ring(nside, ipix) 687 688 689def ang2vec(theta, phi, lonlat=False): 690 """ang2vec : convert angles to 3D position vector 691 692 Parameters 693 ---------- 694 theta : float, scalar or arry-like 695 colatitude in radians measured southward from north pole (in [0,pi]). 696 phi : float, scalar or array-like 697 longitude in radians measured eastward (in [0, 2*pi]). 698 lonlat : bool 699 If True, input angles are assumed to be longitude and latitude in degree, 700 otherwise, they are co-latitude and longitude in radians. 701 702 Returns 703 ------- 704 vec : float, array 705 if theta and phi are vectors, the result is a 2D array with a vector per row 706 otherwise, it is a 1D array of shape (3,) 707 708 See Also 709 -------- 710 vec2ang, rotator.dir2vec, rotator.vec2dir 711 """ 712 if lonlat: 713 theta, phi = lonlat2thetaphi(theta, phi) 714 check_theta_valid(theta) 715 sintheta = np.sin(theta) 716 return np.array([sintheta * np.cos(phi), sintheta * np.sin(phi), np.cos(theta)]).T 717 718 719def vec2ang(vectors, lonlat=False): 720 """vec2ang: vectors [x, y, z] -> theta[rad], phi[rad] 721 722 Parameters 723 ---------- 724 vectors : float, array-like 725 the vector(s) to convert, shape is (3,) or (N, 3) 726 lonlat : bool, optional 727 If True, return angles will be longitude and latitude in degree, 728 otherwise, angles will be longitude and co-latitude in radians (default) 729 730 Returns 731 ------- 732 theta, phi : float, tuple of two arrays 733 the colatitude and longitude in radians 734 735 See Also 736 -------- 737 ang2vec, rotator.vec2dir, rotator.dir2vec 738 """ 739 vectors = vectors.reshape(-1, 3) 740 dnorm = np.sqrt(np.sum(np.square(vectors), axis=1)) 741 theta = np.arccos(vectors[:, 2] / dnorm) 742 phi = np.arctan2(vectors[:, 1], vectors[:, 0]) 743 phi[phi < 0] += 2 * np.pi 744 if lonlat: 745 return thetaphi2lonlat(theta, phi) 746 else: 747 return theta, phi 748 749 750def ring2nest(nside, ipix): 751 """Convert pixel number from RING ordering to NESTED ordering. 752 753 Parameters 754 ---------- 755 nside : int, scalar or array-like 756 the healpix nside parameter 757 ipix : int, scalar or array-like 758 the pixel number in RING scheme 759 760 Returns 761 ------- 762 ipix : int, scalar or array-like 763 the pixel number in NESTED scheme 764 765 See Also 766 -------- 767 nest2ring, reorder 768 769 Examples 770 -------- 771 >>> import healpy as hp 772 >>> hp.ring2nest(16, 1504) 773 1130 774 775 >>> hp.ring2nest(2, np.arange(10)) 776 array([ 3, 7, 11, 15, 2, 1, 6, 5, 10, 9]) 777 778 >>> hp.ring2nest([1, 2, 4, 8], 11) 779 array([ 11, 13, 61, 253]) 780 """ 781 check_nside(nside, nest=True) 782 return pixlib._ring2nest(nside, ipix) 783 784 785def nest2ring(nside, ipix): 786 """Convert pixel number from NESTED ordering to RING ordering. 787 788 Parameters 789 ---------- 790 nside : int, scalar or array-like 791 the healpix nside parameter 792 ipix : int, scalar or array-like 793 the pixel number in NESTED scheme 794 795 Returns 796 ------- 797 ipix : int, scalar or array-like 798 the pixel number in RING scheme 799 800 See Also 801 -------- 802 ring2nest, reorder 803 804 Examples 805 -------- 806 >>> import healpy as hp 807 >>> hp.nest2ring(16, 1130) 808 1504 809 810 >>> hp.nest2ring(2, np.arange(10)) 811 array([13, 5, 4, 0, 15, 7, 6, 1, 17, 9]) 812 813 >>> hp.nest2ring([1, 2, 4, 8], 11) 814 array([ 11, 2, 12, 211]) 815 """ 816 check_nside(nside, nest=True) 817 return pixlib._nest2ring(nside, ipix) 818 819 820@accept_ma 821def reorder(map_in, inp=None, out=None, r2n=None, n2r=None): 822 """Reorder a healpix map from RING/NESTED ordering to NESTED/RING 823 824 Parameters 825 ---------- 826 map_in : array-like 827 the input map to reorder, accepts masked arrays 828 inp, out : ``'RING'`` or ``'NESTED'`` 829 define the input and output ordering 830 r2n : bool 831 if True, reorder from RING to NESTED 832 n2r : bool 833 if True, reorder from NESTED to RING 834 835 Returns 836 ------- 837 map_out : array-like 838 the reordered map, as masked array if the input was a 839 masked array 840 841 Notes 842 ----- 843 if ``r2n`` or ``n2r`` is defined, override ``inp`` and ``out``. 844 845 See Also 846 -------- 847 nest2ring, ring2nest 848 849 Examples 850 -------- 851 >>> import healpy as hp 852 >>> hp.reorder(np.arange(48), r2n = True) 853 array([13, 5, 4, 0, 15, 7, 6, 1, 17, 9, 8, 2, 19, 11, 10, 3, 28, 854 20, 27, 12, 30, 22, 21, 14, 32, 24, 23, 16, 34, 26, 25, 18, 44, 37, 855 36, 29, 45, 39, 38, 31, 46, 41, 40, 33, 47, 43, 42, 35]) 856 >>> hp.reorder(np.arange(12), n2r = True) 857 array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]) 858 >>> hp.reorder(hp.ma(np.arange(12.)), n2r = True) 859 masked_array(data = [ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11.], 860 mask = False, 861 fill_value = -1.6375e+30) 862 <BLANKLINE> 863 >>> m = [np.arange(12.), np.arange(12.), np.arange(12.)] 864 >>> m[0][2] = hp.UNSEEN 865 >>> m[1][2] = hp.UNSEEN 866 >>> m[2][2] = hp.UNSEEN 867 >>> m = hp.ma(m) 868 >>> hp.reorder(m, n2r = True) 869 masked_array(data = 870 [[0.0 1.0 -- 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0] 871 [0.0 1.0 -- 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0] 872 [0.0 1.0 -- 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0]], 873 mask = 874 [[False False True False False False False False False False False False] 875 [False False True False False False False False False False False False] 876 [False False True False False False False False False False False False]], 877 fill_value = -1.6375e+30) 878 <BLANKLINE> 879 """ 880 typ = maptype(map_in) 881 if typ == 0: 882 npix = len(map_in) 883 else: 884 npix = len(map_in[0]) 885 nside = npix2nside(npix) 886 if nside > 128: 887 bunchsize = npix // 24 888 else: 889 bunchsize = npix 890 if r2n: 891 inp = "RING" 892 out = "NEST" 893 if n2r: 894 inp = "NEST" 895 out = "RING" 896 inp = str(inp).upper()[0:4] 897 out = str(out).upper()[0:4] 898 if inp not in ["RING", "NEST"] or out not in ["RING", "NEST"]: 899 raise ValueError("inp and out must be either RING or NEST") 900 if typ == 0: 901 mapin = [map_in] 902 else: 903 mapin = map_in 904 mapout = [] 905 for m_in in mapin: 906 if inp == out: 907 mapout.append(m_in) 908 elif inp == "RING": 909 m_out = np.zeros(npix, dtype=type(m_in[0])) 910 for ibunch in range(npix // bunchsize): 911 ipix_n = np.arange(ibunch * bunchsize, (ibunch + 1) * bunchsize) 912 ipix_r = nest2ring(nside, ipix_n) 913 m_out[ipix_n] = np.asarray(m_in)[ipix_r] 914 mapout.append(m_out) 915 elif inp == "NEST": 916 m_out = np.zeros(npix, dtype=type(m_in[0])) 917 for ibunch in range(npix // bunchsize): 918 ipix_r = np.arange(ibunch * bunchsize, (ibunch + 1) * bunchsize) 919 ipix_n = ring2nest(nside, ipix_r) 920 m_out[ipix_r] = np.asarray(m_in)[ipix_n] 921 mapout.append(m_out) 922 if typ == 0: 923 return mapout[0] 924 else: 925 return np.array(mapout) 926 927 928def nside2npix(nside): 929 """Give the number of pixels for the given nside. 930 931 Parameters 932 ---------- 933 nside : int 934 healpix nside parameter; an exception is raised if nside is not valid 935 (nside must be a power of 2, less than 2**30) 936 937 Returns 938 ------- 939 npix : int 940 corresponding number of pixels 941 942 Notes 943 ----- 944 Raise a ValueError exception if nside is not valid. 945 946 Examples 947 -------- 948 >>> import healpy as hp 949 >>> import numpy as np 950 >>> hp.nside2npix(8) 951 768 952 953 >>> np.all([hp.nside2npix(nside) == 12 * nside**2 for nside in [2**n for n in range(12)]]) 954 True 955 956 >>> hp.nside2npix(7) 957 588 958 """ 959 return 12 * nside * nside 960 961 962def nside2order(nside): 963 """Give the resolution order for a given nside. 964 965 Parameters 966 ---------- 967 nside : int 968 healpix nside parameter; an exception is raised if nside is not valid 969 (nside must be a power of 2, less than 2**30) 970 971 Returns 972 ------- 973 order : int 974 corresponding order where nside = 2**(order) 975 976 Notes 977 ----- 978 Raise a ValueError exception if nside is not valid. 979 980 Examples 981 -------- 982 >>> import healpy as hp 983 >>> import numpy as np 984 >>> hp.nside2order(128) 985 7 986 987 >>> np.all([hp.nside2order(2**o) == o for o in range(30)]) 988 True 989 990 >>> hp.nside2order(7) 991 Traceback (most recent call last): 992 ... 993 ValueError: 7 is not a valid nside parameter (must be a power of 2, less than 2**30) 994 """ 995 check_nside(nside, nest=True) 996 return len("{0:b}".format(nside)) - 1 997 998 999def nside2resol(nside, arcmin=False): 1000 """Give approximate resolution (pixel size in radian or arcmin) for nside. 1001 1002 Resolution is just the square root of the pixel area, which is a gross 1003 approximation given the different pixel shapes 1004 1005 Parameters 1006 ---------- 1007 nside : int 1008 healpix nside parameter, must be a power of 2, less than 2**30 1009 arcmin : bool 1010 if True, return resolution in arcmin, otherwise in radian 1011 1012 Returns 1013 ------- 1014 resol : float 1015 approximate pixel size in radians or arcmin 1016 1017 Notes 1018 ----- 1019 Raise a ValueError exception if nside is not valid. 1020 1021 Examples 1022 -------- 1023 >>> import healpy as hp 1024 >>> hp.nside2resol(128, arcmin = True) 1025 27.483891294539248 1026 1027 >>> hp.nside2resol(256) 1028 0.0039973699529159707 1029 1030 >>> hp.nside2resol(7) 1031 0.1461895297066412 1032 """ 1033 1034 resol = np.sqrt(nside2pixarea(nside)) 1035 1036 if arcmin: 1037 resol = np.rad2deg(resol) * 60 1038 1039 return resol 1040 1041 1042def nside2pixarea(nside, degrees=False): 1043 """Give pixel area given nside in square radians or square degrees. 1044 1045 Parameters 1046 ---------- 1047 nside : int 1048 healpix nside parameter, must be a power of 2, less than 2**30 1049 degrees : bool 1050 if True, returns pixel area in square degrees, in square radians otherwise 1051 1052 Returns 1053 ------- 1054 pixarea : float 1055 pixel area in square radian or square degree 1056 1057 Notes 1058 ----- 1059 Raise a ValueError exception if nside is not valid. 1060 1061 Examples 1062 -------- 1063 >>> import healpy as hp 1064 >>> hp.nside2pixarea(128, degrees = True) 1065 0.2098234113027917 1066 1067 >>> hp.nside2pixarea(256) 1068 1.5978966540475428e-05 1069 1070 >>> hp.nside2pixarea(7) 1071 0.021371378595848933 1072 """ 1073 1074 pixarea = 4 * np.pi / nside2npix(nside) 1075 1076 if degrees: 1077 pixarea = np.rad2deg(np.rad2deg(pixarea)) 1078 1079 return pixarea 1080 1081 1082def npix2nside(npix): 1083 """Give the nside parameter for the given number of pixels. 1084 1085 Parameters 1086 ---------- 1087 npix : int 1088 the number of pixels 1089 1090 Returns 1091 ------- 1092 nside : int 1093 the nside parameter corresponding to npix 1094 1095 Notes 1096 ----- 1097 Raise a ValueError exception if number of pixel does not correspond to 1098 the number of pixel of a healpix map. 1099 1100 Examples 1101 -------- 1102 >>> import healpy as hp 1103 >>> hp.npix2nside(768) 1104 8 1105 1106 >>> np.all([hp.npix2nside(12 * nside**2) == nside for nside in [2**n for n in range(12)]]) 1107 True 1108 1109 >>> hp.npix2nside(1000) 1110 Traceback (most recent call last): 1111 ... 1112 ValueError: Wrong pixel number (it is not 12*nside**2) 1113 """ 1114 if not isnpixok(npix): 1115 raise ValueError("Wrong pixel number (it is not 12*nside**2)") 1116 return int(np.sqrt(npix / 12.)) 1117 1118 1119def order2nside(order): 1120 """Give the nside parameter for the given resolution order. 1121 1122 Parameters 1123 ---------- 1124 order : int 1125 the resolution order 1126 1127 Returns 1128 ------- 1129 nside : int 1130 the nside parameter corresponding to order 1131 1132 Notes 1133 ----- 1134 Raise a ValueError exception if order produces an nside out of range. 1135 1136 Examples 1137 -------- 1138 >>> import healpy as hp 1139 >>> hp.order2nside(7) 1140 128 1141 1142 >>> hp.order2nside(np.arange(8)) 1143 array([ 1, 2, 4, 8, 16, 32, 64, 128]) 1144 1145 >>> hp.order2nside(31) 1146 Traceback (most recent call last): 1147 ... 1148 ValueError: 2147483648 is not a valid nside parameter (must be a power of 2, less than 2**30) 1149 """ 1150 nside = 1 << order 1151 check_nside(nside, nest=True) 1152 return nside 1153 1154 1155def isnsideok(nside, nest=False): 1156 """Returns :const:`True` if nside is a valid nside parameter, :const:`False` otherwise. 1157 1158 NSIDE needs to be a power of 2 only for nested ordering 1159 1160 Parameters 1161 ---------- 1162 nside : int, scalar or array-like 1163 integer value to be tested 1164 1165 Returns 1166 ------- 1167 ok : bool, scalar or array-like 1168 :const:`True` if given value is a valid nside, :const:`False` otherwise. 1169 1170 Examples 1171 -------- 1172 >>> import healpy as hp 1173 >>> hp.isnsideok(13, nest=True) 1174 False 1175 1176 >>> hp.isnsideok(13, nest=False) 1177 True 1178 1179 >>> hp.isnsideok(32) 1180 True 1181 1182 >>> hp.isnsideok([1, 2, 3, 4, 8, 16], nest=True) 1183 array([ True, True, False, True, True, True], dtype=bool) 1184 """ 1185 # we use standard bithacks from http://graphics.stanford.edu/~seander/bithacks.html#DetermineIfPowerOf2 1186 if hasattr(nside, "__len__"): 1187 if not isinstance(nside, np.ndarray): 1188 nside = np.asarray(nside) 1189 is_nside_ok = ( 1190 (nside == nside.astype(np.int)) & (nside > 0) & (nside <= max_nside) 1191 ) 1192 if nest: 1193 is_nside_ok &= (nside.astype(np.int) & (nside.astype(np.int) - 1)) == 0 1194 else: 1195 is_nside_ok = nside == int(nside) and 0 < nside <= max_nside 1196 if nest: 1197 is_nside_ok = is_nside_ok and (int(nside) & (int(nside) - 1)) == 0 1198 return is_nside_ok 1199 1200 1201def check_nside(nside, nest=False): 1202 """Raises exception is nside is not valid""" 1203 if not np.all(isnsideok(nside, nest=nest)): 1204 raise ValueError( 1205 "%s is not a valid nside parameter (must be a power of 2, less than 2**30)" 1206 % str(nside) 1207 ) 1208 1209 1210def isnpixok(npix): 1211 """Return :const:`True` if npix is a valid value for healpix map size, :const:`False` otherwise. 1212 1213 Parameters 1214 ---------- 1215 npix : int, scalar or array-like 1216 integer value to be tested 1217 1218 Returns 1219 ------- 1220 ok : bool, scalar or array-like 1221 :const:`True` if given value is a valid number of pixel, :const:`False` otherwise 1222 1223 Examples 1224 -------- 1225 >>> import healpy as hp 1226 >>> hp.isnpixok(12) 1227 True 1228 1229 >>> hp.isnpixok(768) 1230 True 1231 1232 >>> hp.isnpixok([12, 768, 1002]) 1233 array([ True, True, False], dtype=bool) 1234 """ 1235 nside = np.sqrt(np.asarray(npix) / 12.) 1236 return nside == np.floor(nside) 1237 1238 1239def get_interp_val(m, theta, phi, nest=False, lonlat=False): 1240 """Return the bi-linear interpolation value of a map using 4 nearest neighbours. 1241 1242 Parameters 1243 ---------- 1244 m : array-like 1245 a healpix map, accepts masked arrays 1246 theta, phi : float, scalar or array-like 1247 angular coordinates of point at which to interpolate the map 1248 nest : bool 1249 if True, the is assumed to be in NESTED ordering. 1250 lonlat : bool 1251 If True, input angles are assumed to be longitude and latitude in degree, 1252 otherwise, they are co-latitude and longitude in radians. 1253 1254 Returns 1255 ------- 1256 val : float, scalar or arry-like 1257 the interpolated value(s), usual numpy broadcasting rules apply. 1258 1259 See Also 1260 -------- 1261 get_interp_weights, get_all_neighbours 1262 1263 Examples 1264 -------- 1265 >>> import healpy as hp 1266 >>> hp.get_interp_val(np.arange(12.), np.pi/2, 0) 1267 4.0 1268 >>> hp.get_interp_val(np.arange(12.), np.pi/2, np.pi/2) 1269 5.0 1270 >>> hp.get_interp_val(np.arange(12.), np.pi/2, np.pi/2 + 2*np.pi) 1271 5.0 1272 >>> hp.get_interp_val(np.arange(12.), np.linspace(0, np.pi, 10), 0) 1273 array([ 1.5 , 1.5 , 1.5 , 2.20618428, 3.40206143, 1274 5.31546486, 7.94639458, 9.5 , 9.5 , 9.5 ]) 1275 >>> hp.get_interp_val(np.arange(12.), 0, np.linspace(90, -90, 10), lonlat=True) 1276 array([ 1.5 , 1.5 , 1.5 , 2.20618428, 3.40206143, 1277 5.31546486, 7.94639458, 9.5 , 9.5 , 9.5 ]) 1278 """ 1279 m2 = m.ravel() 1280 nside = npix2nside(m2.size) 1281 if lonlat: 1282 theta, phi = lonlat2thetaphi(theta, phi) 1283 if nest: 1284 r = pixlib._get_interpol_nest(nside, theta, phi) 1285 else: 1286 r = pixlib._get_interpol_ring(nside, theta, phi) 1287 p = np.array(r[0:4]) 1288 w = np.array(r[4:8]) 1289 del r 1290 return np.sum(m2[p] * w, 0) 1291 1292 1293def get_interp_weights(nside, theta, phi=None, nest=False, lonlat=False): 1294 """Return the 4 closest pixels on the two rings above and below the 1295 location and corresponding weights. 1296 Weights are provided for bilinear interpolation along latitude and longitude 1297 1298 Parameters 1299 ---------- 1300 nside : int 1301 the healpix nside 1302 theta, phi : float, scalar or array-like 1303 if phi is not given, theta is interpreted as pixel number, 1304 otherwise theta[rad],phi[rad] are angular coordinates 1305 nest : bool 1306 if ``True``, NESTED ordering, otherwise RING ordering. 1307 lonlat : bool 1308 If True, input angles are assumed to be longitude and latitude in degree, 1309 otherwise, they are co-latitude and longitude in radians. 1310 1311 Returns 1312 ------- 1313 res : tuple of length 2 1314 contains pixel numbers in res[0] and weights in res[1]. 1315 Usual numpy broadcasting rules apply. 1316 1317 See Also 1318 -------- 1319 get_interp_val, get_all_neighbours 1320 1321 Examples 1322 -------- 1323 >>> import healpy as hp 1324 >>> hp.get_interp_weights(1, 0) 1325 (array([0, 1, 4, 5]), array([ 1., 0., 0., 0.])) 1326 1327 >>> hp.get_interp_weights(1, 0, 0) 1328 (array([1, 2, 3, 0]), array([ 0.25, 0.25, 0.25, 0.25])) 1329 1330 >>> hp.get_interp_weights(1, 0, 90, lonlat=True) 1331 (array([1, 2, 3, 0]), array([ 0.25, 0.25, 0.25, 0.25])) 1332 1333 >>> hp.get_interp_weights(1, [0, np.pi/2], 0) 1334 (array([[ 1, 4], 1335 [ 2, 5], 1336 [ 3, 11], 1337 [ 0, 8]]), array([[ 0.25, 1. ], 1338 [ 0.25, 0. ], 1339 [ 0.25, 0. ], 1340 [ 0.25, 0. ]])) 1341 """ 1342 check_nside(nside, nest=nest) 1343 if phi is None: 1344 theta, phi = pix2ang(nside, theta, nest=nest) 1345 elif lonlat: 1346 theta, phi = lonlat2thetaphi(theta, phi) 1347 if nest: 1348 r = pixlib._get_interpol_nest(nside, theta, phi) 1349 else: 1350 r = pixlib._get_interpol_ring(nside, theta, phi) 1351 p = np.array(r[0:4]) 1352 w = np.array(r[4:8]) 1353 return (p, w) 1354 1355 1356def get_all_neighbours(nside, theta, phi=None, nest=False, lonlat=False): 1357 """Return the 8 nearest pixels. 1358 1359 Parameters 1360 ---------- 1361 nside : int 1362 the nside to work with 1363 theta, phi : scalar or array-like 1364 if phi is not given or None, theta is interpreted as pixel number, 1365 otherwise, theta[rad],phi[rad] are angular coordinates 1366 nest : bool 1367 if ``True``, pixel number will be NESTED ordering, otherwise RING ordering. 1368 lonlat : bool 1369 If True, input angles are assumed to be longitude and latitude in degree, 1370 otherwise, they are co-latitude and longitude in radians. 1371 1372 Returns 1373 ------- 1374 ipix : int, array 1375 pixel number of the SW, W, NW, N, NE, E, SE and S neighbours, 1376 shape is (8,) if input is scalar, otherwise shape is (8, N) if input is 1377 of length N. If a neighbor does not exist (it can be the case for W, N, E and S) 1378 the corresponding pixel number will be -1. 1379 1380 See Also 1381 -------- 1382 get_interp_weights, get_interp_val 1383 1384 Examples 1385 -------- 1386 >>> import healpy as hp 1387 >>> hp.get_all_neighbours(1, 4) 1388 array([11, 7, 3, -1, 0, 5, 8, -1]) 1389 1390 >>> hp.get_all_neighbours(1, np.pi/2, np.pi/2) 1391 array([ 8, 4, 0, -1, 1, 6, 9, -1]) 1392 1393 >>> hp.get_all_neighbours(1, 90, 0, lonlat=True) 1394 array([ 8, 4, 0, -1, 1, 6, 9, -1]) 1395 """ 1396 check_nside(nside, nest=nest) 1397 if phi is not None: 1398 theta = ang2pix(nside, theta, phi, nest=nest, lonlat=lonlat) 1399 if nest: 1400 r = pixlib._get_neighbors_nest(nside, theta) 1401 else: 1402 r = pixlib._get_neighbors_ring(nside, theta) 1403 res = np.array(r[0:8]) 1404 return res 1405 1406 1407def max_pixrad(nside, degrees=False): 1408 """Maximum angular distance between any pixel center and its corners 1409 1410 Parameters 1411 ---------- 1412 nside : int 1413 the nside to work with 1414 degrees : bool 1415 if True, returns pixel radius in degrees, in radians otherwise 1416 1417 Returns 1418 ------- 1419 rads: double 1420 angular distance (in radians or degrees) 1421 1422 Examples 1423 -------- 1424 >>> '%.14f' % max_pixrad(1) 1425 '0.84106867056793' 1426 >>> '%.14f' % max_pixrad(16) 1427 '0.06601476143251' 1428 """ 1429 check_nside(nside, nest=False) 1430 1431 if degrees: 1432 return np.rad2deg(pixlib._max_pixrad(nside)) 1433 1434 return pixlib._max_pixrad(nside) 1435 1436 1437def fit_dipole(m, nest=False, bad=UNSEEN, gal_cut=0): 1438 """Fit a dipole and a monopole to the map, excluding bad pixels. 1439 1440 Parameters 1441 ---------- 1442 m : float, array-like 1443 the map to which a dipole is fitted and subtracted, accepts masked maps 1444 nest : bool 1445 if ``False`` m is assumed in RING scheme, otherwise map is NESTED 1446 bad : float 1447 bad values of pixel, default to :const:`UNSEEN`. 1448 gal_cut : float [degrees] 1449 pixels at latitude in [-gal_cut;+gal_cut] degrees are not taken into account 1450 1451 Returns 1452 ------- 1453 res : tuple of length 2 1454 the monopole value in res[0] and the dipole vector (as array) in res[1] 1455 1456 See Also 1457 -------- 1458 remove_dipole, fit_monopole, remove_monopole 1459 """ 1460 m = ma_to_array(m) 1461 m = np.asarray(m) 1462 npix = m.size 1463 nside = npix2nside(npix) 1464 if nside > 128: 1465 bunchsize = npix // 24 1466 else: 1467 bunchsize = npix 1468 aa = np.zeros((4, 4), dtype=np.float64) 1469 v = np.zeros(4, dtype=np.float64) 1470 for ibunch in range(npix // bunchsize): 1471 ipix = np.arange(ibunch * bunchsize, (ibunch + 1) * bunchsize) 1472 ipix = ipix[(m.flat[ipix] != bad) & (np.isfinite(m.flat[ipix]))] 1473 x, y, z = pix2vec(nside, ipix, nest) 1474 if gal_cut > 0: 1475 w = np.abs(z) >= np.sin(gal_cut * np.pi / 180) 1476 ipix = ipix[w] 1477 x = x[w] 1478 y = y[w] 1479 z = z[w] 1480 del w 1481 aa[0, 0] += ipix.size 1482 aa[1, 0] += x.sum() 1483 aa[2, 0] += y.sum() 1484 aa[3, 0] += z.sum() 1485 aa[1, 1] += (x ** 2).sum() 1486 aa[2, 1] += (x * y).sum() 1487 aa[3, 1] += (x * z).sum() 1488 aa[2, 2] += (y ** 2).sum() 1489 aa[3, 2] += (y * z).sum() 1490 aa[3, 3] += (z ** 2).sum() 1491 v[0] += m.flat[ipix].sum() 1492 v[1] += (m.flat[ipix] * x).sum() 1493 v[2] += (m.flat[ipix] * y).sum() 1494 v[3] += (m.flat[ipix] * z).sum() 1495 aa[0, 1] = aa[1, 0] 1496 aa[0, 2] = aa[2, 0] 1497 aa[0, 3] = aa[3, 0] 1498 aa[1, 2] = aa[2, 1] 1499 aa[1, 3] = aa[3, 1] 1500 aa[2, 3] = aa[3, 2] 1501 res = np.dot(np.linalg.inv(aa), v) 1502 mono = res[0] 1503 dipole = res[1:4] 1504 return mono, dipole 1505 1506 1507def remove_dipole( 1508 m, nest=False, bad=UNSEEN, gal_cut=0, fitval=False, copy=True, verbose=True 1509): 1510 """Fit and subtract the dipole and the monopole from the given map m. 1511 1512 Parameters 1513 ---------- 1514 m : float, array-like 1515 the map to which a dipole is fitted and subtracted, accepts masked arrays 1516 nest : bool 1517 if ``False`` m is assumed in RING scheme, otherwise map is NESTED 1518 bad : float 1519 bad values of pixel, default to :const:`UNSEEN`. 1520 gal_cut : float [degrees] 1521 pixels at latitude in [-gal_cut;+gal_cut] are not taken into account 1522 fitval : bool 1523 whether to return or not the fitted values of monopole and dipole 1524 copy : bool 1525 whether to modify input map or not (by default, make a copy) 1526 verbose : bool 1527 print values of monopole and dipole 1528 1529 Returns 1530 ------- 1531 res : array or tuple of length 3 1532 if fitval is False, returns map with monopole and dipole subtracted, 1533 otherwise, returns map (array, in res[0]), monopole (float, in res[1]), 1534 dipole_vector (array, in res[2]) 1535 1536 See Also 1537 -------- 1538 fit_dipole, fit_monopole, remove_monopole 1539 """ 1540 input_ma = is_ma(m) 1541 m = ma_to_array(m) 1542 m = np.array(m, copy=copy) 1543 npix = m.size 1544 nside = npix2nside(npix) 1545 if nside > 128: 1546 bunchsize = npix // 24 1547 else: 1548 bunchsize = npix 1549 mono, dipole = fit_dipole(m, nest=nest, bad=bad, gal_cut=gal_cut) 1550 for ibunch in range(npix // bunchsize): 1551 ipix = np.arange(ibunch * bunchsize, (ibunch + 1) * bunchsize) 1552 ipix = ipix[(m.flat[ipix] != bad) & (np.isfinite(m.flat[ipix]))] 1553 x, y, z = pix2vec(nside, ipix, nest) 1554 m.flat[ipix] -= dipole[0] * x 1555 m.flat[ipix] -= dipole[1] * y 1556 m.flat[ipix] -= dipole[2] * z 1557 m.flat[ipix] -= mono 1558 if verbose: 1559 from . import rotator as R 1560 1561 lon, lat = R.vec2dir(dipole, lonlat=True) 1562 amp = np.sqrt((dipole * dipole).sum()) 1563 print( 1564 "monopole: {0:g} dipole: lon: {1:g}, lat: {2:g}, amp: {3:g}".format( 1565 mono, lon, lat, amp 1566 ) 1567 ) 1568 if is_ma: 1569 m = ma(m) 1570 if fitval: 1571 return m, mono, dipole 1572 else: 1573 return m 1574 1575 1576def fit_monopole(m, nest=False, bad=UNSEEN, gal_cut=0): 1577 """Fit a monopole to the map, excluding unseen pixels. 1578 1579 Parameters 1580 ---------- 1581 m : float, array-like 1582 the map to which a dipole is fitted and subtracted, accepts masked arrays 1583 nest : bool 1584 if ``False`` m is assumed in RING scheme, otherwise map is NESTED 1585 bad : float 1586 bad values of pixel, default to :const:`UNSEEN`. 1587 gal_cut : float [degrees] 1588 pixels at latitude in [-gal_cut;+gal_cut] degrees are not taken into account 1589 1590 Returns 1591 ------- 1592 res: float 1593 fitted monopole value 1594 1595 See Also 1596 -------- 1597 fit_dipole, remove_monopole, remove_monopole 1598 """ 1599 m = ma_to_array(m) 1600 m = np.asarray(m) 1601 npix = m.size 1602 nside = npix2nside(npix) 1603 if nside > 128: 1604 bunchsize = npix // 24 1605 else: 1606 bunchsize = npix 1607 aa = v = 0.0 1608 for ibunch in range(npix // bunchsize): 1609 ipix = np.arange(ibunch * bunchsize, (ibunch + 1) * bunchsize) 1610 ipix = ipix[(m.flat[ipix] != bad) & (np.isfinite(m.flat[ipix]))] 1611 x, y, z = pix2vec(nside, ipix, nest) 1612 if gal_cut > 0: 1613 w = np.abs(z) >= np.sin(gal_cut * np.pi / 180) 1614 ipix = ipix[w] 1615 x = x[w] 1616 y = y[w] 1617 z = z[w] 1618 del w 1619 aa += ipix.size 1620 v += m.flat[ipix].sum() 1621 mono = v / aa 1622 return mono 1623 1624 1625def remove_monopole( 1626 m, nest=False, bad=UNSEEN, gal_cut=0, fitval=False, copy=True, verbose=True 1627): 1628 """Fit and subtract the monopole from the given map m. 1629 1630 Parameters 1631 ---------- 1632 m : float, array-like 1633 the map to which a monopole is fitted and subtracted 1634 nest : bool 1635 if ``False`` m is assumed in RING scheme, otherwise map is NESTED 1636 bad : float 1637 bad values of pixel, default to :const:`UNSEEN`. 1638 gal_cut : float [degrees] 1639 pixels at latitude in [-gal_cut;+gal_cut] are not taken into account 1640 fitval : bool 1641 whether to return or not the fitted value of monopole 1642 copy : bool 1643 whether to modify input map or not (by default, make a copy) 1644 verbose: bool 1645 whether to print values of monopole 1646 1647 Returns 1648 ------- 1649 res : array or tuple of length 3 1650 if fitval is False, returns map with monopole subtracted, 1651 otherwise, returns map (array, in res[0]) and monopole (float, in res[1]) 1652 1653 See Also 1654 -------- 1655 fit_dipole, fit_monopole, remove_dipole 1656 """ 1657 input_ma = is_ma(m) 1658 m = ma_to_array(m) 1659 m = np.array(m, copy=copy) 1660 npix = m.size 1661 nside = npix2nside(npix) 1662 if nside > 128: 1663 bunchsize = npix // 24 1664 else: 1665 bunchsize = npix 1666 mono = fit_monopole(m, nest=nest, bad=bad, gal_cut=gal_cut) 1667 for ibunch in range(npix // bunchsize): 1668 ipix = np.arange(ibunch * bunchsize, (ibunch + 1) * bunchsize) 1669 ipix = ipix[(m.flat[ipix] != bad) & (np.isfinite(m.flat[ipix]))] 1670 x, y, z = pix2vec(nside, ipix, nest) 1671 m.flat[ipix] -= mono 1672 if verbose: 1673 print("monopole: {0:g}".format(mono)) 1674 if input_ma: 1675 m = ma(m) 1676 if fitval: 1677 return m, mono 1678 else: 1679 return m 1680 1681 1682def get_map_size(m): 1683 """Returns the npix of a given map (implicit or explicit pixelization). 1684 1685 If map is a dict type, assumes explicit pixelization: use nside key if present, or use 1686 nside attribute if present, otherwise use the smallest valid npix given the maximum key value. 1687 otherwise assumes implicit pixelization and returns len(m). 1688 1689 Parameters 1690 ---------- 1691 m : array-like or dict-like 1692 a map with implicit (array-like) or explicit (dict-like) pixellization 1693 1694 Returns 1695 ------- 1696 npix : int 1697 a valid number of pixel 1698 1699 Notes 1700 ----- 1701 In implicit pixellization, raise a ValueError exception if the size of the input 1702 is not a valid pixel number. 1703 1704 Examples 1705 -------- 1706 >>> import healpy as hp 1707 >>> m = {0: 1, 1: 1, 2: 1, 'nside': 1} 1708 >>> print(hp.get_map_size(m)) 1709 12 1710 1711 >>> m = {0: 1, 767: 1} 1712 >>> print(hp.get_map_size(m)) 1713 768 1714 1715 >>> print(hp.get_map_size(np.zeros(12 * 8 ** 2))) 1716 768 1717 """ 1718 if isinstance(m, dict): 1719 if "nside" in m: 1720 return nside2npix(m["nside"]) 1721 elif hasattr(ma, "nside"): 1722 return nside2npix(m.nside) 1723 else: 1724 nside = get_min_valid_nside(max(m.keys()) + 1) 1725 return nside2npix(nside) 1726 else: 1727 if isnpixok(len(m)): 1728 return len(m) 1729 else: 1730 raise ValueError("Wrong pixel number (it is not 12*nside**2)") 1731 1732 1733def get_min_valid_nside(npix): 1734 """Returns the minimum acceptable nside so that npix <= nside2npix(nside). 1735 1736 Parameters 1737 ---------- 1738 npix : int 1739 a minimal number of pixel 1740 1741 Returns 1742 ------- 1743 nside : int 1744 a valid healpix nside so that 12 * nside ** 2 >= npix 1745 1746 Examples 1747 -------- 1748 >>> import healpy as hp 1749 >>> hp.pixelfunc.get_min_valid_nside(355) 1750 8 1751 >>> hp.pixelfunc.get_min_valid_nside(768) 1752 8 1753 """ 1754 order = 0.5 * np.log2(npix / 12.) 1755 return 1 << int(np.ceil(order)) 1756 1757 1758def get_nside(m): 1759 """Return the nside of the given map. 1760 1761 Parameters 1762 ---------- 1763 m : sequence 1764 the map to get the nside from. 1765 1766 Returns 1767 ------- 1768 nside : int 1769 the healpix nside parameter of the map (or sequence of maps) 1770 1771 Notes 1772 ----- 1773 If the input is a sequence of maps, all of them must have same size. 1774 If the input is not a valid map (not a sequence, unvalid number of pixels), 1775 a TypeError exception is raised. 1776 """ 1777 typ = maptype(m) 1778 if typ == 0: 1779 return npix2nside(len(m)) 1780 else: 1781 return npix2nside(len(m[0])) 1782 1783 1784@accept_ma 1785def ud_grade( 1786 map_in, 1787 nside_out, 1788 pess=False, 1789 order_in="RING", 1790 order_out=None, 1791 power=None, 1792 dtype=None, 1793): 1794 """Upgrade or degrade resolution of a map (or list of maps). 1795 1796 in degrading the resolution, ud_grade sets the value of the superpixel 1797 as the mean of the children pixels. 1798 1799 Parameters 1800 ---------- 1801 map_in : array-like or sequence of array-like 1802 the input map(s) (if a sequence of maps, all must have same size) 1803 nside_out : int 1804 the desired nside of the output map(s) 1805 pess : bool 1806 if ``True``, in degrading, reject pixels which contains 1807 a bad sub_pixel. Otherwise, estimate average with good pixels 1808 order_in, order_out : str 1809 pixel ordering of input and output ('RING' or 'NESTED') 1810 power : float 1811 if non-zero, divide the result by (nside_in/nside_out)**power 1812 Examples: 1813 power=-2 keeps the sum of the map invariant (useful for hitmaps), 1814 power=2 divides the mean by another factor of (nside_in/nside_out)**2 1815 (useful for variance maps) 1816 dtype : type 1817 the type of the output map 1818 1819 Returns 1820 ------- 1821 map_out : array-like or sequence of array-like 1822 the upgraded or degraded map(s) 1823 1824 Examples 1825 -------- 1826 >>> import healpy as hp 1827 >>> hp.ud_grade(np.arange(48.), 1) 1828 array([ 5.5 , 7.25, 9. , 10.75, 21.75, 21.75, 23.75, 25.75, 1829 36.5 , 38.25, 40. , 41.75]) 1830 """ 1831 check_nside(nside_out, nest=order_in != "RING") 1832 typ = maptype(map_in) 1833 if typ < 0: 1834 raise TypeError("Invalid map") 1835 if typ == 0: 1836 m_in = [map_in] 1837 else: 1838 m_in = map_in 1839 mapout = [] 1840 if order_out is None: 1841 order_out = order_in 1842 for m in m_in: 1843 if str(order_in).upper()[0:4] == "RING": 1844 m = reorder(m, r2n=True) 1845 mout = _ud_grade_core(m, nside_out, pess=pess, power=power, dtype=dtype) 1846 if str(order_out).upper()[0:4] == "RING": 1847 mout = reorder(mout, n2r=True) 1848 mapout.append(mout) 1849 if typ == 0: 1850 return mapout[0] 1851 else: 1852 return np.array(mapout) 1853 1854 1855def _ud_grade_core(m, nside_out, pess=False, power=None, dtype=None): 1856 """Internal routine used by ud_grade. It assumes that the map is NESTED 1857 and single (not a list of maps) 1858 """ 1859 nside_in = get_nside(m) 1860 if dtype: 1861 type_out = dtype 1862 else: 1863 type_out = type(m[0]) 1864 check_nside(nside_out, nest=True) 1865 npix_in = nside2npix(nside_in) 1866 npix_out = nside2npix(nside_out) 1867 1868 if power: 1869 power = float(power) 1870 ratio = (float(nside_out) / float(nside_in)) ** power 1871 else: 1872 ratio = 1 1873 1874 if nside_out > nside_in: 1875 rat2 = npix_out // npix_in 1876 fact = np.ones(rat2, dtype=type_out) * ratio 1877 map_out = np.outer(m, fact).reshape(npix_out) 1878 elif nside_out < nside_in: 1879 rat2 = npix_in // npix_out 1880 mr = m.reshape(npix_out, rat2) 1881 goods = ~(mask_bad(mr) | (~np.isfinite(mr))) 1882 map_out = np.sum(mr * goods, axis=1).astype(type_out) 1883 nhit = goods.sum(axis=1) 1884 if pess: 1885 badout = np.where(nhit != rat2) 1886 else: 1887 badout = np.where(nhit == 0) 1888 if power: 1889 nhit = nhit / ratio 1890 map_out[nhit != 0] = map_out[nhit != 0] / nhit[nhit != 0] 1891 try: 1892 map_out[badout] = UNSEEN 1893 except OverflowError: 1894 pass 1895 else: 1896 map_out = m 1897 return map_out.astype(type_out) 1898