aacgmv2.wrapper.convert_latlon() - Code Metrics - Inspection of "Un-deprecation of utilities" - aburrell/aacgmv2 - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Pull Request — develop (#57)

by Angeline

created 2020-09-03 15:32 UTC

aacgmv2.wrapper.convert_latlon() B

↳ Parent: aacgmv2.wrapper

Complexity

Conditions

Size

Total Lines	91
Code Lines	34

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
eloc	34
dl	0
loc	91
rs	7.1973
c	0
b	0
f	0
cc	8
nop	5

How to fix Long Method

# 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.
8
9			"""
10
11			from __future__ import division, absolute_import, unicode_literals
12			import datetime as dt
13			import numpy as np
14			import os
15			import sys
16			import warnings
17
18			import aacgmv2
19			import aacgmv2._aacgmv2 as c_aacgmv2
20			from aacgmv2._aacgmv2 import TRACE, ALLOWTRACE, BADIDEA
21
22			if sys.version_info.major == 2:
23			warnings.filterwarnings('error')
24
25
26			def test_time(dtime):
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):
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):
49			dtime = dt.datetime.combine(dtime, dt.time(0))
50			elif not isinstance(dtime, dt.datetime):
51			raise ValueError('time variable (dtime) must be a datetime object')
52
53			return dtime
54
55
56			def test_height(height, bit_code):
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:
90			aacgmv2.logger.warning('conversion not intended for altitudes < 0 km')
91
92			# Test the conditions for using the coefficient method
93			if(height > aacgmv2.high_alt_coeff
94			and not (bit_code & (TRACE \| ALLOWTRACE \| BADIDEA))):
95			estr = ''.join(['coefficients are not valid for altitudes above ',
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)
101			return False
102
103			# Test the conditions for using the tracing method
104			if height > aacgmv2.high_alt_trace and not (bit_code & BADIDEA):
105			estr = ''.join(['these coordinates are not intended for the ',
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)
110			return False
111
112			return True
113
114
115			def set_coeff_path(igrf_file=False, coeff_prefix=False):
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:
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:
134			coeff_prefix = aacgmv2.AACGM_v2_DAT_PREFIX
135
136			if hasattr(os, "unsetenv"):
137			os.unsetenv('AACGM_v2_DAT_PREFIX')
138			else:
139			del os.environ['AACGM_v2_DAT_PREFIX']
140			os.environ['AACGM_v2_DAT_PREFIX'] = coeff_prefix
141
142			# Define IGRF file if requested
143			if igrf_file is not False:
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:
147			igrf_file = aacgmv2.IGRF_COEFFS
148
149			if hasattr(os, "unsetenv"):
150			os.unsetenv('IGRF_COEFFS')
151			else:
152			del os.environ['IGRF_COEFFS']
153			os.environ['IGRF_COEFFS'] = igrf_file
154
155			return
156
157
158			def convert_latlon(in_lat, in_lon, height, dtime, method_code="G2A"):
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)
199
200			# Initialise output
201			lat_out = np.nan
202			lon_out = np.nan
203			r_out = np.nan
204
205			# Set the coordinate coversion method code in bits
206			try:
207			bit_code = convert_str_to_bit(method_code.upper())
208			except AttributeError:
209			bit_code = method_code
210
211			if not isinstance(bit_code, int):
212			raise ValueError("unknown method code {:}".format(method_code))
213
214			# Test height that may or may not cause failure
215			if not test_height(height, bit_code):
216			return lat_out, lon_out, r_out
217
218			# Test latitude range
219			if abs(in_lat) > 90.0:
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:
223			raise ValueError('unrealistic latitude')
224			in_lat = np.sign(in_lat) * 90.0
225
226			# Constrain longitudes between -180 and 180
227			in_lon = ((in_lon + 180.0) % 360.0) - 180.0
228
229			# Set current date and time
230			try:
231			c_aacgmv2.set_datetime(dtime.year, dtime.month, dtime.day, dtime.hour,
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:
238			lat_out, lon_out, r_out = c_aacgmv2.convert(in_lat, in_lon, height,
239			bit_code)
240			except Exception:
241			err = sys.exc_info()[0]
242			estr = "unable to perform conversion at {:.1f},".format(in_lat)
243			estr = "{:s}{:.1f} {:.1f} km, {:} ".format(estr, in_lon, height, dtime)
244			estr = "{:s}using method {:}: {:}".format(estr, bit_code, err)
245			aacgmv2.logger.warning(estr)
246			pass
247
248			return lat_out, lon_out, r_out
249
250
251			def convert_latlon_arr(in_lat, in_lon, height, dtime, method_code="G2A"):
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)
302			in_lon = np.array(in_lon)
303			height = np.array(height)
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),
308			len(height.shape)])
309
310			if test_array.max() > 1:
311			raise ValueError("unable to process multi-dimensional arrays")
312			else:
313			if test_array.max() == 0:
314			aacgmv2.logger.info("".join(["for a single location, consider ",
315			"using convert_latlon or ",
316			"get_aacgm_coord"]))
317			in_lat = np.array([in_lat])
318			in_lon = np.array([in_lon])
319			height = np.array([height])
320			else:
321			max_len = max([len(arr) for i, arr in enumerate([in_lat, in_lon,
322			height])
323			if test_array[i] > 0])
324
325			if not test_array[0] or (len(in_lat) == 1 and max_len > 1):
326			in_lat = np.full(shape=(max_len,), fill_value=in_lat)
327			if not test_array[1] or (len(in_lon) == 1 and max_len > 1):
328			in_lon = np.full(shape=(max_len,), fill_value=in_lon)
329			if not test_array[2] or (len(height) == 1 and max_len > 1):
330			height = np.full(shape=(max_len,), fill_value=height)
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):
335			raise ValueError('lat, lon, and height arrays are mismatched')
336
337			# Test time
338			dtime = test_time(dtime)
339
340			# Initialise output
341			lat_out = np.full(shape=in_lat.shape, fill_value=np.nan)
342			lon_out = np.full(shape=in_lon.shape, fill_value=np.nan)
343			r_out = np.full(shape=height.shape, fill_value=np.nan)
344
345			# Test and set the conversion method code
346			try:
347			bit_code = convert_str_to_bit(method_code.upper())
348			except AttributeError:
349			bit_code = method_code
350
351			if not isinstance(bit_code, int):
352			raise ValueError("unknown method code {:}".format(method_code))
353
354			# Test height
355			if not test_height(np.nanmax(height), bit_code):
356			return lat_out, lon_out, r_out
357
358			# Test latitude range
359			if np.abs(in_lat).max() > 90.0:
360			if np.abs(in_lat).max() > 90.1:
361			raise ValueError('unrealistic latitude')
362			in_lat = np.clip(in_lat, -90.0, 90.0)
363
364			# Constrain longitudes between -180 and 180
365			in_lon = ((in_lon + 180.0) % 360.0) - 180.0
366
367			# Set current date and time
368			try:
369			c_aacgmv2.set_datetime(dtime.year, dtime.month, dtime.day, dtime.hour,
370			dtime.minute, dtime.second)
371			except (TypeError, RuntimeError) as err:
372			raise RuntimeError("cannot set time for {:}: {:}".format(dtime, err))
373
374			try:
375			lat_out, lon_out, r_out, bad_ind = c_aacgmv2.convert_arr(list(in_lat),
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)
382			lon_out = np.array(lon_out)
383			r_out = np.array(r_out)
384			bad_ind = np.array(bad_ind) >= 0
385
386			# Replace any bad indices with NaN, casting output as numpy arrays
387			if np.any(bad_ind):
388			lat_out[bad_ind] = np.nan
389			lon_out[bad_ind] = np.nan
390			r_out[bad_ind] = np.nan
391			except SystemError as serr:
392			aacgmv2.logger.warning('C Error encountered: {:}'.format(serr))
393
394			return lat_out, lon_out, r_out
395
396
397			def get_aacgm_coord(glat, glon, height, dtime, method="ALLOWTRACE"):
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)
430
431			# Get magnetic lat and lon.
432			mlat, mlon, _ = convert_latlon(glat, glon, height, dtime,
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]
437
438			return mlat, mlon, mlt
439
440
441			def get_aacgm_coord_arr(glat, glon, height, dtime, method="ALLOWTRACE"):
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)
475
476			# Get magnetic lat and lon.
477			mlat, mlon, _ = convert_latlon_arr(glat, glon, height, dtime,
478			method_code=method_code)
479
480			if np.any(np.isfinite(mlon)):
481			# Get magnetic local time
482			mlt = convert_mlt(mlon, dtime, m2a=False)
483			else:
484			mlt = np.full(shape=len(mlat), fill_value=np.nan)
485
486			return mlat, mlon, mlt
487
488
489			def convert_str_to_bit(method_code):
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,
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("\|")
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
525			if k in convert_code.keys()])
526
527			return bit_code
528
529
530			def convert_bool_to_bit(a2g=False, trace=False, allowtrace=False,
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
556
557			if trace:
558			bit_code += c_aacgmv2.TRACE
559			if allowtrace:
560			bit_code += c_aacgmv2.ALLOWTRACE
561			if badidea:
562			bit_code += c_aacgmv2.BADIDEA
563			if geocentric:
564			bit_code += c_aacgmv2.GEOCENTRIC
565
566			return bit_code
567
568
569			def convert_mlt(arr, dtime, m2a=False):
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)
595			if arr.shape == ():
596			arr = np.array([arr])
597
598			if len(arr.shape) > 1:
599			raise ValueError("unable to process multi-dimensional arrays")
600
601			# Test time
602			try:
603			dtime = test_time(dtime)
604			years = [dtime.year for dd in arr]
605			months = [dtime.month for dd in arr]
606			days = [dtime.day for dd in arr]
607			hours = [dtime.hour for dd in arr]
608			minutes = [dtime.minute for dd in arr]
609			seconds = [dtime.second for dd in arr]
610			except ValueError as verr:
611			dtime = np.asarray(dtime)
612			if dtime.shape == ():
613			raise ValueError(verr)
614			elif dtime.shape != arr.shape:
615			raise ValueError("array input for datetime and MLon/MLT must match")
616
617			years = [dd.year for dd in dtime]
618			months = [dd.month for dd in dtime]
619			days = [dd.day for dd in dtime]
620			hours = [dd.hour for dd in dtime]
621			minutes = [dd.minute for dd in dtime]
622			seconds = [dd.second for dd in dtime]
623
624			arr = list(arr)
625
626			# Calculate desired location, C routines set date and time
627			if m2a:
628			# Get the magnetic longitude
629			if len(arr) == 1:
630			out = c_aacgmv2.inv_mlt_convert(years[0], months[0], days[0],
631			hours[0], minutes[0], seconds[0],
632			arr[0])
633			else:
634			out = c_aacgmv2.inv_mlt_convert_arr(years, months, days, hours,
635			minutes, seconds, arr)
636			else:
637			# Get magnetic local time
638			if len(arr) == 1:
639			out = c_aacgmv2.mlt_convert(years[0], months[0], days[0], hours[0],
640			minutes[0], seconds[0], arr[0])
641			out = np.array([out])
642			else:
643			out = np.array(c_aacgmv2.mlt_convert_arr(years, months, days, hours,
644			minutes, seconds, arr))
645
646			return out
647

aburrell / aacgmv2

Pull Request — develop (#57)

aacgmv2.wrapper.convert_latlon() B

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like