aacgmv2.wrapper.convert_str_to_bit() - Code Metrics - aburrell/aacgmv2 - Measure and Improve Code Quality continuously with Scrutinizer

aacgmv2.wrapper.convert_str_to_bit() A
last analyzed 2024-12-27 21:19 UTC

↳ Parent: aacgmv2.wrapper

Complexity

Conditions

Size

Total Lines	45
Code Lines	9

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
eloc	9
dl	0
loc	45
rs	9.95
c	0
b	0
f	0
cc	1
nop	1

# 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."""

import datetime as dt
import numpy as np
import os

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 : any
        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 : bool
        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:
    1. Coefficients: 0-2000 km
    2. 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
    trace_opt = (TRACE | ALLOWTRACE | BADIDEA)
    if height > aacgmv2.high_alt_coeff and not (bit_code & trace_opt):
        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):
    """Set 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"):
    """Convert 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 : dt.datetime
        Datetime for magnetic field
    method_code : str or int
        Bit code or string denoting which type(s) of conversion to perform
        (default="G2A")

        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
    -------
    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 as err:
        estr = "".join(["unable to perform conversion at {:.1f}".format(in_lat),
                        ", {:.1f} {:.1f} km, {:}".format(in_lon, height, dtime),
                        " using method {:} <{:}>. Recall".format(bit_code, err),
                        " that AACGMV2 is undefined near the equator."])
        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"):
    """Convert between geomagnetic coordinates and AACGM coordinates.

    Parameters
    ----------
    in_lat : np.ndarray, 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 : dt.datetime
        Single datetime object for magnetic field
    method_code : int or str
        Bit code or string denoting which type(s) of conversion to perform
        (default="G2A")

        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
    -------
    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 : dt.datetime
        Date and time to calculate magnetic location
    method : str
        The type(s) of conversion to perform (default="ALLOWTRACE")

        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
    -------
    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 : dt.datetime
        Date and time to calculate magnetic location
    method : str
        The type(s) of conversion to perform (default="ALLOWTRACE")

        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
    -------
    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:

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

aburrell / aacgmv2

aacgmv2.wrapper.convert_str_to_bit() A last analyzed 2024-12-27 21:19 UTC

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like

aacgmv2.wrapper.convert_str_to_bit() A
last analyzed 2024-12-27 21:19 UTC