Source code for spacepy.coordinates

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Implementation of Coords class functions for coordinate transformations

The coordinate systems supported by this module cover the most commonly
used geophysical and magnetospheric systems. The naming conventions can
follow the names used by the popular IRBEM library, but for inertial
systems we use a more consistent, fine-grained naming convention that
clarifies the different systems.

    Earth-centered Inertial Systems
    -------------------------------
    * **ECI2000** Earth-centered Inertial, J2000 epoch
    * **ECIMOD** Earth-centered Inertial, mean-of-date
    * **ECITOD** Earth-centered Inertial, true-of-date
    * **GEI** Geocentric Equatorial Inertial (IRBEM approximation of TOD)

    Magnetospheric Systems
    ----------------------
    * **GSM** Geocentric Solar Magnetospheric
    * **GSE** Geocentric Solar Ecliptic
    * **SM** Solar Magnetic
    * **MAG** Geomagnetic Coordinate System (aka CDMAG)

    Earth-fixed Systems
    -------------------
    * **GEO** Geocentric geographic, aka Earth-centered Earth-fixed
    * **GDZ** Geodetic coordinates

By convention *all* systems are treated as natively Cartesian except
geodetic (GDZ), which is defined in [altitude, latitude, longitude]
where altitude is relative to a reference ellipsoid. Similarly, distance
units are assumed to be Earth radii (Re) in all systems except GDZ, where
altitude is given in km. Conversions to GDZ will output altitude in km
regardless of the input distance units and conversions from GDZ will
output in Re regardless of input units. In all other cases, the distance
units will be preserved.

    .. versionchanged:: 0.3.0

    The new CTrans backend was added, which includes support for the names
    ``ECI2000``, ``ECIMOD``, ``ECITOD``, and ``CDMAG``. With the exception
    of ``ECIMOD``, these can be used with the existing IRBEM backend, and
    will be converted to their closest equivalents.

Notes on differences between representations
--------------------------------------------
IRBEM's coordinate transformations are low-accuracy and were written for
a library with a driving philosophy of speed and robustness as priorities.
The coordinate transformations are therefore approximate. Further, most of
the geophysical systems (e.g., GSE, SM) are derived from an inertial
system. It is standard practice to use ECIMOD as this system. However,
IRBEM does not currently make ECIMOD available as one of its inertial
systems. IRBEM's default inertial system (called GEI) is consistent with
an approximation of ECITOD. Hence there will be small differences between
IRBEM's transformations and those using SpacePy's CTrans backend.
Further sources of difference include: IRBEM uses a low-order approximation
to the sidereal time and other parameters; the calculation of the Earth-Sun
vector differs between the representations; the definitions of an Earth
radius differ (SpacePy = 6378.137km; IRBEM = 6371.2 km). SpacePy's in-built
representation is higher accuracy and is comprehensively tested, including
tests for consistency with other high accuracy packages such as LANLGeoMag
and AstroPy. However, for use cases where the required precision is of order
1 percent the output can be considered equivalent.

Setting options for coordinate transformation
---------------------------------------------
The backend for coordinate transformations can be provided at
instantiation of a :class:`Coords` object using a keyword
argument. However, for convenience and flexibility the options can be
set at the module level. Configurable options include the backend used
(:mod:`~spacepy.irbempy` or SpacePy's :mod:`~spacepy.ctrans`) and the
reference ellipsoid (only configurable for the SpacePy backend). A
warning will be raised if the backend is not set (either through the
defaults or the keyword argument). The final configurable option
(``itol``) is the maximum separation, in seconds, for which the
coordinate transformations will not be recalculated. To force all
transformations to use an exact transform for the time, set ``itol``
to zero. Values between 10s and 60s are recommended for speed while
also preserving accuracy, though different applications will require
different accuracies.  For example, assuming this module has been
imported as ``spc``, to set the SpacePy backend as the default and set
``itol`` to 5 seconds:

    >>> spc.DEFAULTS.set_values(use_irbem=False, itol=5)

Authors: Steven Morley and Josef Koller
Institution: Los ALamos National Laboratory
Contact: smorley@lanl.gov

Copyright 2010-2016 Los Alamos National Security, LLC.

"""

from collections import namedtuple
import numbers
import warnings
import numpy as np
from spacepy import help
import spacepy
import spacepy.datamodel as dm
import spacepy.time
from spacepy import ctrans

__contact__ = 'Steven Morley, smorley@lanl.gov'


# -----------------------------------------------
# space coordinate class
# -----------------------------------------------

SYSAXES_TYPES = {'GDZ': {'sph': 0, 'car': None},
                 'GEO': {'sph': None, 'car': 1}, 'GSM': {'sph': None, 'car': 2},
                 'GSE': {'sph': None, 'car': 3}, 'SM': {'sph': None, 'car': 4},
                 'GEI': {'sph': None, 'car': 5}, 'CDMAG': {'sph': None, 'car': 6},
                 'MAG': {'sph': None, 'car': 6}, 'SPH': {'sph': 7, 'car': None},
                 'RLL': {'sph': 8, 'car': None}, 'TOD': {'sph': None, 'car': 12},
                 'ECITOD': {'sph': None, 'car': 12}, 'J2000': {'sph': None, 'car': 13},
                 'ECI2000': {'sph': None, 'car': 13}, 'ECIMOD': {'sph': None, 'car': 14}}

SYS_EQUIV = {'GEI': 'ECITOD', 'TOD': 'ECITOD', 'J2000': 'ECI2000', 'MAG': 'CDMAG'}

IRBEM_RE = 6371.2  # kilometers


class __Defaults(object):
    '''Configuration factory for coordinates module default settings
    '''
    def __init__(self):
        self._fac = namedtuple('values', 'use_irbem, ellipsoid, itol')
        self.set_values()

    def __repr__(self):
        full = '{}'.format(self.values)
        b1 = full.index('(')
        return 'DEFAULTS{}'.format(full[b1:])

    def set_values(self, use_irbem=None, ellipsoid=ctrans.WGS84, itol=30):
        '''Set options for coordinate transforms

        Parameters
        ----------
        use_irbem : bool
            If True, use IRBEM as coordinate transform backend. If False, use SpacePy's
            implementations. Default is to use IRBEM (and warn) but this will change.
        ellipsoid : spacepy.ctrans.Ellipsoid
            A reference ellipsoid to use for the Earth radius definition and for geodetic
            coordinate conversion. SpacePy defaults to the WGS84 ellipsoid and the semi-major
            axis as the Earth radius (6738.137 km). The IRBEM backend uses WGS84, but expresses
            an Earth radius using 6371.2 km.
        '''
        self.values = self._fac(use_irbem=use_irbem,
                                ellipsoid=ellipsoid,
                                itol=itol)

DEFAULTS = __Defaults()


[docs]class Coords(object): ''' a = Coords( data, dtype, carsph, [units, ticks, use_irbem]) A class holding spatial coordinates and enabling transformation between coordinate systems. Coordinates can be stored as Cartesian or spherical and units are assumed to be Re (distance) and degrees (angle) .. note:: Although other units may be specified and will be carried through, most functions throughout SpacePy assume distances in Re and angles in degrees, regardless of specified units. By default, coordinate transforms are based on the IRBEM library; `its manual <http://svn.code.sf.net/p/irbem/code/trunk/manual/user_guide.html>`_ may prove useful. SpacePy also provides a framework for accurate coordinate transformations. This can be used by setting the `use_irbem` flag to False. In a future release of SpacePy this will become the default method. For a good reference on heliospheric and magnetospheric coordinate systems, see Franz & Harper, "Heliospheric Coordinate Systems", Planet. Space Sci., 50, pp 217-233, 2002 (https://doi.org/10.1016/S0032-0633(01)00119-2). Parameters ---------- data : list or ndarray, dim = (n,3) coordinate points [X,Y,Z] or [rad, lat, lon] dtype : string coordinate system; supported systems are defined in module-level documentation. Common systems include GEO, GSE, GSM, SM, MAG, ECIMOD carsph : string Cartesian or spherical, 'car' or 'sph' units : list of strings, optional standard are ['Re', 'Re', 'Re'] or ['Re', 'deg', 'deg'] depending on the carsph content. See note. ticks : Ticktock instance, optional used for coordinate transformations (see a.convert) use_irbem : bool .. versionadded:: 0.3.0 Set to True to use IRBEM for coordinate transforms (current default). Otherwise use SpacePy's coordinate transform library. Returns ------- out : Coords instance instance with a.data, a.carsph, etc. See Also -------- spacepy.coordinates.DEFAULTS spacepy.time.Ticktock Examples -------- >>> from spacepy import coordinates as coord >>> cvals = coord.Coords([[1,2,4],[1,2,2]], 'GEO', 'car') >>> cvals.x # returns all x coordinates array([1, 1]) >>> from spacepy.time import Ticktock >>> cvals.ticks = Ticktock(['2002-02-02T12:00:00', '2002-02-02T12:00:00'], 'ISO') # add ticks >>> newcoord = cvals.convert('GSM', 'sph') >>> newcoord .. currentmodule:: spacepy.coordinates .. autosummary:: ~Coords.append ~Coords.convert ~Coords.from_skycoord ~Coord.to_skycoord .. automethod:: append .. automethod:: convert .. automethod:: from_skycoord .. automethod:: to_skycoord ''' def __init__(self, data, dtype, carsph, units=None, ticks=None, use_irbem=None): if use_irbem is None: use_irbem = DEFAULTS.values.use_irbem if use_irbem is None: use_irbem = True warnings.warn('No coordinate backend specified; using IRBEM.' ' This default will change in the future.', DeprecationWarning) if use_irbem: from . import irbempy as op self.use_irbem = use_irbem # setup units self.Re = IRBEM_RE if use_irbem else DEFAULTS.values.ellipsoid['A'] # kilometers # Make sure that inputs are all formed correctly self.data = np.atleast_2d(data).astype(np.float_) if len(self.data.shape) != 2 or self.data.shape[-1] != 3: raise ValueError('Input position vectors must be Nx3') dtype = dtype.upper() carsph = carsph.lower() if dtype not in list(SYSAXES_TYPES.keys()): raise ValueError('This dtype=' + dtype + ' is not supported. Only ' + str(list(SYSAXES_TYPES.keys())) ) if carsph not in ['car', 'sph']: raise ValueError('This carsph=' + str(carsph) + ' is not supported. Only "car" or "sph"') if dtype in ['GDZ', 'RLL'] and carsph == 'car': raise ValueError("Geodetic coordinates are referenced to an ellipsoid " + "and can not be properly represented as Cartesian") # add ticks if ticks: assert len(ticks) == len(self.data), 'Ticktock dimensions seem off' self.ticks = ticks # GEO and SPH are cartesian/spherical versions of GEO, so use GEO for everything if dtype == 'SPH': dtype = 'GEO' # Make sure we're using the SpacePy name if we're using the SpacePy backend self.dtype = SYS_EQUIV[dtype] if (not use_irbem and dtype in SYS_EQUIV) else dtype self.sysaxes = SYSAXES_TYPES[dtype][carsph] self.carsph = carsph # setup units if units is None and carsph == 'car': # use standard units self.units = ['Re', 'Re', 'Re'] elif units is None and carsph == 'sph': self.units = ['Re', 'deg', 'deg'] else: self.units = units if dtype == 'GDZ' and carsph == 'sph': self.units = ['km', 'deg', 'deg'] # setup x,y,z etc if carsph == 'car': self.x = spacepy.datamodel.dmcopy(self.data[:, 0]) self.y = spacepy.datamodel.dmcopy(self.data[:, 1]) self.z = spacepy.datamodel.dmcopy(self.data[:, 2]) else: self.radi = spacepy.datamodel.dmcopy(self.data[:, 0]) self.lati = spacepy.datamodel.dmcopy(self.data[:, 1]) self.long = spacepy.datamodel.dmcopy(self.data[:, 2]) self.shape = np.shape(self.data) return None # ----------------------------------------------- def __str__(self): ''' a.__str__() or a Will be called when printing Coords instance a Returns ------- out : string string represenation of the instance Examples -------- >>> from spacepy.coordinates import Coords >>> y = Coords([[1,2,4],[1,2,2]], 'GEO', 'car') >>> y Coords( [[1 2 4] [1 2 2]] ), dtype=GEO,car, units=['Re', 'Re', 'Re'] ''' rstr = "Coords( {0} , '{1}', '{2}')".format(self.data.tolist(), self.dtype, self.carsph) return rstr __repr__ = __str__ # ----------------------------------------------- def __getitem__(self, idx): ''' a.__getitem__(idx) or a[idx] Will be called when requesting items in this instance Parameters ---------- idx : int integer numbers as index Returns ------- out : numpy array new values Examples -------- >>> from spacepy.coordinates import Coords >>> y = Coords([[1,2,4],[1,2,2]], 'GEO', 'car') >>> y[0] Coords( [[1 2 4]] ), dtype=GEO,car, units=['Re', 'Re', 'Re'] ''' arr = np.array(self.data) t_select = self.ticks[idx] if self.ticks else self.ticks return Coords(arr[idx].tolist(), self.dtype, self.carsph, self.units, t_select, use_irbem=self.use_irbem) # ----------------------------------------------- def __setitem__(self, idx, vals): ''' a.__setitem__(idx, vals) or a[idx] = vals Will be called setting items in this instance Parameters ---------- idx : int integer numbers as index vals : numpy array or list new values Examples -------- >>> from spacepy.coordinates import Coords >>> y = Coords([[1,2,4],[1,2,2]], 'GEO', 'car') >>> y[1] = [9,9,9] >>> y Coords( [[1 2 4] [9 9 9]] ), dtype=GEO,car, units=['Re', 'Re', 'Re'] ''' self.data[idx] = vals return # ----------------------------------------------- def __len__(self): ''' a.__len__() or len(a) Will be called when requesting the length, i.e. number of items Returns ------- out : int length Examples -------- >>> from spacepy.coordinates import Coords >>> y = Coords([[1,2,4],[1,2,2]], 'GEO', 'car') >>> len(y) 2 ''' if isinstance(self.data, (list, np.ndarray)): return len(self.data) else: return 1 return # -----------------------------------------------
[docs] def convert(self, returntype, returncarsph): '''Create a new Coords instance with new coordinate types Parameters ---------- returntype : string coordinate system, see module level documentation for supported systems returncarsph : string coordinate type, possible 'car' for Cartesian and 'sph' for spherical Returns ------- out : Coords object Coords object in the new coordinate system Examples -------- >>> from spacepy.coordinates import Coords >>> y = Coords([[1,2,4],[1,2,2]], 'GEO', 'car') >>> from spacepy.time import Ticktock >>> y.ticks = Ticktock(['2002-02-02T12:00:00', '2002-02-02T12:00:00'], 'ISO') >>> x = y.convert('SM','car') >>> x Coords( [[ 0.81134097 2.6493305 3.6500375 ] [ 0.92060408 2.30678864 1.68262126]] ), dtype=SM,car, units=['Re', 'Re', 'Re'] ''' # Check whether the name is an equivalent IRBEM name returnname = SYS_EQUIV[returntype] if returntype in SYS_EQUIV else returntype is_dtype_returntype = (self.dtype == returntype) or (self.dtype == returnname) if returntype in ['GDZ', 'RLL'] and returncarsph == 'car': raise ValueError("Geodetic coordinates are referenced to an ellipsoid " + "and can not be properly represented as Cartesian") # Then check whether a full coordinate transformation is required if is_dtype_returntype: if self.carsph == returncarsph: # no change necessary return self if self.carsph != returncarsph: # only car2sph or sph2car is needed if returncarsph == 'car': carsph = 'car' units = [self.units[0]]*3 data = sph2car(self.data) else: carsph = 'sph' units = [self.units[0], 'deg', 'deg'] data = car2sph(self.data) return Coords(data, self.dtype, carsph, units, self.ticks, use_irbem=self.use_irbem) # check the length of ticks and do the more complex conversions if self.ticks and (len(self.ticks) != len(self)): raise ValueError('Ticktock dimension does not match Coords dimensions') # check if car2sph is needed first for oneralib compatibility if (self.sysaxes is None): # a car2sph or sph2car is needed if self.carsph == 'sph': carsph = 'car' units = [self.units[0]]*3 data = sph2car(self.data) else: carsph = 'sph' units = [self.units[0], 'deg', 'deg'] data = car2sph(self.data) else: data = self.data.copy() units = dm.dmcopy(self.units) carsph = dm.dmcopy(self.carsph) if self.use_irbem: from . import irbempy as op # irbempy.coord_trans needs the passed Coords object to have the "from" dtype NewCoords = Coords(data, self.dtype, carsph, units, self.ticks, use_irbem=self.use_irbem) else: # SpacePy conversion takes input/output type separately, so set to desired output NewCoords = Coords(data, returntype if self.use_irbem else returnname, returncarsph, units, self.ticks, use_irbem=self.use_irbem) # now convert to other coordinate system if not is_dtype_returntype: assert self.ticks, "Time information required; add a .ticks attribute" if self.use_irbem: if returnname in ['GDZ', 'RLL'] and units[0] == 'km': NewCoords.data /= self.Re # IRBEM geodetic calculation requires Re input NewCoords.data = op.coord_trans(NewCoords, returntype, returncarsph) NewCoords.dtype = returntype if returnname == 'RLL' and units[0] == 'km': NewCoords.data[:, 0] *= self.Re # Preserve units else: if returnname in ['GDZ', 'RLL'] and self.units[0] == 'Re': # converting from a Cartesian referenced system to an # ellipsoid referenced system, input must be in km data *= self.Re # geodetic is defined in [km, deg, deg] NewCoords.data = np.atleast_2d(ctrans.convert_multitime(data, self.ticks, self.dtype, returnname, DEFAULTS.values)) if returnname == 'RLL' and units[0] == 'Re': # RLL is defined in Re, deg, deg, but all geodetic calcs are in km NewCoords.data[:, 0] = NewCoords.data[:, 0]/self.Re if (self.dtype in ['RLL', 'GDZ'] and returnname not in ['RLL', 'GDZ'] and units[0] == 'km'): # input was in km, needs to be Earth radii NewCoords.data /= self.Re if returncarsph == 'sph' and returnname not in ['GDZ', 'RLL']: NewCoords.data = car2sph(NewCoords.data) if self.dtype == 'GDZ': # GDZ is defined in km and everything else is strictly Re units[0] = 'Re' # Converting back from GDZ set units to Re NewCoords.carsph = returncarsph NewCoords.sysaxes = SYSAXES_TYPES[NewCoords.dtype][NewCoords.carsph] # Correct corresponding attributes and units if returncarsph == 'sph': if returnname == 'GDZ': NewCoords.units = ['km', 'deg', 'deg'] else: NewCoords.units = [units[0], 'deg', 'deg'] for k in ('x', 'y', 'z'): if hasattr(NewCoords, k): delattr(NewCoords, k) NewCoords.radi = NewCoords.data[:, 0] NewCoords.lati = NewCoords.data[:, 1] NewCoords.long = NewCoords.data[:, 2] else: # 'car' NewCoords.units = [units[0]]*3 for k in ('radi', 'lati', 'long'): if hasattr(NewCoords, k): delattr(NewCoords, k) NewCoords.x = NewCoords.data[:, 0] NewCoords.y = NewCoords.data[:, 1] NewCoords.z = NewCoords.data[:, 2] return NewCoords
# ----------------------------------------------- def __getstate__(self): ''' Is called when pickling See Also http://docs.python.org/library/pickle.html ''' odict = self.__dict__.copy() # copy the dict since we change it return odict def __setstate__(self, dict): ''' Is called when unpickling See Also http://docs.python.org/library/pickle.html ''' self.__dict__.update(dict) return # -----------------------------------------------
[docs] def append(self, other): '''Append another Coords instance to the current one Parameters ---------- other : Coords instance Coords instance to append Examples -------- ''' data = list(self.data) otherdata = other.convert(self.dtype, self.carsph) data.extend(list(otherdata.data)) with warnings.catch_warnings(): warnings.filterwarnings('ignore', message=r'Use of IRBEM to perform', category=DeprecationWarning, module=r'spacepy.coordinates$') newobj = Coords(data, dtype=self.dtype, carsph=self.carsph, use_irbem=self.use_irbem) return newobj
# -----------------------------------------------
[docs] def to_skycoord(self): ''' Create an Astropy SkyCoord instance based on this instance Returns ======= out : `astropy.coordinates.SkyCoord` This coordinate as an Astropy SkyCoord Notes ===== This method requires Astropy to be installed. This method uses the GEO coordinate frame as the common frame between the two libraries. ''' if self.ticks is None: raise ValueError("This method requires the attribute `ticks` to be specified.") try: import astropy.coordinates import astropy.time except ImportError as e: raise e.__class__("This method requires Astropy to be installed.") coord = self.convert('GEO', 'car') # GEO is Astropy ITRS data = coord.data * self.Re * 1000 # convert to meters obstime = coord.ticks.APT return astropy.coordinates.SkyCoord(x=data[:, 0], y=data[:, 1], z=data[:, 2], unit='m', frame='itrs', obstime=obstime)
# -----------------------------------------------
[docs] @classmethod def from_skycoord(cls, skycoord, use_irbem=None): ''' Create a Coords instance from an Astropy SkyCoord instance Parameters ========== skycoord : `astropy.coordinates.SkyCoord` The coordinate to be converted Returns ======= out : Coords instance The converted coordinate Notes ===== This method requires Astropy to be installed. This method uses the GEO coordinate frame as the common frame between the two libraries. ''' try: import astropy.coordinates except ImportError as e: raise e.__class__("This method requires Astropy to be installed.") skycoord = astropy.coordinates.SkyCoord(skycoord).itrs # Astropy ITRS is GEO # Ticktock constructor warns if use_irbem is None; just do right thing use_Re = IRBEM_RE if (use_irbem or use_irbem is None)\ else DEFAULTS.values.ellipsoid['A'] # kilometers data = (skycoord.cartesian.xyz.to('m').value / (1000 * use_Re)).T # convert Cartesian to Re units ticks = spacepy.time.Ticktock(skycoord.obstime) return cls(data, 'GEO', 'car', ticks=ticks, use_irbem=use_irbem)
[docs]def car2sph(car_in): """ Coordinate transformation from Cartesian to spherical Parameters ---------- - car_in (list or ndarray) : coordinate points in (n,3) shape with n coordinate points in units of [Re, Re, Re] = [x,y,z] Returns ------- - results (ndarray) : values after conversion to spherical coordinates in radius, latitude, longitude in units of [Re, deg, deg] Examples -------- >>> sph2car([1,45,0]) array([ 0.70710678, 0. , 0.70710678]) See Also -------- sph2car """ if isinstance(car_in[0], numbers.Number): car = np.array([car_in]) else: car = np.asanyarray(car_in) res = np.zeros(np.shape(car)) for i in np.arange(len(car)): x, y, z = car[i, 0], car[i, 1], car[i, 2] r = np.sqrt(x*x+y*y+z*z) sq = np.sqrt(x*x+y*y) if (x == 0) & (y == 0): # on the poles longi = 0. if z < 0: lati = -90. else: lati = 90.0 else: longi = np.arctan2(y, x)*180./np.pi lati = 90. - np.arctan2(sq, z)*180./np.pi res[i, :] = [r, lati, longi] if isinstance(car_in[0], numbers.Number): return res[0] else: return res
[docs]def sph2car(sph_in): """ Coordinate transformation from spherical to Cartesian Parameters ---------- - sph_in (list or ndarray) : coordinate points in (n,3) shape with n coordinate points in units of [Re, deg, deg] = [r, latitude, longitude] Returns ------- - results (ndarray) : values after conversion to cartesian coordinates x,y,z Examples -------- >>> sph2car([1,45,45]) array([ 0.5 , 0.5 , 0.70710678]) See Also -------- car2sph """ if isinstance(sph_in[0], numbers.Number): sph = np.array([sph_in]) else: sph = np.asanyarray(sph_in) res = np.zeros(np.shape(sph)) for i in np.arange(len(sph)): r, lati, longi = sph[i, 0], sph[i, 1], sph[i, 2] colat = np.pi/2. - lati*np.pi/180. x = r*np.sin(colat)*np.cos(longi*np.pi/180.) y = r*np.sin(colat)*np.sin(longi*np.pi/180.) z = r*np.cos(colat) res[i, :] = [x, y, z] if isinstance(sph_in[0], numbers.Number): return res[0] else: return res
[docs]def quaternionNormalize(Qin, scalarPos='last'): ''' Given an input quaternion (or array of quaternions), return the unit quaternion Parameters ---------- vec : array_like input quaternion to normalize Returns ------- out : array_like normalized quaternion Examples -------- >>> import spacepy.coordinates >>> spacepy.coordinates.quaternionNormalize([0.707, 0, 0.707, 0.2]) array([ 0.69337122, 0. , 0.69337122, 0.19614462]) ''' w = {'first': 0, 'last': 3}.get(scalarPos.lower()) if w is None: errmsg = 'quaternionNormalize: scalarPos must be set to "First" or "Last"' raise NotImplementedError(errmsg) Quse = np.asanyarray(Qin).astype(float) squeeze = Quse.ndim < 2 if squeeze: Quse = np.expand_dims(Quse, axis=0) magn = np.linalg.norm(Quse, axis=-1) # Find places where the magnitude is tiny, convert to unit real idx = np.where(magn <= 1e-12) magn[idx] = 1 outarr = Quse / np.expand_dims(magn, axis=magn.ndim) outarr[idx + (None,)] = 0 outarr[idx + (w,)] = 1 if squeeze: outarr = outarr[0, ...] return outarr
[docs]def quaternionRotateVector(Qin, Vin, scalarPos='last', normalize=True): ''' Given quaternions and vectors, return the vectors rotated by the quaternions Parameters ---------- Qin : array_like input quaternion to rotate by Vin : array-like input vector to rotate Returns ------- out : array_like rotated vector Examples -------- >>> import spacepy.coordinates >>> import numpy as np >>> vec = [1, 0, 0] >>> quat_wijk = [np.sin(np.pi/4), 0, np.sin(np.pi/4), 0.0] >>> quat_ijkw = [0.0, np.sin(np.pi/4), 0, np.sin(np.pi/4)] >>> spacepy.coordinates.quaternionRotateVector(quat_ijkw, vec) array([ 0., 0., -1.]) >>> spacepy.coordinates.quaternionRotateVector( ... quat_wijk, vec, scalarPos='first') array([ 0., 0., -1.]) See Also -------- quaternionMultiply ''' if scalarPos.lower() == 'last': i, j, k = 0, 1, 2 w = 3 elif scalarPos.lower() == 'first': i, j, k = 1, 2, 3 w = 0 else: errmsg = 'quaternionRotateVector: scalarPos must be set to "First" or "Last"' raise NotImplementedError(errmsg) Quse = np.asanyarray(Qin).astype(float) if normalize: Quse = quaternionNormalize(Quse, scalarPos=scalarPos) Vuse = np.asanyarray(Vin).astype(float) try: Quse.shape[1] except IndexError: Quse = np.asanyarray([Quse]) try: Vuse.shape[1] except IndexError: Vuse = np.asanyarray([Vuse]) try: assert Vuse.shape[0] == Quse.shape[0] except AssertionError: errmsg = ' '.join(['quaternionRotateVector:', 'Input vector array must have same length', 'as input quaternion array']) raise ValueError(errmsg) outarr = np.empty((Quse.shape[0], 3)) for idx, row in enumerate(Quse): Vrow = Vuse[idx] ii, jj, kk, ww = row[i]**2, row[j]**2, row[k]**2, row[w]**2 iw, jw, kw = row[i]*row[w], row[j]*row[w], row[k]*row[w] ij, ik, jk = row[i]*row[j], row[i]*row[k], row[j]*row[k] outarr[idx, 0] = ww*Vrow[0] + 2*jw*Vrow[2] - 2*kw*Vrow[1] + \ ii*Vrow[0] + 2*ij*Vrow[1] + 2*ik*Vrow[2] - \ kk*Vrow[0] - jj*Vrow[0] outarr[idx, 1] = 2*ij*Vrow[0] + jj*Vrow[1] + 2*jk*Vrow[2] + \ 2*kw*Vrow[0] - kk*Vrow[1] + ww*Vrow[1] - \ 2*iw*Vrow[2] - ii*Vrow[1] outarr[idx, 2] = 2*ik*Vrow[0] + 2*jk*Vrow[1] + kk*Vrow[2] - \ 2.0*jw*Vrow[0] - jj*Vrow[2] + 2.0*iw*Vrow[1] - \ ii*Vrow[2] + ww*Vrow[2] return outarr.squeeze()
[docs]def quaternionMultiply(Qin1, Qin2, scalarPos='last'): ''' Given quaternions, return the product, i.e. Qin1*Qin2 Parameters ---------- Qin1 : array_like input quaternion, first position Qin2 : array-like input quaternion, second position Returns ------- out : array_like quaternion product Examples -------- >>> import spacepy.coordinates >>> import numpy as np >>> vecX = [1, 0, 0] #shared X-axis >>> vecZ = [0, 0, 1] #unshared, but similar, Z-axis >>> quat_eci_to_gsm = [-0.05395384, 0.07589845, -0.15172533, 0.98402634] >>> quat_eci_to_gse = [ 0.20016056, 0.03445775, -0.16611386, 0.96496352] >>> quat_gsm_to_eci = spacepy.coordinates.quaternionConjugate( ... quat_eci_to_gsm) >>> quat_gse_to_gsm = spacepy.coordinates.quaternionMultiply( ... quat_gsm_to_eci, quat_eci_to_gse) >>> spacepy.coordinates.quaternionRotateVector(quat_gse_to_gsm, vecX) array([ 1.00000000e+00, 1.06536725e-09, -1.16892107e-08]) >>> spacepy.coordinates.quaternionRotateVector(quat_gse_to_gsm, vecZ) array([ 1.06802834e-08, -4.95669027e-01, 8.68511494e-01]) ''' if scalarPos.lower() == 'last': i, j, k = 0, 1, 2 w = 3 elif scalarPos.lower() == 'first': i, j, k = 1, 2, 3 w = 0 else: raise NotImplementedError('quaternionMultiply: scalarPos must be set to "First" or "Last"') Quse1 = np.asanyarray(Qin1).astype(float) Quse2 = np.asanyarray(Qin2).astype(float) try: Quse1.shape[1] except IndexError: Quse1 = np.asanyarray([Quse1]) try: Quse2.shape[1] except IndexError: Quse2 = np.asanyarray([Quse2]) try: assert Quse2.shape == Quse1.shape except AssertionError: raise ValueError('quaternionMultiply: Input quaternion arrays must have same length') outarr = np.empty_like(Quse1) for idx, row in enumerate(Quse1): row1 = row row2 = Quse2[idx] # vector components outarr[idx, i] = row1[w]*row2[i] + row1[i]*row2[w] + row1[j]*row2[k] - row1[k]*row2[j] outarr[idx, j] = row1[w]*row2[j] - row1[i]*row2[k] + row1[j]*row2[w] + row1[k]*row2[i] outarr[idx, k] = row1[w]*row2[k] + row1[i]*row2[j] - row1[j]*row2[i] + row1[k]*row2[w] # real part outarr[idx, w] = row1[w]*row2[w] - row1[i]*row2[i] - row1[j]*row2[j] - row1[k]*row2[k] return outarr.squeeze()
[docs]def quaternionConjugate(Qin, scalarPos='last'): ''' Given an input quaternion (or array of quaternions), return the conjugate Parameters ---------- Qin : array_like input quaternion to conjugate Returns ------- out : array_like conjugate quaternion Examples -------- >>> import spacepy.coordinates >>> spacepy.coordinates.quaternionConjugate( ... [0.707, 0, 0.707, 0.2], scalarPos='last') array([-0.707, -0. , -0.707, 0.2 ]) See Also -------- quaternionMultiply ''' if scalarPos.lower() == 'last': i, j, k = 0, 1, 2 w = 3 elif scalarPos.lower() == 'first': i, j, k = 1, 2, 3 w = 0 else: errmsg = 'quaternionConjugate: scalarPos must be set to "First" or "Last"' raise NotImplementedError(errmsg) Quse = np.asanyarray(Qin).astype(float) try: Quse.shape[1] except IndexError: Quse = np.asanyarray([Quse]) outarr = np.empty_like(Quse) for idx, row in enumerate(Quse): outarr[idx, i] = -row[i] outarr[idx, j] = -row[j] outarr[idx, k] = -row[k] outarr[idx, w] = row[w] return outarr.squeeze()
[docs]def quaternionFromMatrix(matrix, scalarPos='last'): ''' Given an input rotation matrix, return the equivalent quaternion The output has one fewer axis than the input (the last axis) and the shape is otherwise unchanged, allowing multi-dimensional matrix input. Parameters ---------- matrix : array_like input rotation matrix or array of matrices Returns ------- out : array_like Quaternions representing the same rotation as the input rotation matrices. Other Parameters ---------------- scalarPos : str Location of the scalar component of the output quaternion, either 'last' (default) or 'first'. Raises ------ NotImplementedError for invalid values of ``scalarPos`` ValueError for inputs which are obviously not valid 3D rotation matrices or arrays thereof: if the size doesn't end in (3, 3), if the matrix is not orthogonal, or not a proper rotation. See Also -------- quaternionToMatrix Notes ----- .. versionadded:: 0.2.2 No attempt is made to resolve the sign ambiguity; in particular, conversions of very similar matrices may result in equivalent quaternions with the opposite sign. This may have implications for interpolating a sequence of quaternions. The conversion of a rotation matrix to a quaternion suffers from some of the same disadvantages inherent to rotation matrices, such as potential numerical instabilities. Working in quaternion space as much as possible is recommended. There are several algorithms; the most well-known algorithm for this conversion is Shepperd's [#Shepperd]_, although the many "rediscoveries" indicate it is not sufficiently well-known. This function uses the method of Bar-Itzhack [#BarItzhack]_ (version 3), which should be resistant to small errors in the rotation matrix. As a result, the input checking is quite coarse and will likely accept many matrices that do not represent valid rotations. Also potentially of interest, although not implemented here, is Sarabandi and Thomas [#Sarabandi]_. References ---------- .. [#Shepperd] S.W. Shepperd, "Quaternion from rotation matrix," Journal of Guidance and Control, Vol. 1, No. 3, pp. 223-224, 1978, `doi:10.2514/3.55767b <https://doi.org/10.2514/3.55767b>`_ .. [#BarItzhack] I. Y. Bar-Itzhack, "New method for extracting the quaternion from a rotation matrix", AIAA Journal of Guidance, Control and Dynamics, 23 (6): 1085–1087, `doi:10.2514/2.4654 <https://doi.org/10.2514/2.4654>`_ .. [#Sarabandi] S. Sarabandi and F. Thomas, "Accurate Computation of Quaternions from Rotation Matrices", In: Lenarcic J., Parenti-Castelli V. (eds) Advances in Robot Kinematics 2018, Springer. `doi:10.1007/978-3-319-93188-3_5 <https://doi.org/10.1007/978-3-319-93188-3_5>`_ Examples -------- >>> import spacepy.coordinates >>> spacepy.coordinates.quaternionFromMatrix( ... [[ 0., 0., 1.], ... [ 1., 0., 0.], ... [ 0., 1., 0.]]) array([0.5, 0.5, 0.5, 0.5]) ''' if scalarPos.lower() not in ('last', 'first'): raise NotImplementedError( 'quaternionFromMatrix: scalarPos must be set to "First" or "Last"') matrix = np.asanyarray(matrix) if len(matrix.shape) < 2 or matrix.shape[-2:] != (3, 3): raise ValueError( 'Input does not appear to be 3D rotation matrix, wrong size.') for i in np.ndindex(matrix.shape[:-2]): m = matrix[i + (Ellipsis,)] if not np.allclose(np.dot(m, m.transpose()), np.identity(3), atol=0.2): raise ValueError('Input rotation matrix{} not orthogonal.'.format( '' if len(matrix.shape) == 2 else ' at ' + str(i))) det = np.linalg.det(m) if det < 0: raise ValueError('Input rotation matrix at {} not proper.' .format(str(i))) inshape = matrix.shape # Flatten out most dimensions, easier to index matrix = np.reshape(matrix, (-1, 3, 3)) # Indexing in Bar-Itzhack is reversed relative to numpy k = (matrix[..., 0, 0] - matrix[..., 1, 1] - matrix[..., 2, 2], matrix[..., 0, 1] + matrix[..., 1, 0], matrix[..., 0, 2] + matrix[..., 2, 0], matrix[..., 2, 1] - matrix[..., 1, 2], # row 0 matrix[..., 0, 1] + matrix[..., 1, 0], matrix[..., 1, 1] - matrix[..., 0, 0] - matrix[..., 2, 2], matrix[..., 1, 2] + matrix[..., 2, 1], matrix[..., 0, 2] - matrix[..., 2, 0], # row 1 matrix[..., 0, 2] + matrix[..., 2, 0], matrix[..., 1, 2] + matrix[..., 2, 1], matrix[..., 2, 2] - matrix[..., 0, 0] - matrix[..., 1, 1], matrix[..., 1, 0] - matrix[..., 0, 1], # row 2 matrix[..., 2, 1] - matrix[..., 1, 2], matrix[..., 0, 2] - matrix[..., 2, 0], matrix[..., 1, 0] - matrix[..., 0, 1], matrix[..., 0, 0] + matrix[..., 1, 1] + matrix[..., 2, 2] # row 3 ) # Stack together on a new axis after the input sample, then # split that axis into the 4x4 k arrays k = np.reshape(np.stack(k, axis=1), (-1, 4, 4)) k = k / 3. evals, evects = np.linalg.eig(k) evects = np.real(evects) idx = np.argmax(evals, axis=-1) Qout = evects[np.mgrid[0:len(idx)], :, idx] if scalarPos.lower() == 'first': Qout = np.roll(Qout, 1, axis=-1) # Recover original shape Qout = np.reshape(Qout, inshape[:-2] + (4,)) return Qout
[docs]def quaternionToMatrix(Qin, scalarPos='last', normalize=True): ''' Given an input quaternion, return the equivalent rotation matrix. The output has one more axis than the input (the last axis) and the shape is otherwise unchanged, allowing multi-dimensional quaternion input. Parameters ---------- Qin : array_like input quaternion or array of quaternions, must be normalized. Returns ------- out : array_like Rotation matrix Other Parameters ---------------- scalarPos : str Location of the scalar component of the input quaternion, either 'last' (default) or 'first'. normalize : True Normalize input quaternions before conversion (default). If False, raises error for non-normalized. Raises ------ NotImplementedError for invalid values of ``scalarPos``. ValueError for inputs which are not valid normalized quaternions or arrays thereof: if the size doesn't end in (4), if the quaternion is not normalized and ``normalize`` is False. See Also -------- quaternionFromMatrix Notes ----- .. versionadded:: 0.2.2 Implementation of the Euler–Rodrigues formula. Examples -------- >>> import spacepy.coordinates >>> spacepy.coordinates.quaternionToMatrix([0.5, 0.5, 0.5, 0.5]) array([[ 0., 0., 1.], [ 1., 0., 0.], [ 0., 1., 0.]]) ''' if scalarPos.lower() not in ('last', 'first'): raise NotImplementedError( 'quaternionToMatrix: scalarPos must be set to "First" or "Last"') Qin = np.asanyarray(Qin) if Qin.shape[-1] != 4: raise ValueError('Input does not appear to be quaternion, wrong size.') if normalize: Qin = quaternionNormalize(Qin, scalarPos=scalarPos) if scalarPos.lower() == 'first': Qin = np.roll(Qin, -1, axis=-1) if not np.allclose(np.sum(Qin ** 2, axis=-1), 1): raise ValueError('Input quaternion not normalized.') # Maintain dimensions at end for stacking a, b, c, d = Qin[..., -1:], Qin[..., 0:1], Qin[..., 1:2], Qin[..., 2:3] # Output array, "flattened" # This can probably be written with a clever vector notation... out = np.concatenate((a ** 2 + b ** 2 - c ** 2 - d ** 2, 2 * (b * c - a * d), 2 * (b * d + a * c), 2 * (b * c + a * d), a ** 2 + c ** 2 - b ** 2 - d ** 2, 2 * (c * d - a * b), 2 * (b * d - a * c), 2 * (c * d + a * b), a ** 2 + d ** 2 - b ** 2 - c ** 2, ), axis=-1) # And last two dims are 3x3 array for matrix out = np.reshape(out, out.shape[:-1] + (3, 3)) return out