aacgmv2.wrapper.test_height() - Code Metrics - Inspection of "Very high altitudes with TRACE" - aburrell/aacgmv2 - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Pull Request — develop (#36)

by Angeline

created 2019-08-16 17:35 UTC

aacgmv2.wrapper.test_height() B

↳ Parent: aacgmv2.wrapper

Complexity

Conditions

Size

Total Lines	57
Code Lines	20

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
eloc	20
dl	0
loc	57
rs	8.4666
c	0
b	0
f	0
cc	6
nop	2

How to fix Long Method

# -*- coding: utf-8 -*-
"""Pythonic wrappers for AACGM-V2 C functions.

Functions
--------------
test_time : Test the time and ensure it is a datetime.datetime object
test_height : Test the height and see if it is appropriate for the method
set_coeff_path : Set the coefficient paths using default or supplied values
convert_latlon : Converts scalar location
convert_latlon_arr : Converts array location
get_aacgm_coord : Get scalar magnetic lat, lon, mlt from geographic location
get_aacgm_coord_arr : Get array magnetic lat, lon, mlt from geographic location
convert_str_to_bit : Convert human readible AACGM flag to bits
convert_bool_to_bit : Convert boolian flags to bits
convert_mlt : Get array mlt
--------------

"""

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

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

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
    TypeError or 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 as terr:
        raise TypeError("unable to set time for {:}: {:}".format(dtime, terr))
    except RuntimeError as rerr:
        raise RuntimeError("unable to set time for {:}: {:}".format(dtime,
                                                                    rerr))

    # convert location
    try:
        lat_out, lon_out, r_out = c_aacgmv2.convert(in_lat, in_lon, height,
                                                    bit_code)
    except:
        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
    TypeError or RuntimeError if unable to set AACGMV2 datetime

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

    """

    # 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 or int, create an array
    test_array = np.array([len(in_lat.shape), len(in_lon.shape),
                           len(height.shape)])
    if test_array.min() == 0:
        if test_array.max() == 0:
            aacgmv2.logger.warning("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:
            imax = test_array.argmax()
            max_shape = in_lat.shape if imax == 0 else (in_lon.shape \
                                            if imax == 1 else height.shape) 
            if not test_array[0]:
                in_lat = np.full(shape=max_shape, fill_value=in_lat)
            if not test_array[1]:
                in_lon = np.full(shape=max_shape, fill_value=in_lon)
            if not test_array[2]:
                height = np.full(shape=max_shape, 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):
        shape_dict = {'lat': in_lat.shape, 'lon': in_lon.shape,
                      'height': height.shape}
        ulen = np.unique(shape_dict.values())
        array_key = [kk for i, kk in enumerate(shape_dict.keys())
                     if shape_dict[kk] != (1,)]
        if len(array_key) == 3:
            raise ValueError('lat, lon, and height arrays are mismatched')
        elif len(array_key) == 2:
            if shape_dict[array_key[0]] == shape_dict[array_dict[1]]:

                raise ValueError('{:s} and {:s} arrays are mismatched'.format(\
                                                                *array_key))

    # 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 as terr:
        raise TypeError("unable to set time for {:}: {:}".format(dtime, terr))
    except RuntimeError as rerr:
        raise RuntimeError("unable to set time for {:}: {:}".format(dtime,
                                                                    rerr))

    # Vectorise the AACGM C routine
    convert_vectorised = np.vectorize(c_aacgmv2.convert)

    # convert
    try:
        lat_out, lon_out, r_out = convert_vectorised(in_lat, in_lon, height,
                                                     bit_code)

    except:
        err = sys.exc_info()[0]
        estr = "unable to perform vector conversion at {:} using ".format(dtime)
        estr = "{:s}method {:}: {:}".format(estr, bit_code, err)
        aacgmv2.logger.warning(estr)
        pass

    return lat_out, lon_out, r_out

def get_aacgm_coord(glat, glon, height, dtime, method="TRACE"):
    """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
    mlt = np.nan if np.isnan(mlon) else convert_mlt(mlon, dtime, m2a=False)

    return mlat, mlon, mlt


def get_aacgm_coord_arr(glat, glon, height, dtime, method="TRACE"):
    """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.all(np.isnan(mlon)):
        mlt = np.full(shape=mlat.shape, fill_value=np.nan)
    else:
        # Get magnetic local time
        mlt = convert_mlt(mlon, dtime, m2a=False)

    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_line or float)
        Magnetic longitudes (degrees E) or MLTs (hours) to convert
    dtime : (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.

    """

    # Test time
    dtime = test_time(dtime)
    
    # Calculate desired location, C routines set date and time
    if m2a:
        # Get the magnetic longitude
        inv_vectorised = np.vectorize(c_aacgmv2.inv_mlt_convert)
        out = inv_vectorised(dtime.year, dtime.month, dtime.day, dtime.hour,
                             dtime.minute, dtime.second, arr)
    else:
        # Get magnetic local time
        mlt_vectorised = np.vectorize(c_aacgmv2.mlt_convert)
        out = mlt_vectorised(dtime.year, dtime.month, dtime.day, dtime.hour,
                             dtime.minute, dtime.second, arr)

    if hasattr(out, "shape") and out.shape == ():
        out = float(out)

    return out


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

aburrell / aacgmv2

Pull Request — develop (#36)

aacgmv2.wrapper.test_height() B

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like