aacgmv2.wrapper.get_aacgm_coord_arr() - 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-07-22 17:54 UTC

aacgmv2.wrapper.get_aacgm_coord_arr() A

↳ Parent: aacgmv2.wrapper

Complexity

Conditions

Size

Total Lines	44
Code Lines	7

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
eloc	7
dl	0
loc	44
rs	10
c	0
b	0
f	0
cc	2
nop	5

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

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

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

    """
    import aacgmv2
    from aacgmv2._aacgmv2 import TRACE, ALLOWTRACE, BADIDEA

    # Test for heights that are allowed but not smart
    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)

    Returns
    ---------
    Void
    """
    import aacgmv2
    import os

    # 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, 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
    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)
    """
    import aacgmv2._aacgmv2 as c_aacgmv2
    import aacgmv2

    # Test time
    if isinstance(dtime, dt.date):
        dtime = dt.datetime.combine(dtime, dt.time(0))

    if not isinstance(dtime, dt.datetime):
        raise ValueError('time must be specified as datetime object')

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

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

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

    # Test height that may or may not cause failure
    good_height = test_height(height, bit_code)

    if not good_height:
        return lat_out, lon_out, r_out

    # Test latitude range
    if abs(in_lat) > 90.0:
        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:
        raise ValueError("unable to set time for {:}".format(dtime))

    # convert location
    try:
        lat_out, lon_out, r_out = c_aacgmv2.convert(in_lat, in_lon, height,
                                                    bit_code)
    except:
        pass

    return lat_out, lon_out, r_out

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

    Parameters
    ------------
    in_lat : (np.ndarray or list or float)
        Input latitude in degrees N (code specifies type of latitude)
    in_lon : (np.ndarray or list or float)
        Input longitude in degrees E (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
    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)

    Notes
    -------
    At least one of in_lat, in_lon, and height must be a list or array.
    """
    import aacgmv2._aacgmv2 as c_aacgmv2
    import aacgmv2

    # 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 test_array[0] == 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):
        ulen = np.unique([in_lat.shape, in_lon.shape, height.shape])
        if ulen.min() != (1,):
            raise ValueError('mismatched input arrays')

    # Test time
    if isinstance(dtime, dt.date):
        dtime = dt.datetime.combine(dtime, dt.time(0))

    if not isinstance(dtime, dt.datetime):
        raise ValueError('time must be specified as datetime object')

    # 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 code
    try:
        bit_code = convert_str_to_bit(code.upper())
    except AttributeError:
        bit_code = code

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

    # Test height
    good_height = test_height(np.nanmax(height), bit_code)

    if not good_height:
        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:
        raise RuntimeError("unable to set time for {:}".format(dtime))

    # Vectorise the AACGM code
    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:
        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 code
    code = "G2A|{:s}".format(method)

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

    # Get magnetic local time
    if np.isnan(mlon):
        mlt = np.nan
    else:
        mlt = 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 code
    code = "G2A|{:s}".format(method)

    # Get magnetic lat and lon.
    mlat, mlon, _ = convert_latlon_arr(glat, glon, height, dtime, code=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(code):
    """convert string code specification to bit code specification

    Parameters
    ------------
    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)
        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'.
    """
    import aacgmv2._aacgmv2 as c_aacgmv2

    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
    codes = code.upper().replace(" ", "").split("|")

    # Add the valid parts of the code, invalid elements are ignored
    bit_code = sum([convert_code[k] for k in 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
    """
    import aacgmv2._aacgmv2 as c_aacgmv2

    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.
    """
    import aacgmv2._aacgmv2 as c_aacgmv2

    # Test time
    if isinstance(dtime, dt.date):
        if not isinstance(dtime, dt.datetime):
            dtime = dt.datetime.combine(dtime, dt.time(0))
    else:
        raise ValueError('time must be specified as datetime object')

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

aburrell / aacgmv2

Pull Request — develop (#36)

aacgmv2.wrapper.get_aacgm_coord_arr() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like