diff --git a/Changelog.rst b/Changelog.rst index 42805bca..2a37d76b 100644 --- a/Changelog.rst +++ b/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 diff --git a/docs/citing.rst b/docs/citing.rst index 9a60f6a2..f0da9fd5 100644 --- a/docs/citing.rst +++ b/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. diff --git a/docs/ocb_datasets.rst b/docs/ocb_datasets.rst index 56fdffb8..08138a24 100644 --- a/docs/ocb_datasets.rst +++ b/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: diff --git a/ocbpy/boundaries/README.md b/ocbpy/boundaries/README.md index ad99f4b1..9d42e781 100644 --- a/ocbpy/boundaries/README.md +++ b/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 ------------------------ diff --git a/ocbpy/ocb_correction.py b/ocbpy/ocb_correction.py index ba7cf8a0..b73a0e6a 100644 --- a/ocbpy/ocb_correction.py +++ b/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 diff --git a/ocbpy/tests/test_ocb_correction.py b/ocbpy/tests/test_ocb_correction.py index fc2f6487..787cc62c 100644 --- a/ocbpy/tests/test_ocb_correction.py +++ b/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