Skip to content

Commit

Permalink
Merge pull request #10 from aburrell/scrutinizer-suggestions
Browse files Browse the repository at this point in the history
Scrutinizer suggestions.  Failures come from code quality checking, and are not all appropriate.
  • Loading branch information
Angeline Burrell committed Mar 20, 2018
2 parents 281488e + 7ac335b commit 391d80f
Show file tree
Hide file tree
Showing 16 changed files with 1,160 additions and 798 deletions.
2 changes: 1 addition & 1 deletion .bumpversion.cfg
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
[bumpversion]
current_version = 2.0.1
current_version = 2.4.0
commit = True
tag = True

Expand Down
4 changes: 3 additions & 1 deletion .travis.yml
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,9 @@ install:
- pip install coveralls
- "python setup.py install"
- pip install tox-travis
script: tox
script:
- tox
- coverage run --source aacgmv2 -m py.test
after_sucess: coveralls
notifications:
email:
Expand Down
23 changes: 12 additions & 11 deletions README.rst
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ Overview
|docs| |version|

This is a Python wrapper for the `AACGM-v2 C library
<https://engineering.dartmouth.edu/superdarn/aacgm.html>`_, which allows
<http://superdarn.thayer.dartmouth.edu/aacgm.html>`_, which allows
converting between geographic and magnetic coordinates. The currently included
version of the C library is 2.4. The package is free software
(MIT license).
Expand All @@ -21,25 +21,26 @@ Convert between AACGM and geographic coordinates::

>>> import aacgmv2
>>> import datetime as dt
>>> import numpy as np
>>> np.set_printoptions(formatter={'float_kind': lambda x:'{:.4f}'.format(x)})
>>> # geo to AACGM, single numbers
>>> dtime = dt.datetime(2013, 11, 3)
>>> mlat, mlon, mlt = aacgmv2.get_aacgm_coord(60, 15, 300, dtime)
>>> "{:.4f} {:.4f} {:.4f}".format(mlat, mlon, mlt)
'57.4698 93.6300 1.4822'
>>> np.array(aacgmv2.get_aacgm_coord(60, 15, 300, dtime))
array([57.4698, 93.6300, 1.4822])
>>> # AACGM to geo, mix arrays/numbers
>>> glat, glon, alt = aacgmv2.convert_latlon_arr([90, -90], 0, 0, dtime, code="A2G")
>>> ["{:.4f} {:.4f} {:.4f}".format(lat, glon[i], alt[i]) for i,lat in enumerate(glat)]
['82.9666 -84.6652 14.1244', '-74.3385 125.8401 12.8771']
>>> aacgmv2.convert_latlon_arr([90, -90], 0, 0, dtime, code="A2G")
(array([82.9666, -74.3385]), array([-84.6652, 125.8401]), array([14.1244, 12.8771]))

Convert between AACGM and MLT::

>>> import aacgmv2
>>> import datetime as dt
>>> import numpy as np
>>> np.set_printoptions(formatter={'float_kind': lambda x:'{:.4f}'.format(x)})
>>> # MLT to AACGM
>>> dtime = dt.datetime(2013, 11, 3, 0, 0, 0)
>>> mlon_check = aacgmv2.convert_mlt([1.4822189, 12], dtime, m2a=True)
>>> ["{:.4f}".format(lon) for lon in mlon_check]
['93.6300', '-108.6033']
>>> aacgmv2.convert_mlt([1.4822189, 12], dtime, m2a=True)
array([93.6300, -108.6033])

If you don't know or use Python, you can also use the command line. See details
in the full documentation.
Expand Down Expand Up @@ -96,7 +97,7 @@ Badges
:alt: Code Quality Status

.. |codacy| image:: https://img.shields.io/codacy/af7fdf6be28841f283dfdbc1c01fa82a.svg?style=flat
:target: https://www.codacy.com/app/aburrell/aacgmv2
:target: https://app.codacy.com/aburrell/aacgmv2
:alt: Codacy Code Quality Status

.. |codeclimate| image:: https://codeclimate.com/github/aburrell/aacgmv2/badges/gpa.svg
Expand Down
17 changes: 9 additions & 8 deletions aacgmv2/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
Parameters
---------------------------------------------------------------------------
AACGM_v2_DAT_PREFIX
AACGM_V2_DAT_PREFIX
IGRF_12_COEFFS
_aacgmv2.G2A
_aacgmv2.A2G
Expand All @@ -27,8 +27,8 @@
get_aacgm_coord_arr
convert
convert_mlt
subsol
set_coeff_path
wrapper.set_coeff_path
deprecated.subsol
_aacgmv2.convert
_aacgmv2.set_datetime
_aacgmv2.mlt_convert
Expand All @@ -39,16 +39,17 @@
Modules
---------------------------------------------------------------------------
depricated
deprecated
wrapper
_aacgmv2
---------------------------------------------------------------------------
"""
import os.path as _path
import logbook as logging
__version__ = "2.0.1"
__version__ = "2.4.0"

# path and filename prefix for the IGRF coefficients
AACGM_v2_DAT_PREFIX = _path.join(_path.realpath(_path.dirname(__file__)),
AACGM_V2_DAT_PREFIX = _path.join(_path.realpath(_path.dirname(__file__)),
'aacgm_coeffs', 'aacgm_coeffs-12-')
IGRF_12_COEFFS = _path.join(_path.realpath(_path.dirname(__file__)),
'igrf12coeffs.txt')
Expand All @@ -64,8 +65,8 @@
logging.exception(__file__ + ' -> aacgmv2: ' + str(err))

try:
from aacgmv2 import (depricated)
from aacgmv2.depricated import (convert, subsol, set_coeff_path)
from aacgmv2 import (deprecated)
from aacgmv2.deprecated import (convert)
except Exception as err:
logging.exception(__file__ + ' -> aacgmv2: ' + str(err))

Expand Down
29 changes: 15 additions & 14 deletions aacgmv2/__main__.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,27 +26,27 @@ def main():

desc = 'Converts between geographical coordinates, AACGM-v2, and MLT'
parser = argparse.ArgumentParser(description=desc)

desc = 'for help, run %(prog)s SUBCOMMAND -h'
subparsers = parser.add_subparsers(title='Subcommands', prog='aacgmv2',
dest='subcommand', description=desc)
subparsers.required = True

desc = 'convert to/from geomagnetic coordinates. Input file must have lines'
desc += 'of the form "LAT LON ALT".'
parser_convert = subparsers.add_parser('convert', help=(desc))

desc = 'convert between magnetic local time (MLT) and AACGM-v2 longitude. '
desc += 'Input file must have a single number on each line.'
parser_convert_mlt = subparsers.add_parser('convert_mlt', help=(desc))

desc = 'input file (stdin if none specified)'
for p in [parser_convert, parser_convert_mlt]:
p.add_argument('-i', '--input', dest='file_in', metavar='FILE_IN',
type=argparse.FileType('r'), default=STDIN, help=desc)
p.add_argument('-o', '--output', dest='file_out', metavar='FILE_OUT',
type=argparse.FileType('wb'), default=STDOUT,
help='output file (stdout if none specified)')
for pp in [parser_convert, parser_convert_mlt]:
pp.add_argument('-i', '--input', dest='file_in', metavar='FILE_IN',
type=argparse.FileType('r'), default=STDIN, help=desc)
pp.add_argument('-o', '--output', dest='file_out', metavar='FILE_OUT',
type=argparse.FileType('wb'), default=STDOUT,
help='output file (stdout if none specified)')

desc = 'date for magnetic field model (1900-2020, default: today)'
parser_convert.add_argument('-d', '--date', dest='date', metavar='YYYYMMDD',
Expand Down Expand Up @@ -91,13 +91,14 @@ def main():
allowtrace=args.allowtrace,
badidea=args.badidea,
geocentric=args.geocentric)
lats, lons, rs = aacgmv2.convert_latlon_arr(array[:,0], array[:,1],
array[:,2], dtime=date,
code=code)
np.savetxt(args.file_out, np.column_stack((lats, lons, rs)), fmt='%.8f')
lats, lons, alts = aacgmv2.convert_latlon_arr(array[:, 0], array[:, 1],
array[:, 2], dtime=date,
code=code)
np.savetxt(args.file_out, np.column_stack((lats, lons, alts)),
fmt='%.8f')
elif args.subcommand == 'convert_mlt':
dtime = dt.datetime.strptime(args.datetime, '%Y%m%d%H%M%S')
out = aacgmv2.convert_mlt(array[:,0], dtime, m2a=args.m2a)
out = aacgmv2.convert_mlt(array[:, 0], dtime, m2a=args.m2a)
np.savetxt(args.file_out, out, fmt='%.8f')


Expand Down
90 changes: 36 additions & 54 deletions aacgmv2/depricated.py → aacgmv2/deprecated.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@
Functions
-------------------------------------------------------------------------------
convert : Converts array location
set_coeff_path : Previously set environment variables, no longer used
subsol : finds subsolar geocentric longitude and latitude
gc2gd_lat : Convert between geocentric and geodetic coordinates
igrf_dipole_axis : Get Cartesian unit vector pointing at the IGRF north dipole
Expand All @@ -21,7 +20,6 @@
from __future__ import division, absolute_import, unicode_literals
import numpy as np
import logbook as logging
import aacgmv2

def convert(lat, lon, alt, date=None, a2g=False, trace=False, allowtrace=False,
badidea=False, geocentric=False):
Expand Down Expand Up @@ -56,35 +54,29 @@ def convert(lat, lon, alt, date=None, a2g=False, trace=False, allowtrace=False,
lon_out : (float)
Output longitude in degrees E
"""
import aacgmv2

if(np.array(alt).max() > 2000 and not trace and not allowtrace and
badidea):
estr = 'coefficients are not valid for altitudes above 2000 km. You'
estr += ' must either use field-line tracing (trace=True '
estr += 'or allowtrace=True) or indicate you know this is a bad idea'
logging.error(estr)
raise ValueError

# construct a code from the boolian flags
bit_code = aacgmv2.convert_bool_to_bit(a2g=a2g, trace=trace,
allowtrace=allowtrace,
badidea=badidea,
geocentric=geocentric)

# convert location
lat_out, lon_out, r_out = aacgmv2.convert_latlon_arr(lat, lon, alt, date,
code=bit_code)
lat_out, lon_out, _ = aacgmv2.convert_latlon_arr(lat, lon, alt, date,
code=bit_code)

return lat_out, lon_out

def set_coeff_path():
"""This depricated routine used to set environment variables, and now is
not needed.
"""

logging.warning("this routine is no longer needed")
return

def subsol(year, doy, ut):
def subsol(year, doy, utime):
"""Finds subsolar geocentric longitude and latitude.
Parameters
Expand All @@ -93,7 +85,7 @@ def subsol(year, doy, ut):
Calendar year between 1601 and 2100
doy : (int)
Day of year between 1-365/366
ut : (float)
utime : (float)
Seconds since midnight on the specified day
Returns
Expand All @@ -116,7 +108,7 @@ def subsol(year, doy, ut):
After Fortran code by A. D. Richmond, NCAR. Translated from IDL
by K. Laundal.
"""
yr = year - 2000
yr2 = year - 2000

if year >= 2101:
logging.error('subsol invalid after 2100. Input year is:', year)
Expand All @@ -130,55 +122,44 @@ def subsol(year, doy, ut):
ncent = 3 - ncent
nleap = nleap + ncent

l0 = -79.549 + (-0.238699 * (yr - 4 * nleap) + 3.08514e-2 * nleap)
g0 = -2.472 + (-0.2558905 * (yr - 4 * nleap) - 3.79617e-2 * nleap)
l_0 = -79.549 + (-0.238699 * (yr2 - 4 * nleap) + 3.08514e-2 * nleap)
g_0 = -2.472 + (-0.2558905 * (yr2 - 4 * nleap) - 3.79617e-2 * nleap)

# Days (including fraction) since 12 UT on January 1 of IYR:
df = (ut / 86400 - 1.5) + doy

# Addition to Mean longitude of Sun since January 1 of IYR:
lf = 0.9856474 * df

# Addition to Mean anomaly since January 1 of IYR:
gf = 0.9856003 * df
# Days (including fraction) since 12 UT on January 1 of IYR2:
dfrac = (utime / 86400 - 1.5) + doy

# Mean longitude of Sun:
l = l0 + lf
l_sun = l_0 + 0.9856474 * dfrac

# Mean anomaly:
grad = np.radians(g0 + gf)
grad = np.radians(g_0 + 0.9856003 * dfrac)

# Ecliptic longitude:
lmrad = np.radians(l + 1.915 * np.sin(grad) + 0.020 * np.sin(2 * grad))
lmrad = np.radians(l_sun + 1.915 * np.sin(grad) + 0.020 * np.sin(2 * grad))
sinlm = np.sin(lmrad)

# Days (including fraction) since 12 UT on January 1 of 2000:
n = df + 365.0 * yr + nleap
epoch_day = dfrac + 365.0 * yr2 + nleap

# Obliquity of ecliptic:
epsrad = np.radians(23.439 - 4.0e-7 * n)
epsrad = np.radians(23.439 - 4.0e-7 * epoch_day)

# Right ascension:
alpha = np.degrees(np.arctan2(np.cos(epsrad) * sinlm, np.cos(lmrad)))

# Declination:
delta = np.degrees(np.arcsin(np.sin(epsrad) * sinlm))

# Subsolar latitude:
sbsllat = delta
# Declination, which is the subsolar latitude:
sbsllat = np.degrees(np.arcsin(np.sin(epsrad) * sinlm))

# Equation of time (degrees):
etdeg = l - alpha
nrot = np.round(etdeg / 360.0)
etdeg = etdeg - 360.0 * nrot
etdeg = l_sun - alpha
etdeg = etdeg - 360.0 * np.round(etdeg / 360.0)

# Apparent time (degrees):
aptime = ut / 240.0 + etdeg # Earth rotates one degree every 240 s.
aptime = utime / 240.0 + etdeg # Earth rotates one degree every 240 s.

# Subsolar longitude:
sbsllon = 180.0 - aptime
nrot = np.round(sbsllon / 360.0)
sbsllon = sbsllon - 360.0 * nrot
sbsllon = sbsllon - 360.0 * np.round(sbsllon / 360.0)

return sbsllon, sbsllat

Expand All @@ -195,8 +176,8 @@ def gc2gd_lat(gc_lat):
gd_lat : (same as input)
Geodetic latitude in degrees N
"""
WGS84_e2 = 0.006694379990141317 - 1.0
return np.rad2deg(-np.arctan(np.tan(np.deg2rad(gc_lat)) / WGS84_e2))
wgs84_e2 = 0.006694379990141317 - 1.0
return np.rad2deg(-np.arctan(np.tan(np.deg2rad(gc_lat)) / wgs84_e2))

def igrf_dipole_axis(date):
"""Get Cartesian unit vector pointing at dipole pole in the north,
Expand All @@ -220,6 +201,7 @@ def igrf_dipole_axis(date):
date, or extrapolated if date > latest IGRF model
"""
import datetime as dt
import aacgmv2

# get time in years, as float:
year = date.year
Expand All @@ -228,8 +210,8 @@ def igrf_dipole_axis(date):
year = year + doy / year_days

# read the IGRF coefficients
with open(aacgmv2.IGRF_12_COEFFS, 'r') as f:
lines = f.readlines()
with open(aacgmv2.IGRF_12_COEFFS, 'r') as f_igrf:
lines = f_igrf.readlines()

years = lines[3].split()[3:][:-1]
years = np.array(years, dtype=float) # time array
Expand All @@ -256,15 +238,15 @@ def igrf_dipole_axis(date):
h11 = np.interp(year, years, h11)
else:
# extrapolation
dt = year - years[-1]
g10 = g10[-1] + g10sv * dt
g11 = g11[-1] + g11sv * dt
h11 = h11[-1] + h11sv * dt
dyear = year - years[-1]
g10 = g10[-1] + g10sv * dyear
g11 = g11[-1] + g11sv * dyear
h11 = h11[-1] + h11sv * dyear

# calculate pole position
B0 = np.sqrt(g10**2 + g11**2 + h11**2)
B_0 = np.sqrt(g10**2 + g11**2 + h11**2)

# Calculate output
m = -np.array([g11, h11, g10]) / B0
return m
m_0 = -np.array([g11, h11, g10]) / B_0

return m_0

0 comments on commit 391d80f

Please sign in to comment.