aburrell / aacgmv2 / 491

pending completion

Build # 491

Build Type

Pull #57

travis-ci

Committed by

web-flow

Commit Message

<a href="https://github.com/aburrell/aacgmv2/commit/<a class=hub.com/aburrell/aacgmv2/commit/<a class="double-link" href="https://git"><a class=hub.com/aburrell/aacgmv2/commit/<a class="double-link" href="https://git"><a class=hub.com/aburrell/aacgmv2/commit/<a class="double-link" href="https://git"><a class=hub.com/aburrell/aacgmv2/commit/98f4d12ce209e9ae2c75995a9d14ab1b35ee5e57">98f4d12ce">&lt;a href=&quot;https://github.com/aburrell/aacgmv2/commit/</a><a class="double-link" href="https://github.com/aburrell/aacgmv2/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/aburrell/aacgmv2/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/aburrell/aacgmv2/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/aburrell/aacgmv2/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/aburrell/aacgmv2/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/aburrell/aacgmv2/commit/98f4d12ce209e9ae2c75995a9d14ab1b35ee5e57">98f4d12ce</a><a href="https://github.com/aburrell/aacgmv2/commit/98f4d12ce209e9ae2c75995a9d14ab1b35ee5e57">&quot;&gt;&amp;lt;a href=&amp;quot;https://github.com/aburrell/aacgmv2/commit/&lt;/a&gt;&lt;a class=&quot;double-link&quot; href=&quot;https://github.com/aburrell/aacgmv2/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://git&quot;&gt;&amp;lt;a class=&lt;/a&gt;hub.com/aburrell/aacgmv2/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://git&quot;&gt;&amp;lt;a class=&lt;/a&gt;hub.com/aburrell/aacgmv2/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://git&quot;&gt;&amp;lt;a class=&lt;/a&gt;hub.com/aburrell/aacgmv2/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://git&quot;&gt;&amp;lt;a class=&lt;/a&gt;hub.com/aburrell/aacgmv2/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://git&quot;&gt;&amp;lt;a class=&lt;/a&gt;hub.com/aburrell/aacgmv2/commit/98f4d12ce209e9ae2c75995a9d14ab1b35ee5e57&quot;&gt;98f4d12ce&lt;/a&gt;&lt;a href=&quot;https://github.com/aburrell/aacgmv2/commit/98f4d12ce209e9ae2c75995a9d14ab1b35ee5e57&quot;&gt... (continued)

Pull Request Pull Request #57: Un-deprecation of utilities

Run Details

232 of 254 branches covered (91.34%)

Branch coverage included in aggregate %.

163 of 163 new or added lines in 8 files covered. (100.0%)

1156 of 1186 relevant lines covered (97.47%)

2.88 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

96.42

/aacgmv2/wrapper.py

# Copyright (C) 2019 NRL
# Author: Angeline Burrell
# Disclaimer: This code is under the MIT license, whose details can be found at
# the root in the LICENSE file
#
# -*- coding: utf-8 -*-
"""Pythonic wrappers for AACGM-V2 C functions.

"""

from __future__ import division, absolute_import, unicode_literals
import datetime as dt
import numpy as np
import os
import sys
import warnings

import aacgmv2
import aacgmv2._aacgmv2 as c_aacgmv2
from aacgmv2._aacgmv2 import TRACE, ALLOWTRACE, BADIDEA

if sys.version_info.major == 2:
    warnings.filterwarnings('error')


def test_time(dtime):
    """ Test the time input and ensure it is a dt.datetime object

    Parameters
    ----------
    dtime : (unknown)
        Time input in an untested format

    Returns
    -------
    dtime : (dt.datetime)
        Time as a datetime object

    Raises
    ------
    ValueError if time is not a dt.date or dt.datetime object

    """
    if isinstance(dtime, dt.date):
        # Because datetime objects identify as both dt.date and dt.datetime,
        # you need an extra test here to ensure you don't lose the time
        # attributes
        if not isinstance(dtime, dt.datetime):
            dtime = dt.datetime.combine(dtime, dt.time(0))
    elif not isinstance(dtime, dt.datetime):
        raise ValueError('time variable (dtime) must be a datetime object')

    return dtime


def test_height(height, bit_code):
    """ Test the input height and ensure it is appropriate for the method

    Parameters
    ----------
    height : (float)
        Height to test in km
    bit_code : (int)
        Code string denoting method to use

    Returns
    -------
    good_height : (boolean)
        True if height and method are appropriate, False if not

    Notes
    -----
    Appropriate altitude ranges for the different methods are explored in
    Shepherd (2014).  Summarized, they are:
    Coefficients: 0-2000 km
    Tracing: 0-1 Earth Radius

    Altitudes below zero will work, but will not provide a good representation
    of the magnetic field because it goes beyond the intended scope of these
    coordiantes.

    If you use the 'BADIDEA' code, you can bypass all constraints, but it
    is a Bad Idea!  If you include a high enough altiutde, the code may hang.

    """
    # Test for heights that are allowed but not within the intended scope
    # of the coordinate system.  The routine will work, but the user should
    # be aware that the results are not as reliable
    if height < 0:
        aacgmv2.logger.warning('conversion not intended for altitudes < 0 km')

    # Test the conditions for using the coefficient method
    if(height > aacgmv2.high_alt_coeff
       and not (bit_code & (TRACE | ALLOWTRACE | BADIDEA))):
        estr = ''.join(['coefficients are not valid for altitudes above ',
                        '{:.0f} km. You '.format(aacgmv2.high_alt_coeff),
                        'must either use field-line tracing (trace=True or',
                        ' allowtrace=True) or indicate you know this is a',
                        ' bad idea'])
        aacgmv2.logger.error(estr)
        return False

    # Test the conditions for using the tracing method
    if height > aacgmv2.high_alt_trace and not (bit_code & BADIDEA):
        estr = ''.join(['these coordinates are not intended for the ',
                        'magnetosphere! You must indicate that you know ',
                        'this is a bad idea.  If you continue, it is ',
                        'possible that the code will hang.'])
        aacgmv2.logger.error(estr)
        return False

    return True


def set_coeff_path(igrf_file=False, coeff_prefix=False):
    """Sets the IGRF_COEFF and AACGMV_V2_DAT_PREFIX environment variables.

    Parameters
    ----------
    igrf_file : (str or bool)
        Full filename of IGRF coefficient file, True to use
        aacgmv2.IGRF_COEFFS, or False to leave as is. (default=False)
    coeff_prefix : (str or bool)
        Location and file prefix for aacgm coefficient files, True to use
        aacgmv2.AACGM_V2_DAT_PREFIX, or False to leave as is. (default=False)

    """

    # Define coefficient file prefix if requested
    if coeff_prefix is not False:
        # Use the default value, if one was not supplied (allow None to
        # comply with depricated behaviour)
        if coeff_prefix is True or coeff_prefix is None:
            coeff_prefix = aacgmv2.AACGM_v2_DAT_PREFIX

        if hasattr(os, "unsetenv"):
            os.unsetenv('AACGM_v2_DAT_PREFIX')
        else:
            del os.environ['AACGM_v2_DAT_PREFIX']
        os.environ['AACGM_v2_DAT_PREFIX'] = coeff_prefix

    # Define IGRF file if requested
    if igrf_file is not False:
        # Use the default value, if one was not supplied (allow None to
        # comply with depricated behaviour)
        if igrf_file is True or igrf_file is None:
            igrf_file = aacgmv2.IGRF_COEFFS

        if hasattr(os, "unsetenv"):
            os.unsetenv('IGRF_COEFFS')
        else:
            del os.environ['IGRF_COEFFS']
        os.environ['IGRF_COEFFS'] = igrf_file

    return


def convert_latlon(in_lat, in_lon, height, dtime, method_code="G2A"):
    """Converts between geomagnetic coordinates and AACGM coordinates

    Parameters
    ----------
    in_lat : (float)
        Input latitude in degrees N (code specifies type of latitude)
    in_lon : (float)
        Input longitude in degrees E (code specifies type of longitude)
    height : (float)
        Altitude above the surface of the earth in km
    dtime : (datetime)
        Datetime for magnetic field
    method_code : (str or int)
        Bit code or string denoting which type(s) of conversion to perform
        G2A        - geographic (geodetic) to AACGM-v2
        A2G        - AACGM-v2 to geographic (geodetic)
        TRACE      - use field-line tracing, not coefficients
        ALLOWTRACE - use trace only above 2000 km
        BADIDEA    - use coefficients above 2000 km
        GEOCENTRIC - assume inputs are geocentric w/ RE=6371.2
        (default is "G2A")

    Returns
    -------
    out_lat : (float)
        Output latitude in degrees N
    out_lon : (float)
        Output longitude in degrees E
    out_r : (float)
        Geocentric radial distance (R_Earth) or altitude above the surface of
        the Earth (km)

    Raises
    ------
    ValueError if input is incorrect
    RuntimeError if unable to set AACGMV2 datetime

    """
    # Test time
    dtime = test_time(dtime)

    # Initialise output
    lat_out = np.nan
    lon_out = np.nan
    r_out = np.nan

    # Set the coordinate coversion method code in bits
    try:
        bit_code = convert_str_to_bit(method_code.upper())
    except AttributeError:
        bit_code = method_code

    if not isinstance(bit_code, int):
        raise ValueError("unknown method code {:}".format(method_code))

    # Test height that may or may not cause failure
    if not test_height(height, bit_code):
        return lat_out, lon_out, r_out

    # Test latitude range
    if abs(in_lat) > 90.0:
        # Allow latitudes with a small deviation from the maximum
        # (+/- 90 degrees) to be set to 90
        if abs(in_lat) > 90.1:
            raise ValueError('unrealistic latitude')
        in_lat = np.sign(in_lat) * 90.0

    # Constrain longitudes between -180 and 180
    in_lon = ((in_lon + 180.0) % 360.0) - 180.0

    # Set current date and time
    try:
        c_aacgmv2.set_datetime(dtime.year, dtime.month, dtime.day, dtime.hour,
                               dtime.minute, dtime.second)
    except (TypeError, RuntimeError) as err:
        raise RuntimeError("cannot set time for {:}: {:}".format(dtime, err))

    # convert location
    try:
        lat_out, lon_out, r_out = c_aacgmv2.convert(in_lat, in_lon, height,
                                                    bit_code)
    except Exception:
        err = sys.exc_info()[0]
        estr = "unable to perform conversion at {:.1f},".format(in_lat)
        estr = "{:s}{:.1f} {:.1f} km, {:} ".format(estr, in_lon, height, dtime)
        estr = "{:s}using method {:}: {:}".format(estr, bit_code, err)
        aacgmv2.logger.warning(estr)
        pass

    return lat_out, lon_out, r_out


def convert_latlon_arr(in_lat, in_lon, height, dtime, method_code="G2A"):
    """Converts between geomagnetic coordinates and AACGM coordinates.

    Parameters
    ----------
    in_lat : (np.ndarray or list or float)
        Input latitude in degrees N (method_code specifies type of latitude)
    in_lon : (np.ndarray or list or float)
        Input longitude in degrees E (method_code specifies type of longitude)
    height : (np.ndarray or list or float)
        Altitude above the surface of the earth in km
    dtime : (datetime)
        Single datetime object for magnetic field
    method_code : (int or str)
        Bit code or string denoting which type(s) of conversion to perform
        G2A        - geographic (geodetic) to AACGM-v2
        A2G        - AACGM-v2 to geographic (geodetic)
        TRACE      - use field-line tracing, not coefficients
        ALLOWTRACE - use trace only above 2000 km
        BADIDEA    - use coefficients above 2000 km
        GEOCENTRIC - assume inputs are geocentric w/ RE=6371.2
        (default = "G2A")

    Returns
    -------
    out_lat : (np.ndarray)
        Output latitudes in degrees N
    out_lon : (np.ndarray)
        Output longitudes in degrees E
    out_r : (np.ndarray)
        Geocentric radial distance (R_Earth) or altitude above the surface of
        the Earth (km)

    Raises
    ------
    ValueError if input is incorrect
    RuntimeError if unable to set AACGMV2 datetime

    Notes
    -----
    At least one of in_lat, in_lon, and height must be a list or array.

    If errors are encountered, NaN or Inf will be included in the input so
    that all successful calculations are returned.  To select only good values
    use a function like `np.isfinite`.

    Multi-dimensional arrays are not allowed.

    """
    # Recast the data as numpy arrays
    in_lat = np.array(in_lat)
    in_lon = np.array(in_lon)
    height = np.array(height)

    # If one or two of these elements is a float, int, or single element array,
    # create an array equal to the length of the longest input
    test_array = np.array([len(in_lat.shape), len(in_lon.shape),
                           len(height.shape)])

    if test_array.max() > 1:
        raise ValueError("unable to process multi-dimensional arrays")
    else:
        if test_array.max() == 0:
            aacgmv2.logger.info("".join(["for a single location, consider ",
                                         "using convert_latlon or ",
                                         "get_aacgm_coord"]))
            in_lat = np.array([in_lat])
            in_lon = np.array([in_lon])
            height = np.array([height])
        else:
            max_len = max([len(arr) for i, arr in enumerate([in_lat, in_lon,
                                                             height])
                           if test_array[i] > 0])

            if not test_array[0] or (len(in_lat) == 1 and max_len > 1):
                in_lat = np.full(shape=(max_len,), fill_value=in_lat)
            if not test_array[1] or (len(in_lon) == 1 and max_len > 1):
                in_lon = np.full(shape=(max_len,), fill_value=in_lon)
            if not test_array[2] or (len(height) == 1 and max_len > 1):
                height = np.full(shape=(max_len,), fill_value=height)

    # Ensure that lat, lon, and height are the same length or if the lengths
    # differ that the different ones contain only a single value
    if not (in_lat.shape == in_lon.shape and in_lat.shape == height.shape):
        raise ValueError('lat, lon, and height arrays are mismatched')

    # Test time
    dtime = test_time(dtime)

    # Initialise output
    lat_out = np.full(shape=in_lat.shape, fill_value=np.nan)
    lon_out = np.full(shape=in_lon.shape, fill_value=np.nan)
    r_out = np.full(shape=height.shape, fill_value=np.nan)

    # Test and set the conversion method code
    try:
        bit_code = convert_str_to_bit(method_code.upper())
    except AttributeError:
        bit_code = method_code

    if not isinstance(bit_code, int):
        raise ValueError("unknown method code {:}".format(method_code))

    # Test height
    if not test_height(np.nanmax(height), bit_code):
        return lat_out, lon_out, r_out

    # Test latitude range
    if np.abs(in_lat).max() > 90.0:
        if np.abs(in_lat).max() > 90.1:
            raise ValueError('unrealistic latitude')
        in_lat = np.clip(in_lat, -90.0, 90.0)

    # Constrain longitudes between -180 and 180
    in_lon = ((in_lon + 180.0) % 360.0) - 180.0

    # Set current date and time
    try:
        c_aacgmv2.set_datetime(dtime.year, dtime.month, dtime.day, dtime.hour,
                               dtime.minute, dtime.second)
    except (TypeError, RuntimeError) as err:
        raise RuntimeError("cannot set time for {:}: {:}".format(dtime, err))

    try:
        lat_out, lon_out, r_out, bad_ind = c_aacgmv2.convert_arr(list(in_lat),
                                                                 list(in_lon),
                                                                 list(height),
                                                                 bit_code)

        # Cast the output as numpy arrays or masks
        lat_out = np.array(lat_out)
        lon_out = np.array(lon_out)
        r_out = np.array(r_out)
        bad_ind = np.array(bad_ind) >= 0

        # Replace any bad indices with NaN, casting output as numpy arrays
        if np.any(bad_ind):
            lat_out[bad_ind] = np.nan
            lon_out[bad_ind] = np.nan
            r_out[bad_ind] = np.nan
    except SystemError as serr:
        aacgmv2.logger.warning('C Error encountered: {:}'.format(serr))

    return lat_out, lon_out, r_out


def get_aacgm_coord(glat, glon, height, dtime, method="ALLOWTRACE"):
    """Get AACGM latitude, longitude, and magnetic local time

    Parameters
    ----------
    glat : (float)
        Geodetic latitude in degrees N
    glon : (float)
        Geodetic longitude in degrees E
    height : (float)
        Altitude above the surface of the earth in km
    dtime : (datetime)
        Date and time to calculate magnetic location
    method : (str)
        String denoting which type(s) of conversion to perform
        TRACE      - use field-line tracing, not coefficients
        ALLOWTRACE - use trace only above 2000 km
        BADIDEA    - use coefficients above 2000 km
        GEOCENTRIC - assume inputs are geocentric w/ RE=6371.2
        (default = "TRACE")

    Returns
    -------
    mlat : (float)
        magnetic latitude in degrees N
    mlon : (float)
        magnetic longitude in degrees E
    mlt : (float)
        magnetic local time in hours

    """
    # Initialize method code
    method_code = "G2A|{:s}".format(method)

    # Get magnetic lat and lon.
    mlat, mlon, _ = convert_latlon(glat, glon, height, dtime,
                                   method_code=method_code)

    # Get magnetic local time (output is always an array, so extract value)
    mlt = np.nan if np.isnan(mlon) else convert_mlt(mlon, dtime, m2a=False)[0]

    return mlat, mlon, mlt


def get_aacgm_coord_arr(glat, glon, height, dtime, method="ALLOWTRACE"):
    """Get AACGM latitude, longitude, and magnetic local time

    Parameters
    ----------
    glat : (np.array or list)
        Geodetic latitude in degrees N
    glon : (np.array or list)
        Geodetic longitude in degrees E
    height : (np.array or list)
        Altitude above the surface of the earth in km
    dtime : (datetime)
        Date and time to calculate magnetic location
    method : (str)
        String denoting which type(s) of conversion to perform
        TRACE      - use field-line tracing, not coefficients
        ALLOWTRACE - use trace only above 2000 km
        BADIDEA    - use coefficients above 2000 km
        GEOCENTRIC - assume inputs are geocentric w/ RE=6371.2
        (default = "TRACE")
        (default = "TRACE")

    Returns
    -------
    mlat : (float)
        magnetic latitude in degrees N
    mlon : (float)
        magnetic longitude in degrees E
    mlt : (float)
        magnetic local time in hours

    """
    # Initialize method code
    method_code = "G2A|{:s}".format(method)

    # Get magnetic lat and lon.
    mlat, mlon, _ = convert_latlon_arr(glat, glon, height, dtime,
                                       method_code=method_code)

    if np.any(np.isfinite(mlon)):
        # Get magnetic local time
        mlt = convert_mlt(mlon, dtime, m2a=False)
    else:
        mlt = np.full(shape=len(mlat), fill_value=np.nan)

    return mlat, mlon, mlt


def convert_str_to_bit(method_code):
    """convert string code specification to bit code specification

    Parameters
    ----------
    method_code : (str)
        Bitwise code for passing options into converter (default=0)
        G2A        - geographic (geodetic) to AACGM-v2
        A2G        - AACGM-v2 to geographic (geodetic)
        TRACE      - use field-line tracing, not coefficients
        ALLOWTRACE - use trace only above 2000 km
        BADIDEA    - use coefficients above 2000 km
        GEOCENTRIC - assume inputs are geocentric w/ RE=6371.2

    Returns
    -------
    bit_code : (int)
        Method code specification in bits

    Notes
    -----
    Multiple codes should be seperated by pipes '|'.  Invalid parts of the code
    are ignored and no code defaults to 'G2A'.

    """

    convert_code = {"G2A": c_aacgmv2.G2A, "A2G": c_aacgmv2.A2G,
                    "TRACE": c_aacgmv2.TRACE, "BADIDEA": c_aacgmv2.BADIDEA,
                    "GEOCENTRIC": c_aacgmv2.GEOCENTRIC,
                    "ALLOWTRACE": c_aacgmv2.ALLOWTRACE}

    # Force upper case, remove any spaces, and split along pipes
    method_codes = method_code.upper().replace(" ", "").split("|")

    # Add the valid parts of the code, invalid elements are ignored
    bit_code = sum([convert_code[k] for k in method_codes
                    if k in convert_code.keys()])

    return bit_code


def convert_bool_to_bit(a2g=False, trace=False, allowtrace=False,
                        badidea=False, geocentric=False):
    """convert boolian flags to bit code specification

    Parameters
    ----------
    a2g : (bool)
        True for AACGM-v2 to geographic (geodetic), False otherwise
        (default=False)
    trace : (bool)
        If True, use field-line tracing, not coefficients (default=False)
    allowtrace : (bool)
        If True, use trace only above 2000 km (default=False)
    badidea : (bool)
        If True, use coefficients above 2000 km (default=False)
    geocentric : (bool)
        True for geodetic, False for geocentric w/RE=6371.2 (default=False)

    Returns
    -------
    bit_code : (int)
        code specification in bits

    """

    bit_code = c_aacgmv2.A2G if a2g else c_aacgmv2.G2A

    if trace:
        bit_code += c_aacgmv2.TRACE
    if allowtrace:
        bit_code += c_aacgmv2.ALLOWTRACE
    if badidea:
        bit_code += c_aacgmv2.BADIDEA
    if geocentric:
        bit_code += c_aacgmv2.GEOCENTRIC

    return bit_code


def convert_mlt(arr, dtime, m2a=False):
    """Converts between magnetic local time (MLT) and AACGM-v2 longitude

    Parameters
    ----------
    arr : (array-like or float)
        Magnetic longitudes (degrees E) or MLTs (hours) to convert
    dtime : (array-like or datetime.datetime)
        Date and time for MLT conversion in Universal Time (UT).
    m2a : (bool)
        Convert MLT to AACGM-v2 longitude (True) or magnetic longitude to MLT
        (False).  (default=False)

    Returns
    -------
    out : (np.ndarray)
        Converted coordinates/MLT in degrees E or hours (as appropriate)

    Notes
    -----
    This routine previously based on Laundal et al. 2016, but now uses the
    improved calculation available in AACGM-V2.4.

    """

    arr = np.asarray(arr)
    if arr.shape == ():
        arr = np.array([arr])

    if len(arr.shape) > 1:
        raise ValueError("unable to process multi-dimensional arrays")

    # Test time
    try:
        dtime = test_time(dtime)
        years = [dtime.year for dd in arr]
        months = [dtime.month for dd in arr]
        days = [dtime.day for dd in arr]
        hours = [dtime.hour for dd in arr]
        minutes = [dtime.minute for dd in arr]
        seconds = [dtime.second for dd in arr]
    except ValueError as verr:
        dtime = np.asarray(dtime)
        if dtime.shape == ():
            raise ValueError(verr)
        elif dtime.shape != arr.shape:
            raise ValueError("array input for datetime and MLon/MLT must match")

        years = [dd.year for dd in dtime]
        months = [dd.month for dd in dtime]
        days = [dd.day for dd in dtime]
        hours = [dd.hour for dd in dtime]
        minutes = [dd.minute for dd in dtime]
        seconds = [dd.second for dd in dtime]

    arr = list(arr)

    # Calculate desired location, C routines set date and time
    if m2a:
        # Get the magnetic longitude
        if len(arr) == 1:
            out = c_aacgmv2.inv_mlt_convert(years[0], months[0], days[0],
                                            hours[0], minutes[0], seconds[0],
                                            arr[0])
        else:
            out = c_aacgmv2.inv_mlt_convert_arr(years, months, days, hours,
                                                minutes, seconds, arr)
    else:
        # Get magnetic local time
        if len(arr) == 1:
            out = c_aacgmv2.mlt_convert(years[0], months[0], days[0], hours[0],
                                        minutes[0], seconds[0], arr[0])
            out = np.array([out])
        else:
            out = np.array(c_aacgmv2.mlt_convert_arr(years, months, days, hours,
                                                     minutes, seconds, arr))

    return out

1	# Copyright (C) 2019 NRL
2	# Author: Angeline Burrell
3	# Disclaimer: This code is under the MIT license, whose details can be found at
4	# the root in the LICENSE file
5	#
6	# -- coding: utf-8 --
7	"""Pythonic wrappers for AACGM-V2 C functions.	3✔
8
9	"""
10
11	from __future__ import division, absolute_import, unicode_literals	3✔
12	import datetime as dt	3✔
13	import numpy as np	3✔
14	import os	3✔
15	import sys	3✔
16	import warnings	3✔
17
18	import aacgmv2	3✔
19	import aacgmv2._aacgmv2 as c_aacgmv2	3✔
20	from aacgmv2._aacgmv2 import TRACE, ALLOWTRACE, BADIDEA	3✔
21
22	if sys.version_info.major == 2:	3✔
23	warnings.filterwarnings('error')	1✔
24
25
26	def test_time(dtime):	3✔
27	""" Test the time input and ensure it is a dt.datetime object
28
29	Parameters
30	----------
31	dtime : (unknown)
32	Time input in an untested format
33
34	Returns
35	-------
36	dtime : (dt.datetime)
37	Time as a datetime object
38
39	Raises
40	------
41	ValueError if time is not a dt.date or dt.datetime object
42
43	"""
44	if isinstance(dtime, dt.date):	3✔
45	# Because datetime objects identify as both dt.date and dt.datetime,
46	# you need an extra test here to ensure you don't lose the time
47	# attributes
48	if not isinstance(dtime, dt.datetime):	3✔
49	dtime = dt.datetime.combine(dtime, dt.time(0))	3✔
50	elif not isinstance(dtime, dt.datetime):	3!
51	raise ValueError('time variable (dtime) must be a datetime object')	3✔
52
53	return dtime	3✔
54
55
56	def test_height(height, bit_code):	3✔
57	""" Test the input height and ensure it is appropriate for the method
58
59	Parameters
60	----------
61	height : (float)
62	Height to test in km
63	bit_code : (int)
64	Code string denoting method to use
65
66	Returns
67	-------
68	good_height : (boolean)
69	True if height and method are appropriate, False if not
70
71	Notes
72	-----
73	Appropriate altitude ranges for the different methods are explored in
74	Shepherd (2014). Summarized, they are:
75	Coefficients: 0-2000 km
76	Tracing: 0-1 Earth Radius
77
78	Altitudes below zero will work, but will not provide a good representation
79	of the magnetic field because it goes beyond the intended scope of these
80	coordiantes.
81
82	If you use the 'BADIDEA' code, you can bypass all constraints, but it
83	is a Bad Idea! If you include a high enough altiutde, the code may hang.
84
85	"""
86	# Test for heights that are allowed but not within the intended scope
87	# of the coordinate system. The routine will work, but the user should
88	# be aware that the results are not as reliable
89	if height < 0:	3✔
90	aacgmv2.logger.warning('conversion not intended for altitudes < 0 km')	3✔
91
92	# Test the conditions for using the coefficient method
93	if(height > aacgmv2.high_alt_coeff	3✔
94	and not (bit_code & (TRACE \| ALLOWTRACE \| BADIDEA))):
95	estr = ''.join(['coefficients are not valid for altitudes above ',	3✔
96	'{:.0f} km. You '.format(aacgmv2.high_alt_coeff),
97	'must either use field-line tracing (trace=True or',
98	' allowtrace=True) or indicate you know this is a',
99	' bad idea'])
100	aacgmv2.logger.error(estr)	3✔
101	return False	3✔
102
103	# Test the conditions for using the tracing method
104	if height > aacgmv2.high_alt_trace and not (bit_code & BADIDEA):	3✔
105	estr = ''.join(['these coordinates are not intended for the ',	3✔
106	'magnetosphere! You must indicate that you know ',
107	'this is a bad idea. If you continue, it is ',
108	'possible that the code will hang.'])
109	aacgmv2.logger.error(estr)	3✔
110	return False	3✔
111
112	return True	3✔
113
114
115	def set_coeff_path(igrf_file=False, coeff_prefix=False):	3✔
116	"""Sets the IGRF_COEFF and AACGMV_V2_DAT_PREFIX environment variables.
117
118	Parameters
119	----------
120	igrf_file : (str or bool)
121	Full filename of IGRF coefficient file, True to use
122	aacgmv2.IGRF_COEFFS, or False to leave as is. (default=False)
123	coeff_prefix : (str or bool)
124	Location and file prefix for aacgm coefficient files, True to use
125	aacgmv2.AACGM_V2_DAT_PREFIX, or False to leave as is. (default=False)
126
127	"""
128
129	# Define coefficient file prefix if requested
130	if coeff_prefix is not False:	3✔
131	# Use the default value, if one was not supplied (allow None to
132	# comply with depricated behaviour)
133	if coeff_prefix is True or coeff_prefix is None:	3✔
134	coeff_prefix = aacgmv2.AACGM_v2_DAT_PREFIX	3✔
135
136	if hasattr(os, "unsetenv"):	3!
137	os.unsetenv('AACGM_v2_DAT_PREFIX')	3✔
138	else:
139	del os.environ['AACGM_v2_DAT_PREFIX']	×
140	os.environ['AACGM_v2_DAT_PREFIX'] = coeff_prefix	3✔
141
142	# Define IGRF file if requested
143	if igrf_file is not False:	3✔
144	# Use the default value, if one was not supplied (allow None to
145	# comply with depricated behaviour)
146	if igrf_file is True or igrf_file is None:	3✔
147	igrf_file = aacgmv2.IGRF_COEFFS	3✔
148
149	if hasattr(os, "unsetenv"):	3!
150	os.unsetenv('IGRF_COEFFS')	3✔
151	else:
152	del os.environ['IGRF_COEFFS']	×
153	os.environ['IGRF_COEFFS'] = igrf_file	3✔
154
155	return	3✔
156
157
158	def convert_latlon(in_lat, in_lon, height, dtime, method_code="G2A"):	3✔
159	"""Converts between geomagnetic coordinates and AACGM coordinates
160
161	Parameters
162	----------
163	in_lat : (float)
164	Input latitude in degrees N (code specifies type of latitude)
165	in_lon : (float)
166	Input longitude in degrees E (code specifies type of longitude)
167	height : (float)
168	Altitude above the surface of the earth in km
169	dtime : (datetime)
170	Datetime for magnetic field
171	method_code : (str or int)
172	Bit code or string denoting which type(s) of conversion to perform
173	G2A - geographic (geodetic) to AACGM-v2
174	A2G - AACGM-v2 to geographic (geodetic)
175	TRACE - use field-line tracing, not coefficients
176	ALLOWTRACE - use trace only above 2000 km
177	BADIDEA - use coefficients above 2000 km
178	GEOCENTRIC - assume inputs are geocentric w/ RE=6371.2
179	(default is "G2A")
180
181	Returns
182	-------
183	out_lat : (float)
184	Output latitude in degrees N
185	out_lon : (float)
186	Output longitude in degrees E
187	out_r : (float)
188	Geocentric radial distance (R_Earth) or altitude above the surface of
189	the Earth (km)
190
191	Raises
192	------
193	ValueError if input is incorrect
194	RuntimeError if unable to set AACGMV2 datetime
195
196	"""
197	# Test time
198	dtime = test_time(dtime)	3✔
199
200	# Initialise output
201	lat_out = np.nan	3✔
202	lon_out = np.nan	3✔
203	r_out = np.nan	3✔
204
205	# Set the coordinate coversion method code in bits
206	try:	3✔
207	bit_code = convert_str_to_bit(method_code.upper())	3✔
208	except AttributeError:	3✔
209	bit_code = method_code	3✔
210
211	if not isinstance(bit_code, int):	3✔
212	raise ValueError("unknown method code {:}".format(method_code))	3✔
213
214	# Test height that may or may not cause failure
215	if not test_height(height, bit_code):	3✔
216	return lat_out, lon_out, r_out	3✔
217
218	# Test latitude range
219	if abs(in_lat) > 90.0:	3✔
220	# Allow latitudes with a small deviation from the maximum
221	# (+/- 90 degrees) to be set to 90
222	if abs(in_lat) > 90.1:	3✔
223	raise ValueError('unrealistic latitude')	3✔
224	in_lat = np.sign(in_lat) * 90.0	3✔
225
226	# Constrain longitudes between -180 and 180
227	in_lon = ((in_lon + 180.0) % 360.0) - 180.0	3✔
228
229	# Set current date and time
230	try:	3✔
231	c_aacgmv2.set_datetime(dtime.year, dtime.month, dtime.day, dtime.hour,	3✔
232	dtime.minute, dtime.second)
233	except (TypeError, RuntimeError) as err:	×
234	raise RuntimeError("cannot set time for {:}: {:}".format(dtime, err))	×
235
236	# convert location
237	try:	3✔
238	lat_out, lon_out, r_out = c_aacgmv2.convert(in_lat, in_lon, height,	3✔
239	bit_code)
240	except Exception:	2✔
241	err = sys.exc_info()[0]	2✔
242	estr = "unable to perform conversion at {:.1f},".format(in_lat)	2✔
243	estr = "{:s}{:.1f} {:.1f} km, {:} ".format(estr, in_lon, height, dtime)	2✔
244	estr = "{:s}using method {:}: {:}".format(estr, bit_code, err)	2✔
245	aacgmv2.logger.warning(estr)	2✔
246	pass	2✔
247
248	return lat_out, lon_out, r_out	3✔
249
250
251	def convert_latlon_arr(in_lat, in_lon, height, dtime, method_code="G2A"):	3✔
252	"""Converts between geomagnetic coordinates and AACGM coordinates.
253
254	Parameters
255	----------
256	in_lat : (np.ndarray or list or float)
257	Input latitude in degrees N (method_code specifies type of latitude)
258	in_lon : (np.ndarray or list or float)
259	Input longitude in degrees E (method_code specifies type of longitude)
260	height : (np.ndarray or list or float)
261	Altitude above the surface of the earth in km
262	dtime : (datetime)
263	Single datetime object for magnetic field
264	method_code : (int or str)
265	Bit code or string denoting which type(s) of conversion to perform
266	G2A - geographic (geodetic) to AACGM-v2
267	A2G - AACGM-v2 to geographic (geodetic)
268	TRACE - use field-line tracing, not coefficients
269	ALLOWTRACE - use trace only above 2000 km
270	BADIDEA - use coefficients above 2000 km
271	GEOCENTRIC - assume inputs are geocentric w/ RE=6371.2
272	(default = "G2A")
273
274	Returns
275	-------
276	out_lat : (np.ndarray)
277	Output latitudes in degrees N
278	out_lon : (np.ndarray)
279	Output longitudes in degrees E
280	out_r : (np.ndarray)
281	Geocentric radial distance (R_Earth) or altitude above the surface of
282	the Earth (km)
283
284	Raises
285	------
286	ValueError if input is incorrect
287	RuntimeError if unable to set AACGMV2 datetime
288
289	Notes
290	-----
291	At least one of in_lat, in_lon, and height must be a list or array.
292
293	If errors are encountered, NaN or Inf will be included in the input so
294	that all successful calculations are returned. To select only good values
295	use a function like `np.isfinite`.
296
297	Multi-dimensional arrays are not allowed.
298
299	"""
300	# Recast the data as numpy arrays
301	in_lat = np.array(in_lat)	3✔
302	in_lon = np.array(in_lon)	3✔
303	height = np.array(height)	3✔
304
305	# If one or two of these elements is a float, int, or single element array,
306	# create an array equal to the length of the longest input
307	test_array = np.array([len(in_lat.shape), len(in_lon.shape),	3✔
308	len(height.shape)])
309
310	if test_array.max() > 1:	3✔
311	raise ValueError("unable to process multi-dimensional arrays")	3✔
312	else:
313	if test_array.max() == 0:	3✔
314	aacgmv2.logger.info("".join(["for a single location, consider ",	3✔
315	"using convert_latlon or ",
316	"get_aacgm_coord"]))
317	in_lat = np.array([in_lat])	3✔
318	in_lon = np.array([in_lon])	3✔
319	height = np.array([height])	3✔
320	else:
321	max_len = max([len(arr) for i, arr in enumerate([in_lat, in_lon,	3✔
322	height])
323	if test_array[i] > 0])
324
325	if not test_array[0] or (len(in_lat) == 1 and max_len > 1):	3✔
326	in_lat = np.full(shape=(max_len,), fill_value=in_lat)	3✔
327	if not test_array[1] or (len(in_lon) == 1 and max_len > 1):	3✔
328	in_lon = np.full(shape=(max_len,), fill_value=in_lon)	3✔
329	if not test_array[2] or (len(height) == 1 and max_len > 1):	3✔
330	height = np.full(shape=(max_len,), fill_value=height)	3✔
331
332	# Ensure that lat, lon, and height are the same length or if the lengths
333	# differ that the different ones contain only a single value
334	if not (in_lat.shape == in_lon.shape and in_lat.shape == height.shape):	3✔
335	raise ValueError('lat, lon, and height arrays are mismatched')	3✔
336
337	# Test time
338	dtime = test_time(dtime)	3✔
339
340	# Initialise output
341	lat_out = np.full(shape=in_lat.shape, fill_value=np.nan)	3✔
342	lon_out = np.full(shape=in_lon.shape, fill_value=np.nan)	3✔
343	r_out = np.full(shape=height.shape, fill_value=np.nan)	3✔
344
345	# Test and set the conversion method code
346	try:	3✔
347	bit_code = convert_str_to_bit(method_code.upper())	3✔
348	except AttributeError:	3✔
349	bit_code = method_code	3✔
350
351	if not isinstance(bit_code, int):	3✔
352	raise ValueError("unknown method code {:}".format(method_code))	3✔
353
354	# Test height
355	if not test_height(np.nanmax(height), bit_code):	3✔
356	return lat_out, lon_out, r_out	3✔
357
358	# Test latitude range
359	if np.abs(in_lat).max() > 90.0:	3✔
360	if np.abs(in_lat).max() > 90.1:	3✔
361	raise ValueError('unrealistic latitude')	3✔
362	in_lat = np.clip(in_lat, -90.0, 90.0)	3✔
363
364	# Constrain longitudes between -180 and 180
365	in_lon = ((in_lon + 180.0) % 360.0) - 180.0	3✔
366
367	# Set current date and time
368	try:	3✔
369	c_aacgmv2.set_datetime(dtime.year, dtime.month, dtime.day, dtime.hour,	3✔
370	dtime.minute, dtime.second)
371	except (TypeError, RuntimeError) as err:	×
372	raise RuntimeError("cannot set time for {:}: {:}".format(dtime, err))	×
373
374	try:	3✔
375	lat_out, lon_out, r_out, bad_ind = c_aacgmv2.convert_arr(list(in_lat),	3✔
376	list(in_lon),
377	list(height),
378	bit_code)
379
380	# Cast the output as numpy arrays or masks
381	lat_out = np.array(lat_out)	3✔
382	lon_out = np.array(lon_out)	3✔
383	r_out = np.array(r_out)	3✔
384	bad_ind = np.array(bad_ind) >= 0	3✔
385
386	# Replace any bad indices with NaN, casting output as numpy arrays
387	if np.any(bad_ind):	3✔
388	lat_out[bad_ind] = np.nan	2✔
389	lon_out[bad_ind] = np.nan	2✔
390	r_out[bad_ind] = np.nan	2✔
391	except SystemError as serr:	×
392	aacgmv2.logger.warning('C Error encountered: {:}'.format(serr))	×
393
394	return lat_out, lon_out, r_out	3✔
395
396
397	def get_aacgm_coord(glat, glon, height, dtime, method="ALLOWTRACE"):	3✔
398	"""Get AACGM latitude, longitude, and magnetic local time
399
400	Parameters
401	----------
402	glat : (float)
403	Geodetic latitude in degrees N
404	glon : (float)
405	Geodetic longitude in degrees E
406	height : (float)
407	Altitude above the surface of the earth in km
408	dtime : (datetime)
409	Date and time to calculate magnetic location
410	method : (str)
411	String denoting which type(s) of conversion to perform
412	TRACE - use field-line tracing, not coefficients
413	ALLOWTRACE - use trace only above 2000 km
414	BADIDEA - use coefficients above 2000 km
415	GEOCENTRIC - assume inputs are geocentric w/ RE=6371.2
416	(default = "TRACE")
417
418	Returns
419	-------
420	mlat : (float)
421	magnetic latitude in degrees N
422	mlon : (float)
423	magnetic longitude in degrees E
424	mlt : (float)
425	magnetic local time in hours
426
427	"""
428	# Initialize method code
429	method_code = "G2A\|{:s}".format(method)	3✔
430
431	# Get magnetic lat and lon.
432	mlat, mlon, _ = convert_latlon(glat, glon, height, dtime,	3✔
433	method_code=method_code)
434
435	# Get magnetic local time (output is always an array, so extract value)
436	mlt = np.nan if np.isnan(mlon) else convert_mlt(mlon, dtime, m2a=False)[0]	3✔
437
438	return mlat, mlon, mlt	3✔
439
440
441	def get_aacgm_coord_arr(glat, glon, height, dtime, method="ALLOWTRACE"):	3✔
442	"""Get AACGM latitude, longitude, and magnetic local time
443
444	Parameters
445	----------
446	glat : (np.array or list)
447	Geodetic latitude in degrees N
448	glon : (np.array or list)
449	Geodetic longitude in degrees E
450	height : (np.array or list)
451	Altitude above the surface of the earth in km
452	dtime : (datetime)
453	Date and time to calculate magnetic location
454	method : (str)
455	String denoting which type(s) of conversion to perform
456	TRACE - use field-line tracing, not coefficients
457	ALLOWTRACE - use trace only above 2000 km
458	BADIDEA - use coefficients above 2000 km
459	GEOCENTRIC - assume inputs are geocentric w/ RE=6371.2
460	(default = "TRACE")
461	(default = "TRACE")
462
463	Returns
464	-------
465	mlat : (float)
466	magnetic latitude in degrees N
467	mlon : (float)
468	magnetic longitude in degrees E
469	mlt : (float)
470	magnetic local time in hours
471
472	"""
473	# Initialize method code
474	method_code = "G2A\|{:s}".format(method)	3✔
475
476	# Get magnetic lat and lon.
477	mlat, mlon, _ = convert_latlon_arr(glat, glon, height, dtime,	3✔
478	method_code=method_code)
479
480	if np.any(np.isfinite(mlon)):	3✔
481	# Get magnetic local time
482	mlt = convert_mlt(mlon, dtime, m2a=False)	3✔
483	else:
484	mlt = np.full(shape=len(mlat), fill_value=np.nan)	3✔
485
486	return mlat, mlon, mlt	3✔
487
488
489	def convert_str_to_bit(method_code):	3✔
490	"""convert string code specification to bit code specification
491
492	Parameters
493	----------
494	method_code : (str)
495	Bitwise code for passing options into converter (default=0)
496	G2A - geographic (geodetic) to AACGM-v2
497	A2G - AACGM-v2 to geographic (geodetic)
498	TRACE - use field-line tracing, not coefficients
499	ALLOWTRACE - use trace only above 2000 km
500	BADIDEA - use coefficients above 2000 km
501	GEOCENTRIC - assume inputs are geocentric w/ RE=6371.2
502
503	Returns
504	-------
505	bit_code : (int)
506	Method code specification in bits
507
508	Notes
509	-----
510	Multiple codes should be seperated by pipes '\|'. Invalid parts of the code
511	are ignored and no code defaults to 'G2A'.
512
513	"""
514
515	convert_code = {"G2A": c_aacgmv2.G2A, "A2G": c_aacgmv2.A2G,	3✔
516	"TRACE": c_aacgmv2.TRACE, "BADIDEA": c_aacgmv2.BADIDEA,
517	"GEOCENTRIC": c_aacgmv2.GEOCENTRIC,
518	"ALLOWTRACE": c_aacgmv2.ALLOWTRACE}
519
520	# Force upper case, remove any spaces, and split along pipes
521	method_codes = method_code.upper().replace(" ", "").split("\|")	3✔
522
523	# Add the valid parts of the code, invalid elements are ignored
524	bit_code = sum([convert_code[k] for k in method_codes	3✔
525	if k in convert_code.keys()])
526
527	return bit_code	3✔
528
529
530	def convert_bool_to_bit(a2g=False, trace=False, allowtrace=False,	3✔
531	badidea=False, geocentric=False):
532	"""convert boolian flags to bit code specification
533
534	Parameters
535	----------
536	a2g : (bool)
537	True for AACGM-v2 to geographic (geodetic), False otherwise
538	(default=False)
539	trace : (bool)
540	If True, use field-line tracing, not coefficients (default=False)
541	allowtrace : (bool)
542	If True, use trace only above 2000 km (default=False)
543	badidea : (bool)
544	If True, use coefficients above 2000 km (default=False)
545	geocentric : (bool)
546	True for geodetic, False for geocentric w/RE=6371.2 (default=False)
547
548	Returns
549	-------
550	bit_code : (int)
551	code specification in bits
552
553	"""
554
555	bit_code = c_aacgmv2.A2G if a2g else c_aacgmv2.G2A	3✔
556
557	if trace:	3✔
558	bit_code += c_aacgmv2.TRACE	3✔
559	if allowtrace:	3✔
560	bit_code += c_aacgmv2.ALLOWTRACE	3✔
561	if badidea:	3✔
562	bit_code += c_aacgmv2.BADIDEA	3✔
563	if geocentric:	3✔
564	bit_code += c_aacgmv2.GEOCENTRIC	3✔
565
566	return bit_code	3✔
567
568
569	def convert_mlt(arr, dtime, m2a=False):	3✔
570	"""Converts between magnetic local time (MLT) and AACGM-v2 longitude
571
572	Parameters
573	----------
574	arr : (array-like or float)
575	Magnetic longitudes (degrees E) or MLTs (hours) to convert
576	dtime : (array-like or datetime.datetime)
577	Date and time for MLT conversion in Universal Time (UT).
578	m2a : (bool)
579	Convert MLT to AACGM-v2 longitude (True) or magnetic longitude to MLT
580	(False). (default=False)
581
582	Returns
583	-------
584	out : (np.ndarray)
585	Converted coordinates/MLT in degrees E or hours (as appropriate)
586
587	Notes
588	-----
589	This routine previously based on Laundal et al. 2016, but now uses the
590	improved calculation available in AACGM-V2.4.
591
592	"""
593
594	arr = np.asarray(arr)	3✔
595	if arr.shape == ():	3✔
596	arr = np.array([arr])	3✔
597
598	if len(arr.shape) > 1:	3✔
599	raise ValueError("unable to process multi-dimensional arrays")	3✔
600
601	# Test time
602	try:	3✔
603	dtime = test_time(dtime)	3✔
604	years = [dtime.year for dd in arr]	3✔
605	months = [dtime.month for dd in arr]	3✔
606	days = [dtime.day for dd in arr]	3✔
607	hours = [dtime.hour for dd in arr]	3✔
608	minutes = [dtime.minute for dd in arr]	3✔
609	seconds = [dtime.second for dd in arr]	3✔
610	except ValueError as verr:	3✔
611	dtime = np.asarray(dtime)	3✔
612	if dtime.shape == ():	3✔
613	raise ValueError(verr)	3✔
614	elif dtime.shape != arr.shape:	3✔
615	raise ValueError("array input for datetime and MLon/MLT must match")	3✔
616
617	years = [dd.year for dd in dtime]	3✔
618	months = [dd.month for dd in dtime]	3✔
619	days = [dd.day for dd in dtime]	3✔
620	hours = [dd.hour for dd in dtime]	3✔
621	minutes = [dd.minute for dd in dtime]	3✔
622	seconds = [dd.second for dd in dtime]	3✔
623
624	arr = list(arr)	3✔
625
626	# Calculate desired location, C routines set date and time
627	if m2a:	3✔
628	# Get the magnetic longitude
629	if len(arr) == 1:	3✔
630	out = c_aacgmv2.inv_mlt_convert(years[0], months[0], days[0],	3✔
631	hours[0], minutes[0], seconds[0],
632	arr[0])
633	else:
634	out = c_aacgmv2.inv_mlt_convert_arr(years, months, days, hours,	3✔
635	minutes, seconds, arr)
636	else:
637	# Get magnetic local time
638	if len(arr) == 1:	3✔
639	out = c_aacgmv2.mlt_convert(years[0], months[0], days[0], hours[0],	3✔
640	minutes[0], seconds[0], arr[0])
641	out = np.array([out])	3✔
642	else:
643	out = np.array(c_aacgmv2.mlt_convert_arr(years, months, days, hours,	3✔
644	minutes, seconds, arr))
645
646	return out	3✔

aburrell / aacgmv2 / 491

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous