Merge branch 'develop' into smag_bugfix

Changelog.rst

@@ -22,6 +22,7 @@ Summary of all changes made since the first stable release
 * ENH: Added function to select data along a satellite track
 * ENH: Changed attributes in VectorData into properties to ensure expected
        behaviour if altering the class data after initialisation
+* ENH: Added IMAGE SI12, SI13, and WIC DMSP corrections to `harmonic`
 * MAINT: Removed support for Python 2.7, 3.5, and 3.6; added support for 3.10
 * MAINT: Improved PEP8 compliance
 * MAINT: Updated pysat routines to v3.0.0 standards

docs/citing.rst

@@ -42,8 +42,8 @@ which may also be cited if a description of the package is desired.
 IMAGE FUV Boundaries
 --------------------
 
-Please cite both the papers discussing both the instrument and the boundary
-retrieval method.
+Please cite both the papers discussing the instrument and the appropriate
+boundary retrieval method.
 
 * **SI12/SI13**: Mende, S., et al. Space Science Reviews (2000) 91: 287-318.
   http://doi.org/10.1023/A:1005292301251.
@@ -53,6 +53,8 @@ retrieval method.
   high‐latitude ionospheric climatologies and empirical models,
   J. Geophys. Res. Space Physics, 122, 932–947,
   http://doi.org/10.1002/2016JA023235.
+* **OCB**: Chisham, G. et al. (2022) Ionospheric Boundaries Derived from Auroral
+  Images. In Prep.
 * **OCB**: Chisham, G. (2017) Auroral Boundary Derived from IMAGE Satellite
   Mission Data (May 2000 - Oct 2002), Version 1.1, Polar Data Centre, Natural
   Environment Research Council, UK.

docs/ocb_datasets.rst

@@ -19,16 +19,24 @@ provided in the :py:mod:`ocbpy.boundaries.files` sub-module.
 IMAGE
 -----
 
-Data from three auroral instruments provide northern hemisphere OCB and EAB
-locations for 3 May 2000 02:41:43 UT - 31 Oct 2002 20:05:16, though not all of
-the times included in these files contain high-quality estimations of the
-boundary locations. Recommended selection criteria are included as defaults in
-the :py:class:`~ocbpy.OCBoundary` class. There are also boundary
-files that combine the information from all instruments to obtain the OCB and
-EAB. You can read more about the OCB determination, EAB determination, this
-selection criteria, and the three auroral instruments (IMAGE Wideband Imaging
-Camera (WIC) and FUV Spectrographic Imagers SI12 and SI13) in the articles
-listed in :ref:`cite-image`.
+Data from three auroral instruments provide northern hemisphere poleward auroral
+boundary (PAB) and EAB locations for 3 May 2000 02:41:43 UT - 31 Oct 2002
+20:05:16, though not all of the times included in these files contain
+high-quality estimations of the boundary locations. Recommended selection
+criteria are included as defaults in the :py:class:`~ocbpy.OCBoundary` class.
+There are also boundary files that combine the information from all instruments
+to obtain the OCB and EAB. These combined files are the default boundaries for
+the IMAGE time period.  You can read more about the OCB determination, EAB
+determination, this selection criteria, and the three auroral instruments
+(IMAGE Wideband Imaging Camera (WIC) and FUV Spectrographic Imagers SI12 and
+SI13) in the articles listed in :ref:`cite-image`.
+
+The most recent corrects for each instrument that add the DMSP particle
+precipitation corrections to the PAB and EAB locations are included in
+:py:mod:`ocbpy.ocb_correction`.  These corrections should be applied to the
+data used to obtain the circle fits included in the instrument files, not the
+circle fits themselves. These data sets may be obtained from the British
+Antarctic Survey.
 
 
 .. _bound-data-ampere:

ocbpy/boundaries/README.md

@@ -1,8 +1,8 @@
 This directory contains files with Open Closed field line Boundaries obtained
 from different instruments
 
-IMAGE (si12/si13/wic) File Format
----------------------------------
+IMAGE (image/si12/si13/wic) File Format
+---------------------------------------
 YR, SOY, NB, PHICENT, RCENT, R, A, R_ERR, R_MERIT
 
 YR      : Year
@@ -20,7 +20,9 @@ R_MERIT : Radial distance from the most typical pole location in degrees
 There are certain ranges for NB, RCENT, and R that you shouldn’t use that can
 be found (and explained) in Chisham (2017), doi:10.1002/2016JA023235.  These
 ranges are the defaults in OCBoundary.get_next_good_ocb_ind.  When using these
-boundaries, remember to cite Chisham (2017).
+boundaries, remember to cite Chisham (2017). From ocbpy version 0.3.0 onward,
+the SI12, SI13, and WIC files contain uncorrected auroral boundary fits, while
+the IMAGE file contains DMSP corrected average boundaries for the OCB and EAB.
 
 AMPERE (amp) File Format
 ------------------------

ocbpy/ocb_correction.py

@@ -9,6 +9,8 @@
 ----------
 .. [4] Burrell, A. G. et al.: AMPERE Polar Cap Boundaries, Ann. Geophys., 38,
    481-490, doi:10.5194/angeo-38-481-2020, 2020.
+.. [6] Chisham, G. et al.: Ionospheric Boundaries Derived from Auroral Images,
+    in prep, 2022.
 
 """
 
@@ -18,7 +20,7 @@
 
 
 def circular(mlt, r_add=0.0):
-    """Return a circular boundary correction
+    """Return a circular boundary correction.
 
     Parameters
     ----------
@@ -42,7 +44,7 @@ def circular(mlt, r_add=0.0):
 
 
 def elliptical(mlt, instrument='ampere', method='median'):
-    """Return the results of an elliptical correction to the data boundary [4]_
+    """Return the results of an elliptical correction to the data boundary.
 
     Parameters
     ----------
@@ -59,6 +61,10 @@ def elliptical(mlt, instrument='ampere', method='median'):
     r_corr : float or array-like
         Radius correction in degrees at this MLT
 
+    References
+    ----------
+    Prefered AMPERE boundary correction validated in [4]_.
+
     """
 
     if instrument.lower() != 'ampere':
@@ -81,7 +87,7 @@ def elliptical(mlt, instrument='ampere', method='median'):
 
 
 def harmonic(mlt, instrument='ampere', method='median'):
-    """Return the results of a harmonic fit correction to the data boundary [4]_
+    """Return the results of a harmonic fit correction to the data boundary.
 
     Parameters
     ----------
@@ -91,38 +97,67 @@ def harmonic(mlt, instrument='ampere', method='median'):
         Data set's instrument name (default='ampere')
     method : str
         Method used to determine coefficients; accepts 'median' or
-        'gaussian' (default='median')
+        'gaussian' when `instrument` is 'ampere'.  Otherwise, accepts 'eab' or
+        'ocb'. (default='median')
 
     Returns
     -------
     r_corr : float or array-like
         Radius correction in degrees at this MLT
 
+    References
+    ----------
+    AMPERE boundaries obtained from [4]_. IMAGE boundaries obtained from [6]_.
+
     """
-    if instrument.lower() != 'ampere':
+    # Define the harmonic coefficients for each instrument and method/location
+    coeff = {'ampere': {'median': [3.31000535, -0.5452934, -1.24389141,
+                                   2.42619653, -0.66677988, -1.03467488,
+                                   -0.30763009, 0.52426756, 0.04359299,
+                                   0.60201848, 0.50618522, 1.04360529,
+                                   0.25186405],
+                        'gaussian': [3.80100827, 0.98555723, -3.43760943,
+                                     1.85084271, -0.36730751, -0.81975654,
+                                     -1.02823832, 1.30637288, -0.53599218,
+                                     0.40380183, -1.22462708, -1.2733629,
+                                     -0.62743381]},
+             'si12': {'ocb': [0.0405, -1.5078, 1.0, 0.5095, 1.0, 0.9371, 1.0,
+                              0.0708, 1.0, 0.0, 1.0, 0.0, 1.0],
+                      'eab': [-0.1447, -1.9779, 1.0, 2.6799, 1.0, 0.5778, 1.0,
+                              -1.2297, 1.0, 0.0, 1.0, 0.0, 1.0]},
+             'si13': {'ocb': [0.5797, -0.6837, 1.0, -0.5020, 1.0, 0.2971, 1.0,
+                              -0.4173, 1.0, 0.0, 1.0, 0.0, 1.0],
+                      'eab': [0.2500, -2.9931, 1.0, 0.8818, 1.0, 0.8511, 1.0,
+                              -0.6300, 1.0, 0.0, 1.0, 0.0, 1.0]},
+             'wic': {'ocb': [1.0298, -1.1249, 1.0, -0.7380, 1.0, 0.1838, 1.0,
+                             -0.6171, 1.0, 0.0, 1.0, 0.0, 1.0],
+                     'eab': [-0.4935, -2.1186, 1.0, 0.3188, 1.0, 0.5749, 1.0,
+                             -0.3118, 1.0, 0.0, 1.0, 0.0, 1.0]}}
+
+    # Check the inputs
+    inst = instrument.lower()
+    if inst not in coeff.keys():
         raise ValueError("no harmonic correction for {:}".format(instrument))
 
     method = method.lower()
-    coeff = {'median': [3.31000535, -0.5452934, -1.24389141, 2.42619653,
-                        -0.66677988, -1.03467488, -0.30763009, 0.52426756,
-                        0.04359299, 0.60201848, 0.50618522, 1.04360529,
-                        0.25186405],
-             'gaussian': [3.80100827, 0.98555723, -3.43760943, 1.85084271,
-                          -0.36730751, -0.81975654, -1.02823832, 1.30637288,
-                          -0.53599218, 0.40380183, -1.22462708, -1.2733629,
-                          -0.62743381]}
-
-    if method not in coeff.keys():
-        raise ValueError("unknown coefficient computation method")
+    if method not in coeff[inst].keys():
+        raise ValueError("".join(["unknown coefficient computation method, ",
+                                  "expects one of: {:}".format(
+                                      coeff[inst].keys())]))
 
+    # Calculate the offset
     rad_mlt = hr2rad(mlt)
-    r_corr = coeff[method][0] \
-        + coeff[method][1] * np.cos(rad_mlt + coeff[method][2]) \
-        + coeff[method][3] * np.sin(rad_mlt + coeff[method][4]) \
-        + coeff[method][5] * np.cos(2.0 * (rad_mlt + coeff[method][6])) \
-        + coeff[method][7] * np.sin(2.0 * (rad_mlt + coeff[method][8])) \
-        + coeff[method][9] * np.cos(3.0 * (rad_mlt + coeff[method][10])) \
-        + coeff[method][11] * np.sin(3.0 * (rad_mlt + coeff[method][12]))
+    r_corr = coeff[inst][method][0] \
+        + coeff[inst][method][1] * np.cos(rad_mlt + coeff[inst][method][2]) \
+        + coeff[inst][method][3] * np.sin(rad_mlt + coeff[inst][method][4]) \
+        + coeff[inst][method][5] * np.cos(2.0 * (
+            rad_mlt + coeff[inst][method][6])) \
+        + coeff[inst][method][7] * np.sin(2.0 * (
+            rad_mlt + coeff[inst][method][8])) \
+        + coeff[inst][method][9] * np.cos(3.0 * (
+            rad_mlt + coeff[inst][method][10])) \
+        + coeff[inst][method][11] * np.sin(3.0 * (
+            rad_mlt + coeff[inst][method][12]))
 
     # Because this is a poleward shift, return the negative of the correction
     return -r_corr

ocbpy/tests/test_ocb_correction.py

@@ -3,8 +3,7 @@
 # Copyright (C) 2017, AGB & GC
 # Full license can be found in License.md
 # -----------------------------------------------------------------------------
-""" Tests the ocboundary class and functions
-"""
+"""Tests the ocboundary class and functions."""
 
 import numpy as np
 import unittest
@@ -13,38 +12,48 @@
 
 
 class TestOCBCorrectionFailure(unittest.TestCase):
+    """Unit tests for correction failure evaluations."""
+
     def setUp(self):
-        """ Set up the test runs """
+        """Set up the test runs."""
         self.mlt = 12.0
         self.bad_kwarg = 'bad_kwarg'
         self.functions = {'elliptical': ocb_cor.elliptical,
                           'harmonic': ocb_cor.harmonic}
         self.bound = None
+        return
 
     def tearDown(self):
-        """ Tear down the test runs """
+        """Tear down the test runs."""
         del self.mlt, self.bad_kwarg, self.functions, self.bound
+        return
 
     def test_instrument_failure(self):
-        """ Test failure when an unknown instrument is provided """
+        """Test failure when an unknown instrument is provided."""
 
         for self.bound in self.functions.keys():
             with self.subTest(bound=self.bound):
-                with self.assertRaises(ValueError):
+                msg = "no {:s} correction for".format(self.bound)
+                with self.assertRaisesRegex(ValueError, msg):
                     self.functions[self.bound](self.mlt,
                                                instrument=self.bad_kwarg)
+        return
 
     def test_method_failure(self):
-        """ Test failure when an unknown method is provided """
+        """Test failure when an unknown method is provided."""
+        msg = "unknown coefficient computation method"
         for self.bound in self.functions.keys():
             with self.subTest(bound=self.bound):
-                with self.assertRaises(ValueError):
+                with self.assertRaisesRegex(ValueError, msg):
                     self.functions[self.bound](self.mlt, method=self.bad_kwarg)
+        return
 
 
 class TestOCBCorrection(unittest.TestCase):
+    """Unit tests for the boundary correction functions."""
+
     def setUp(self):
-        """ Set up test runs """
+        """Set up test runs."""
         self.functions = {'circular': ocb_cor.circular,
                           'elliptical': ocb_cor.elliptical,
                           'harmonic': ocb_cor.harmonic}
@@ -56,36 +65,65 @@ def setUp(self):
                                                   -3.4392638193624325])}
         self.gaus_results = {'elliptical': -2.51643691301747,
                              'harmonic': -2.293294645880221}
+        self.image_results = {'si12': {'ocb': np.array([0.67103129,
+                                                        -0.10084541]),
+                                       'eab': np.array([0.31691853,
+                                                        2.68970685])},
+                              'si13': {'ocb': np.array([0.71521016,
+                                                        -0.86843608]),
+                                       'eab': np.array([1.55220967,
+                                                        -0.19812977])},
+                              'wic': {'ocb': np.array([0.83660688,
+                                                       -1.62097642]),
+                                      'eab': np.array([1.89268527,
+                                                       0.13983824])}}
         self.bound = None
+        return
 
     def tearDown(self):
-        """ Clean up after each test """
+        """Clean up after each test."""
         del self.mlt, self.functions, self.def_results, self.gaus_results
-        del self.bound
+        del self.bound, self.image_results
+        return
 
     def test_default_float(self):
-        """ Test the boundary functions using a float and function defaults"""
+        """Test the boundary functions using a float and function defaults."""
         for self.bound in self.functions.keys():
             with self.subTest(bound=self.bound):
                 self.assertEqual(self.functions[self.bound](self.mlt[0]),
                                  self.def_results[self.bound][0])
+        return
 
     def test_default_arr(self):
-        """ Test the boundary functions using an array and function defaults"""
+        """Test the boundary functions using an array and function defaults."""
         for self.bound in self.functions.keys():
             with self.subTest(bound=self.bound):
                 self.assertTrue(np.all(
                     abs(self.functions[self.bound](self.mlt)
                         - self.def_results[self.bound]) < 1.0e-7))
+        return
 
     def test_circular_offset(self):
-        """ Test the circular boundary function with an offset """
+        """Test the circular boundary function with an offset."""
         self.assertEqual(ocb_cor.circular(self.mlt[0], r_add=1.0), 1.0)
+        return
 
     def test_gauss_method(self):
-        """ Test the boundary functions using an array and function defaults"""
+        """Test the boundary functions using an array and function defaults."""
         for self.bound in self.gaus_results.keys():
             with self.subTest(bound=self.bound):
                 self.assertAlmostEqual(
                     self.functions[self.bound](self.mlt[0], method='gaussian'),
                     self.gaus_results[self.bound])
+        return
+
+    def test_image_harmonic(self):
+        """Test the IMAGE harmonic correction functions."""
+        for self.bound in self.image_results.keys():
+            for method in self.image_results[self.bound].keys():
+                with self.subTest(bound=self.bound, method=method):
+                    self.assertTrue(np.all(
+                        abs(self.functions["harmonic"](
+                            self.mlt, instrument=self.bound, method=method)
+                            - self.image_results[self.bound][method]) < 1.0e-7))
+        return