Source code for cosmo_utils.utils.geometry
#! /usr/bin/env python
# -*- coding: utf-8 -*-
# Victor Calderon
# Created : 2018-05-03
# Last Modified: 2018-05-03
from __future__ import absolute_import, division, print_function
__author__ = ['Victor Calderon']
__copyright__ = ["Copyright 2018 Victor Calderon"]
__email__ = ['victor.calderon@vanderbilt.edu']
__maintainer__ = ['Victor Calderon']
__all__ = [ "flip_angles",
"Ang_Distance",
"Coord_Transformation"]
"""
Set of geometrical definitions for translations, coordinate tranformations,
etc.
"""
## Import modules
import numpy as np
import pandas as pd
from cosmo_utils.utils import file_utils as fd
from cosmo_utils.custom_exceptions import LSSUtils_Error
## Functions
# Restricting angles to be between 0 and 360
[docs]def flip_angles(ang, unit='deg'):
"""
Ensures that an angle is always between 0 and 360 degrees.
Parameters
----------
ang : float, int, `numpy.ndarray`
Angle in units of degrees.
unit : {'deg', 'rad'} str, optional
Unit of the angle. This variable is set to 'deg' by default.
If 'rad', the output will be in units of radians.
Returns
----------
ang_out : float
Convertions of `ang`, ranging from 0 to 360 degrees.
Raises
----------
LSSUtils_Error : Exception
This error gets raised when `ang` is not a digit or a number.
Examples
----------
>>> flip_angles(-50, unit='deg')
310.0
>>> flip_angles(110, unit='deg')
110.0
"""
file_msg = fd.Program_Msg(__file__)
# Checking type of `ang`
valid_types = (float, int, list, np.ndarray)
if not (isinstance(ang, valid_types)):
msg = '{0} `ang` ({1}) is not a number! Exiting!'.format(
file_msg, ang)
raise LSSUtils_Error(msg)
# Checking `unit`
if not (unit.lower() in ['deg', 'rad']):
msg = '{0} `unit` ({1}) is not a valid option! Exiting!'.format(
file_msg, unit)
raise LSSUtils_Error(msg)
##
## Converting angle
if isinstance(ang, float) or isinstance(ang, int):
ang = float(ang)
# Converting from radians to degrees, if applicable
if unit == 'rad':
ang_converted = np.degrees(ang)
elif unit == 'deg':
ang_converted = ang
# Checking the range of `ang`
if ang_converted < 0:
ang_converted += 360
# Converting back to radians, if applicable
if unit == 'rad':
ang_final = np.radians(ang_converted)
elif unit == 'deg':
ang_final = float(ang_converted)
else:
try:
ang = np.asarray(ang)
# Converting to degrees, if applicable
if unit == 'rad':
ang_converted = np.degrees(ang)
elif unit == 'deg':
ang_converted = ang
# Checking the range of `ang`
ang_converted = np.asarray([xx if xx > 0 else xx+360 for xx in
ang_converted])
# Converting back to radians, if applicable
if unit == 'rad':
ang_final = np.radians(ang_converted)
elif unit == 'deg':
ang_final = ang_converted
# Converting to float
ang_final = ang_final.astype(float)
except:
msg = '{0} `ang` could not be converted!'.format(file_msg)
raise LSSUtils_Error(msg)
return ang_final
## Calculates the Angular distance between 2 points
[docs]def Ang_Distance(ra1, ra2, dec1, dec2, unit='deg', method='haversine'):
"""
Calculates angular separation between two sets of points with given
right ascensions and declinations.
Taken from: https://en.wikipedia.org/wiki/Haversine_formula
Parameters
-----------
ra1, ra2 : float
Right Ascension of the 1st and 2nd points. Units in `degrees`
by default.
dec1, dec2 : float
Declination of the 1st and 2nd points. Units in `degrees` by default.
unit : {'dec','rad'} str, optional
Unit of `ra1`, `ra2`, `dec1`, and `dec2`.
This will also determine the final unit that outputs this function.
method : {'haversine', 'astropy'} str, optional
Method to use in order to calculate angular separation.
This variable is to by default to the `haversine` method.
If `astropy`, it will use the astropy framework to determine the
angular separation.
Returns
-----------
ang_sep : float
Angular separation between 1st and 2nd point.
In units of `degrees`.
Notes
-----------
A = 90. - `dec2`
B = 90. - `dec1`
D = `ra1` - `ra2`
c = Angle between two points
"""
file_msg = fd.Program_Msg(__file__)
## Checking input arguments
# Valid options
units_valid = ['deg', 'rad']
methods_valid = ['haversine', 'astropy']
# Units
if not (unit in units_valid):
msg = '{0} `unit` ({1}) is not a valid argument'.format(
file_msg, unit)
raise LSSUtils_Error(msg)
# Method
if not (method in methods_valid):
msg = '{0} `method` ({1}) is not a valid argument'.format(
file_msg, method)
raise LSSUtils_Error(msg)
##
## Flipping angles
ra1 = flip_angles(ra1, unit=unit, )
ra2 = flip_angles(ra2, unit=unit)
##
## Haversine Method
if method == 'haversine':
A = np.radians(90. - dec1)
B = np.radians(90. - dec2)
D = np.radians(ra1 - ra2 )
# Distance
ang_sep = (np.sin((A - B) * .5))**2.
ang_sep += np.sin(A) * np.sin(B) * (np.sin(D * .5))**2.
ang_sep = np.degrees(2 * np.arcsin(ang_sep**0.5))
##
## Astropy Method
if method == 'astropy':
# Imports
from astropy.coordinates import SkyCoord
from astropy import units as u
# Converting to `units`
if unit == 'deg':
unit_opt = u.degree
elif unit == 'rad':
unit_opt = u.radians
# Calculations
p1 = SkyCoord(ra=ra1, dec=dec1, unit=(unit_opt, unit_opt))
p2 = SkyCoord(ra=ra2, dec=dec2, unit=(unit_opt, unit_opt))
# Angular Separation
ang_sep = p1.separation(p2)
# Converting to final units
if unit == 'deg':
ang_sep = ang_sep.degree
elif unit == 'rad':
ang_sep = ang_sep.radians
return ang_sep
## Coordinate Transformation function
[docs]def Coord_Transformation(ra, dec, dist, ra_cen, dec_cen, dist_cen,
trans_opt=4, return_dict=False, unit='deg'):
"""
Transforms spherical coordinates (ra, dec, dist) into cartesian
coordinates.
Parameters
-----------
ra, dec, dist : array_like, shape (N,)
Arrays of Right Ascension, declination, and distance.
Units are ['degrees', 'degrees', 'distance_units']
ra_cen, dec_cen, dist_cen : float, int
Right Ascension, declination, and distance for the center of
the coordinates. These correspond to where the corodinates
`ra`, `dec`, and `dist` will be centered.
trans_opt : {1, 2, 3, 4} int, optional
Option for cartesian translation/transformation for elements.
This variable ist set to `4` by default.
Options:
- 1 : No translation involved
- 2 : Translation to the center point.
- 3 : Translation `and` rotation to the center point.
- 4 : Translation and 2 rotaitons about the center point
return_dict : {True, False}, `bool`, optional
If `True`, this functions returns 2 dictionaries with `spherical`
and `cartesian` coordinates.
If `False`, it returns a `pandas.DataFrame` with the columns.
This variable is set to `False` by default.
unit : {'dec','rad'} str, optional
Unit of `ra1`, `ra2`, `dec1`, and `dec2`.
This will also determine the final unit that outputs this function.
This variable is set to `deg` by default.
Returns
-----------
coord_dict (coord_pd) : python dictionary
Dictionary with spherical and cartesian dictionary of elements
based on `trans_opt` value. This value is returned if
`return_dict` is set to `True`. If not, a `pandas.DataFrame` is
return.
"""
file_msg = fd.Program_Msg(__file__)
## Check types of elements
# Units
unit_arr = ['deg', 'rad']
if not (unit in unit_arr):
'{0} `unit` ({1}) is not a valid input!'.format(
file_msg, unit)
# Valid types
valid_types = (float, int, np.ndarray, list)
# Right Ascension
if not (isinstance(ra, valid_types)):
msg = '{0} `ra` ({1}) is not a valid type!'.format(
file_msg, type(ra))
raise LSSUtils_Error(msg)
# Declination
if not (isinstance(dec, valid_types)):
msg = '{0} `dec` ({1}) is not a valid type!'.format(
file_msg, type(dec))
raise LSSUtils_Error(msg)
# Distance
if not (isinstance(dist, valid_types)):
msg = '{0} `dist` ({1}) is not a valid type!'.format(
file_msg, type(dist))
raise LSSUtils_Error(msg)
# trans_opt
if not (trans_opt in list(range(1, 5))):
msg = '{0} `trans_opt` ({1}) is not within valid range!'.format(
file_msg, trans_opt)
raise LSSUtils_Error(msg)
##
## Centre's RA, DEC, DIST
# Right Ascension
if (isinstance(ra_cen, float)):
ra_cen = flip_angles(ra_cen)
else:
msg = '{0} `ra_cen` ({1}) is not a float!'.format(
file_msg, type(ra_cen))
raise LSSUtils_Error(msg)
# Declination
if (isinstance(dec_cen, float)):
dec_cen = flip_angles(dec_cen)
else:
msg = '{0} `dec_cen` ({1}) is not a float!'.format(
file_msg, type(dec_cen))
raise LSSUtils_Error(msg)
# Distance
if (isinstance(dist_cen, float)):
dist_cen = float(dist_cen)
else:
msg = '{0} `dist_cen` ({1}) is not a float!'.format(
file_msg, type(dist_cen))
raise LSSUtils_Error(msg)
##
## Check type of elements
# Right Ascension
if (isinstance(ra, float) or isinstance(ra, int)):
ra = np.array([ra])
else:
ra = np.array(ra)
# Declination
if (isinstance(dec, float) or isinstance(dec, int)):
dec = np.array([dec])
else:
dec = np.array(dec)
# Distance
if (isinstance(dist, float) or isinstance(dist, int)):
dist = np.array([dist])
else:
dist = np.array(dist)
##
## Converting to desired units
if unit == 'rad':
# Right Ascension
ra_rad = ra
ra_cen_deg = np.degrees(ra_cen)
ra_cen_rad = ra_cen
# Declination
dec_rad = dec
dec_cen_deg = np.degrees(dec_cen)
dec_cen_rad = dec_cen
elif unit == 'deg':
ra_rad = np.radians(ra)
ra_cen_deg = ra_cen
ra_cen_rad = np.radians(ra_cen)
# Declination
dec_rad = np.radians(dec)
dec_cen_deg = dec_cen
dec_cen_rad = np.radians(dec_cen)
##
## Initializing pandas DataFrame
dict_keys = ['ra', 'dec', 'dist']
coord_dict = dict(zip(dict_keys, np.vstack([ra, dec, dist])))
##
## Spherical to Cartesian transformation
## 1st tranformation
# Centre
x_cen = dist_cen * np.cos(ra_cen_rad) * np.cos(dec_cen_rad)
y_cen = dist_cen * np.sin(ra_cen_rad) * np.cos(dec_cen_rad)
z_cen = dist_cen * np.sin(dec_cen_rad)
# All galaxies
x = dist * np.cos(ra_rad) * np.cos(dec_rad)
y = dist * np.sin(ra_rad) * np.cos(dec_rad)
z = dist * np.sin(dec_rad)
##
## Rotations about z- and x-axis by `tau` and `Omega`
## Rotating the points, not the axes
x1 = x - x_cen
y1 = y - y_cen
z1 = z - z_cen
# Angles
omega = np.radians(90. - ra_cen_deg )
tau = np.radians(90. - dec_cen_deg)
# Rotations about z-axis by `omega`
x2 = x1 * np.cos(omega) - y1 * np.sin(omega)
y2 = x1 * np.sin(omega) + y1 * np.cos(omega)
z2 = z1.copy()
# Rotations about X-axis by `tau`
x3 = x2.copy()
y3 = y2 * np.cos(tau) - z2 * np.sin(tau)
z3 = z2 * np.sin(tau) + z2 * np.cos(tau)
##
## Definining which variables to return
# No Translation
if trans_opt == 1:
coord_dict['x'] = x
coord_dict['y'] = y
coord_dict['z'] = z
# Translation
if trans_opt == 2:
coord_dict['x'] = x1
coord_dict['y'] = y1
coord_dict['z'] = z1
# Translation + Rotation
if trans_opt == 3:
coord_dict['x'] = x2
coord_dict['y'] = y2
coord_dict['z'] = z2
# Translation + 2 Rotation (centered about the centre)
if trans_opt == 4:
coord_dict['x'] = x3
coord_dict['y'] = y3
coord_dict['z'] = z3
##
## Checking what object to return, i.e. DataFrame or python dictionary
if return_dict:
return coord_dict
else:
coord_pd = pd.DataFrame(coord_dict)
return coord_pd
## Coordinate Transformation - Cartesian to Spherical coordinates
def cart_to_sph_coords(dist, x_arr, y_arr, z_arr, return_dict=False,
unit='deg'):
"""
Transforms cartesian coordinates (x, y, z) into spherical
coordinates.
Parameters
-----------
dist : `np.ndarray`
Array of the the distance to the N-number of objects.
Shape (N,).
x_arr, y_arr, z_arr : `np.ndarray`
Arrays of the cartesian coordinates (x, y, z) of the
N objects. Shape (N,). Units are in `unit`.
return_dict : {True, False}, `bool`, optional
If `True`, this functions returns 2 dictionaries with `spherical`
and `cartesian` coordinates.
If `False`, it returns a `pandas.DataFrame` with the columns.
This variable is set to `False` by default.
unit : {'dec','rad'} `str`, optional
Unit of the output `ra`, `dec` coordinates.
This variable is set to `deg` by default.
Returns
-----------
coord_dict (coord_pd) : `dict`
Dictionary with spherical coordinates (`ra`, `dec`) and
distance (`dist`) of N elements. This value is returned if
`return_dict` is set to `True`. If not, a `pd.DataFrame` is
returned.
"""
file_msg = fd.Program_Msg(__file__)
## Check types of elements
# Units
unit_arr = ['deg', 'rad']
if not (unit in unit_arr):
'{0} `unit` ({1}) is not a valid input!'.format(
file_msg, unit)
# Valid types
valid_types = (float, int, np.ndarray, list)
# X-coordinate
if not (isinstance(x_arr, valid_types)):
msg = '{0} `x_arr` ({1}) is not a valid type!'.format(
file_msg, type(x_arr))
raise LSSUtils_Error(msg)
# Y-coordinate
if not (isinstance(y_arr, valid_types)):
msg = '{0} `y_arr` ({1}) is not a valid type!'.format(
file_msg, type(y_arr))
raise LSSUtils_Error(msg)
# Z-coordinate
if not (isinstance(z_arr, valid_types)):
msg = '{0} `z_arr` ({1}) is not a valid type!'.format(
file_msg, type(z_arr))
raise LSSUtils_Error(msg)
# Distance
if not (isinstance(dist, valid_types)):
msg = '{0} `dist` ({1}) is not a valid type!'.format(
file_msg, type(dist))
raise LSSUtils_Error(msg)
##
## Check type of elements
# X-Coordinate
if (isinstance(x_arr, float) or isinstance(x_arr, int)):
x_arr = np.array([x_arr])
else:
x_arr = np.array(x_arr)
# Y-Coordinate
if (isinstance(y_arr, float) or isinstance(y_arr, int)):
y_arr = np.array([y_arr])
else:
y_arr = np.array(y_arr)
# Z-Coordinate
if (isinstance(z_arr, float) or isinstance(z_arr, int)):
z_arr = np.array([z_arr])
else:
z_arr = np.array(z_arr)
# Distance
if (isinstance(dist, float) or isinstance(dist, int)):
dist = np.array([dist])
else:
dist = np.array(dist)
##
## Initializing pandas DataFrame
dict_keys = ['x', 'y', 'z', 'dist']
coord_dict = dict(zip(dict_keys, np.vstack([x_arr, y_arr, z_arr, dist])))
##
## Cartesian to Spherical coordinates
##
## - Normalized cartesian coordinates
cart_arr = np.column_stack([x_arr, y_arr, z_arr])
( x_val,
y_val,
z_val) = cart_arr/float(dist)
# Distance to object
dist = float(dist)
## Declination
dec_val = 90. - num.degrees(num.arccos(z_val))
## Right ascension
if x_val == 0:
if y_val > 0.:
ra_val = 90.
elif y_val < 0.:
ra_val = -90.
else:
ra_val = num.degrees(num.arctan(y_val/x_val))
##
## Seeing on which quadrant the point is at
if x_val < 0.:
ra_val += 180.
elif (x_val >= 0.) and (y_val < 0.):
ra_val += 360.
return ra_val, dec_val