Source code for gammapy.utils.coordinates.celestial

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Celestial coordinate utility functions.
"""
from __future__ import absolute_import, division, print_function, unicode_literals
import numpy as np
from astropy.coordinates.angle_utilities import angular_separation

__all__ = ["minimum_separation", "pair_correlation"]


[docs]def minimum_separation(lon1, lat1, lon2, lat2): """Compute minimum distance of each (lon1, lat1) to any (lon2, lat2). Parameters ---------- lon1, lat1 : array_like Primary coordinates of interest lon2, lat2 : array_like Counterpart coordinate array unit : {'deg', 'rad'} Units of input and output coordinates Returns ------- theta_min : array Minimum distance """ lon1 = np.asanyarray(lon1) lat1 = np.asanyarray(lat1) theta_min = np.empty_like(lon1, dtype=np.float64) for i1 in range(lon1.size): thetas = angular_separation(lon1[i1], lat1[i1], lon2, lat2) theta_min[i1] = thetas.min() return theta_min
[docs]def pair_correlation(lon, lat, theta_bins): """Compute pair correlation function for points on the sphere. Parameters ---------- lon, lat : array_like Coordinate arrays theta_bins : array_like Array defining the ``theta`` binning. ``theta`` is the angular offset between positions. unit : {'deg', 'rad'} Units of input and output coordinates Returns ------- counts : array Array of point separations per ``theta`` bin. """ # TODO: Implement speedups: # - use radians # - avoid processing each pair twice (distance a to b and b to a) counts = np.zeros(shape=len(theta_bins) - 1, dtype=int) # If there are many points this should have acceptable performance # because the inner loop is in np.histogram, not in Python for ii in range(len(lon)): theta = angular_separation(lon[ii], lat[ii], lon, lat) hist = np.histogram(theta, theta_bins)[0] counts += hist return counts