diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 36f417fc..1d0bb489 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -50,33 +50,14 @@ jobs: coverage report coverage xml --rcfile=pyproject.toml - - name: Install and run Coveralls Reporter(Linux) - if: startsWith(matrix.os, 'ubuntu') - env: - COVERALLS_REPO_TOKEN: ${{ secrets.COVERALLS_REPO_TOKEN }} - COVERALLS_PARALLEL: true - run: | - curl -sL https://coveralls.io/coveralls-linux.tar.gz | tar -xz - ./coveralls report -f coverage.xml --parallel --repo-token=${{ secrets.COVERALLS_REPO_TOKEN }} --build-number ${{ github.run_number }} - - - name: Install and run Coveralls Reporter (Windows) - if: startsWith(matrix.os, 'windows') - env: - COVERALLS_REPO_TOKEN: ${{ secrets.COVERALLS_REPO_TOKEN }} - COVERALLS_PARALLEL: true - run: | - curl -L https://github.com/coverallsapp/coverage-reporter/releases/latest/download/coveralls-windows.exe -o coveralls.exe - ./coveralls.exe report -f coverage.xml --parallel --repo-token=${{ secrets.COVERALLS_REPO_TOKEN }} --build-number ${{ github.run_number }} - - - name: Report and run Coveralls (macOS) - if: startsWith(matrix.os, 'macos') - env: - COVERALLS_REPO_TOKEN: ${{ secrets.COVERALLS_REPO_TOKEN }} - COVERALLS_PARALLEL: true - run: | - brew tap coverallsapp/coveralls --quiet - brew install coveralls --quiet - coveralls report -f coverage.xml --parallel --repo-token=${{ secrets.COVERALLS_REPO_TOKEN }} --build-number ${{ github.run_number }} + - name: Coveralls Parallel + uses: coverallsapp/github-action@v2 + with: + github-token: ${{ secrets.GITHUB_TOKEN }} + flag-name: run=${{ join(matrix.*, '-') }} + parallel: true + format: cobertura + debug: true finish: name: Finish Coverage Analysis diff --git a/Changelog.rst b/Changelog.rst index 09b9bbb5..66184879 100644 --- a/Changelog.rst +++ b/Changelog.rst @@ -3,6 +3,10 @@ Changelog Summary of all changes made since the first stable release +0.7.0 (XX-XX-2025) +------------------ +* ENH: Added the Gussenhoven (1983) model for the EAB + 0.6.0 (07-07-2025) ------------------ * ENH: Updated the AMPERE boundaries to include 2022-2024, inclusive diff --git a/docs/citing.rst b/docs/citing.rst index fee63c5d..1e4fd674 100644 --- a/docs/citing.rst +++ b/docs/citing.rst @@ -112,3 +112,17 @@ when publishing your work. * Starkov, G. V. (1994) Mathematical model of the auroral boundaries, Geomagnetism and Aeronomy, English Translation, 34(3), 331-336. + + +.. _cite-guss: + +Gussenhoven Model +----------------- + +The Gussenhoven mathematical diffuse auroral boundary model is described in the +following article. If you use this model please cite both it and your Kp data +source when publishing your work. + +* Gussenhoven, M. S. et al. (1983) Systematics of the Equatorward Diffuse + Auroral Boundary, J. Geophys. Res, 88(A7), 5692-5708. + diff --git a/docs/examples/ex_model_boundaries.rst b/docs/examples/ex_model_boundaries.rst index a0731001..f4925b5c 100644 --- a/docs/examples/ex_model_boundaries.rst +++ b/docs/examples/ex_model_boundaries.rst @@ -32,7 +32,7 @@ be calculated for a set of MLTs using the This provides the AACGMV2 magnetic latitude for the entire range of MLT at the specified AL values. We can plot these boundaries, using the formatting function -previously defined in Example :ref:`format-polar-axes`. +previously defined in :ref:`set_up_polar_plot `. :: diff --git a/docs/ocb_models.rst b/docs/ocb_models.rst index fc4b3652..e8e5be27 100644 --- a/docs/ocb_models.rst +++ b/docs/ocb_models.rst @@ -27,6 +27,19 @@ in the code here using the boundary keyworkds: 'ocb', 'eab', and 'diffuse', respectively. +.. _bound-model-guss: + +Gussenhoven +----------- + +The Gussenhoven 1983 model (see :ref:`cite-guss`) uses a mathematical +formulation based on DMSP data and the Kp index. They specify a single boundary: +the equatorward edge of the diffuse aurora. As this model defines the boundary +using separate linear fits for different MLT bins, the results may be returned +at the binned times ('binned'), at the nearest binned time ('closest'), or at +the requested times using a circle fit to all of binned times ('circle'). + + .. _bound-model-module: Boundary Models Module diff --git a/ocbpy/boundaries/models.py b/ocbpy/boundaries/models.py index 4c36d7fa..fafc9457 100644 --- a/ocbpy/boundaries/models.py +++ b/ocbpy/boundaries/models.py @@ -9,11 +9,17 @@ ---------- .. [8] Starkov, G. V. (1994) Mathematical model of the auroral boundaries, Geomagnetism and Aeronomy, English Translation, 34(3), 331-336. +.. [9] Gussenhoven, M. S. et al. (1983) Systematics of the Equatorward Diffuse + Auroral Boundary, J. Geophys. Res, 88(A7), 5692-5708. +.. [10] Umbach and Jones (2003) A Few Methods for Fitting Circles to Data, + IEEE Transactions on Instruments and Measurements, 52(6), pp 1881-1885 """ import numpy as np +from ocbpy import ocb_time + def starkov_auroral_boundary(mlt, al=-1, bnd='ocb'): """Calculate the location of the auroral boundaries. @@ -121,3 +127,234 @@ def starkov_coefficient_values(al, coeff_name, bnd): log_al**2) + coeff_terms[coeff_name][bnd][3] * (log_al**3) return coeff + + +def gussenhoven_equatorward_auroral_boundary(mlt, kp=0, model='circle'): + """Calculate the location of the diffuse EAB. + + Parameters + ---------- + mlt : float or array-like + Magnetic local time in hours + kp : float or int + Decimal Kp geomagnetic index, ranges from 0 to 9 with "+" translating + to plus one third and "-" translating to minus one third (default=0) + model : str + Model flavour, with 'binned' using the value from the available hourly + solutions or NaN if that hour is absent, 'closest' to use the closest + hourly bin, or 'circle' to use the value of the circle fit to the + available hourly solutions (default='circle') + + Returns + ------- + bnd_lat : float or array-like + Location of the boundary in degrees away from the pole in corrected + geomagnetic coordinates for the specified magnetic local times. + + Raises + ------ + ValueError + If an unknown model method is requested. + + References + ---------- + [9]_ + + """ + closest = True if model.lower() == 'closest' else False + + # If desired, calculate the integer hour + if model.lower() in ['binned', 'closest']: + imlt = np.floor(mlt).astype(int) + + if imlt.shape == (): + imlt = np.array([imlt]) + else: + imlt = None + + # Get the desired latitude boundaries + colats, mlts = gussenhoven_colatitudes(kp, mlt_inds=imlt, closest=closest) + + # If desired, fit the co-latitude boundaries and return the locations + # at the exact MLT values + if model.lower() in ['binned', 'closest']: + bnd_lat = colats + elif model.lower() == "circle": + # Fit a circle to the boundaries at this Kp + phi_cent, r_cent, radius, _ = circle_fit(mlts, colats) + + # Calculate the circle location at the desired MLTs. Use the positive + # solution of the quadratic equation, as radius >= r_cent + theta = ocb_time.hr2rad(mlt) - phi_cent + bnd_lat = r_cent * np.cos(theta) + np.sqrt( + radius**2 - r_cent**2 * (np.sin(theta))**2) + else: + raise ValueError('unknown model method: {:}'.format(model)) + + return bnd_lat + + +def gussenhoven_colatitudes(kp, mlt_inds=None, closest=False): + """Calculate the Gussenhoven EAB model co-latitude values. + + Parameters + ---------- + kp : float or int + Decimal Kp geomagnetic index, ranges from 0 to 9 with "+" translating + to plus one third and "-" translating to minus one third + mlt_inds : array-like or NoneType + MLT hourly bins at which boundary locations will be returned or None + to return all available values (default=None) + closest : bool + If False, will return NaN for the request for an MLT value that does not + exist. If True, will find the closest binned value and return that + instead (default=False) + + Returns + ------- + colats : array-like + Co-latitude values in degrees away from the pole for the requested + MLT bins + mlts : array-like + MLT bin values for each co-latitude + + References + ---------- + [9]_ + + """ + # Set the function values for each MLT bin + offset = {0: 66.1, 1: 65.1, 4: 67.7, 5: 67.8, 6: 68.2, 7: 68.9, 8: 69.3, + 9: 69.5, 10: 69.6, 11: 70.1, 12: 69.4, 15: 70.9, 16: 71.6, + 17: 71.1, 18: 71.2, 19: 70.4, 20: 69.4, 21: 68.6, 22: 67.9, + 23: 67.8} + slope = {0: -1.99, 1: -1.55, 4: -1.48, 5: -1.87, 6: -1.90, 7: -1.91, + 8: -1.87, 9: -1.67, 10: -1.41, 11: -1.25, 12: -0.84, 15: -0.81, + 16: -1.28, 17: -1.31, 18: -1.74, 19: -1.83, 20: -1.89, 21: -1.86, + 22: -1.78, 23: -2.07} + + # Ensure the MLT indices are set + mlt_keys = np.array(list(offset.keys())) + if mlt_inds is None: + mlt_inds = mlt_keys + else: + mlt_inds = np.asarray(mlt_inds) + + # Initialize the output + colats = np.full(shape=mlt_inds.shape, fill_value=np.nan) + mlts = np.asarray(mlt_inds) + + # Cycle through each MLT index + for i, imlt in enumerate(mlt_inds): + if imlt in mlt_keys: + # Set the calculation key + colat_key = imlt + elif closest: + # Find the closest time with a boundary solution + iclose = abs(imlt - mlt_keys).argmin() + colat_key = mlt_keys[iclose] + mlts[i] = colat_key + else: + colat_key = None + + if colat_key is not None: + # Find the closest co-latitude solution + colat = 90 - (offset[colat_key] + slope[colat_key] * kp) + + # Ensure there are no values outside of the allowed range + if colat < 0.0: + colat = 0.0 + elif colat > 90.0: + colat = 90.0 + + # Save the output + colats[i] = colat + + return colats, mlts + + +def circle_fit(mlt, colat): + """Fit a circle to boundary estimates. + + Parameters + ---------- + mlt : array-like + Array of MLT values in hours for the boundary locations + colat : array-like + Array of magnetic co-latitude values in degrees denoting the boundary + locations, paired to the `mlt` inputs + + Returns + ------- + phi_cent : float + MLT location of the circle centers in radians + r_cent : float + Co-latitude of the circle center in degrees + radius : float + Radius of the circle in degrees latitude + r_err : float + RMS error of the fit + + Raises + ------ + ValueError + If the inputs are the wrong shape + + References + ---------- + Section D in [10]_ + + """ + # Ensure the arrays are of equal length + phi = ocb_time.hr2rad(np.asarray(mlt)) + colat = np.asarray(colat) + + if phi.shape != colat.shape: + raise ValueError("input MLT and MLat arrays must be the same shape") + + if len(phi.shape) != 1: + raise ValueError('routine only supports 1D arrays') + + # Get the cartesian values of the boundary points. Define the x-axis to + # lie along 0 MLT and the y-axis to lie along 6 MLT, with the origin at the + # magnetic pole. + x_bound = colat * np.cos(phi) + y_bound = colat * np.sin(phi) + + # Apply the circle fitting algorithm, a modified least-squares method, + # from Umbach and Jones section D + sum_x = np.sum(x_bound) + sum_y = np.sum(y_bound) + sum_x2 = np.sum(x_bound**2) + sum_y2 = np.sum(y_bound**2) + sum_xy = np.sum((x_bound * y_bound)) + sum_xy2 = np.sum(x_bound * (y_bound**2)) + sum_yx2 = np.sum(y_bound * (x_bound**2)) + sum_x3 = np.sum(x_bound * x_bound * x_bound) + sum_y3 = np.sum(y_bound * y_bound * y_bound) + + fit_a = (colat.shape[0] * sum_x2) - sum_x**2 + fit_b = (colat.shape[0] * sum_xy) - (sum_x * sum_y) + fit_c = (colat.shape[0] * sum_y2) - sum_y**2 + fit_d = 0.5 * ((colat.shape[0] * sum_xy2) - (sum_x * sum_y2) + + (colat.shape[0] * sum_x3) - (sum_x * sum_x2)) + fit_e = 0.5 * ((colat.shape[0] * sum_yx2) - (sum_y * sum_x2) + + (colat.shape[0] * sum_y3) - (sum_y * sum_y2)) + + # Determine the centre of the fitted circle in Cartesian coordinates + circle_denom = (fit_a * fit_c) - fit_b**2 + major_a = ((fit_d * fit_c) - (fit_b * fit_e)) / circle_denom + minor_b = ((fit_a * fit_e) - (fit_b * fit_d)) / circle_denom + + # Convert the circle centre to co-latitude in degrees + r_cent = np.sqrt(major_a**2 + minor_b**2) + phi_cent = np.arctan2(minor_b, major_a) + + # Determine the radius of the fitted circle + rad_bound = np.sqrt((x_bound - major_a)**2 + (y_bound - minor_b)**2) + radius = rad_bound.mean() + + # Calculate the error of the radius + r_err = np.sqrt(np.mean((rad_bound - radius)**2)) + + return phi_cent, r_cent, radius, r_err diff --git a/ocbpy/tests/test_models.py b/ocbpy/tests/test_models.py index a437d792..a34b9e6f 100644 --- a/ocbpy/tests/test_models.py +++ b/ocbpy/tests/test_models.py @@ -15,7 +15,7 @@ class TestStarkovModel(unittest.TestCase): """"Unit tests for the Starkov 1994 routines.""" def setUp(self): - """Initialize the test case by copying over necessary files.""" + """Initialize the test case by setting some values to test against.""" self.mlt = np.arange(0, 24, 1) self.coeff_out = { 'A0': {'ocb': -.07, 'eab': 1.16, 'diffuse': 3.44}, @@ -97,3 +97,199 @@ def test_bound_loc_float(self): self.assertLessEqual(lat, self.max_lat[bnd][ia]) return + + +class TestGussenhovelModel(unittest.TestCase): + """"Unit tests for the Gussenhoven 1983 routines.""" + + def setUp(self): + """Initialize the test case.""" + self.mlt = np.arange(0, 24, 1) + self.bad_mlt = [2, 3, 13, 14] + self.colat = {0: [23.9, 24.9, 22.3, 22.2, 21.8, 21.1, 20.7, 20.5, 20.4, + 19.9, 20.6, 19.1, 18.4, 18.9, 18.8, 19.6, 20.6, 21.4, + 22.1, 22.2], + 9: [41.81, 38.85, 35.62, 39.03, 38.90, 38.29, 37.53, + 35.53, 33.09, 31.15, 28.16, 26.39, 29.92, 30.69, + 34.46, 36.07, 37.61, 38.14, 38.12, 40.83]} + return + + def tearDown(self): + """Clean up the test environment.""" + del self.mlt, self.bad_mlt, self.colat + return + + def test_gussenhoven_colatitudes_good(self): + """Test the colat calculation for MLTs with solutions""" + + for kp in self.colat.keys(): + with self.subTest(kp=kp): + # Calculate the colatitude values with default kwargs + out_lat, _ = models.gussenhoven_colatitudes(kp) + + # Compare the output + self.assertEqual(len(out_lat), len(self.colat[kp])) + diff_lat = out_lat - np.asarray(self.colat[kp]) + self.assertLess(sum(diff_lat), 1.0e-3) + return + + def test_gussenhoven_colatitudes_bad(self): + """Test the colat calculation for MLTs without solutions.""" + # Cycle through each Kp + for kp in self.colat.keys(): + with self.subTest(kp=kp): + # Calculate the colatitude values with default kwargs + out_lat, out_mlt = models.gussenhoven_colatitudes( + kp, mlt_inds=self.bad_mlt) + + # Compare the output + self.assertTrue(np.isnan(out_lat).all()) + self.assertListEqual(list(out_mlt), self.bad_mlt) + return + + def test_gussenhoven_colatitudes_bad_close(self): + """Test the colat calculation for MLTs at nearest solutions.""" + # Reshape the MLT to exclude the bad MLT values + self.mlt = list(self.mlt) + for imlt in self.bad_mlt: + self.mlt.pop(self.mlt.index(imlt)) + + # Cycle through each Kp + for kp in self.colat.keys(): + with self.subTest(kp=kp): + # Calculate the colatitude values with default kwargs + out_lat, out_mlt = models.gussenhoven_colatitudes( + kp, mlt_inds=self.bad_mlt, closest=True) + + # Compare the output + self.assertFalse(np.isnan(out_lat).any()) + self.assertEqual(len(out_mlt), len(self.bad_mlt)) + + for i, imlt in enumerate(out_mlt): + self.assertFalse(imlt in self.bad_mlt) + self.assertAlmostEqual( + out_lat[i], self.colat[kp][self.mlt.index(imlt)]) + return + + def test_bad_model(self): + """Test a ValueError is raised for an unknown model type.""" + model = "not a model" + with self.assertRaisesRegex(ValueError, model): + models.gussenhoven_equatorward_auroral_boundary( + self.mlt, model=model) + return + + def test_gussenhoven_eab_binned(self): + """Test the EAB for MLTs with binned solutions.""" + # Cycle through each Kp + for kp in self.colat.keys(): + with self.subTest(kp=kp): + # Calculate the colatitude values with default kwargs + out_lat = models.gussenhoven_equatorward_auroral_boundary( + self.mlt, kp=kp, model='binned') + + # Compare the output + self.assertEqual(len(out_lat[~np.isnan(out_lat)]), + len(self.colat[kp])) + self.assertEqual(len(out_lat[np.isnan(out_lat)]), + len(self.bad_mlt)) + diff_lat = out_lat[~np.isnan(out_lat)] - np.asarray( + self.colat[kp]) + self.assertLess(sum(diff_lat), 1.0e-3) + return + + def test_gussenhoven_eab_closest(self): + """Test the EAB for MLTs with closest value solutions.""" + # Cycle through each Kp + for kp in self.colat.keys(): + with self.subTest(kp=kp): + # Calculate the colatitude values with default kwargs + out_lat = models.gussenhoven_equatorward_auroral_boundary( + self.mlt, kp=kp, model='closest') + + # Compare the output + self.assertEqual(len(out_lat), len(self.mlt)) + + j = 0 + for i, imlt in enumerate(self.mlt): + if imlt in self.bad_mlt: + try: + j -= 1 + self.assertAlmostEqual( + out_lat[i], self.colat[kp][j]) + j += 1 + except AssertionError: + j += 1 + self.assertAlmostEqual( + out_lat[i], self.colat[kp][j], + msg="failed at {:} MLT with j={:}".format( + imlt, j)) + else: + self.assertAlmostEqual( + out_lat[i], self.colat[kp][j], + msg="failed at {:} MLT with j={:}".format(imlt, j)) + j += 1 + return + + def test_gussenhoven_eab_circle(self): + """Test the EAB for MLTs with circle fit solutions.""" + # Cycle through each Kp + for kp in self.colat.keys(): + with self.subTest(kp=kp): + # Calculate the colatitude values with default kwargs + out_lat = models.gussenhoven_equatorward_auroral_boundary( + self.mlt, kp=kp, model='circle') + + # Compare the output + self.assertEqual(len(out_lat), len(self.mlt)) + + j = 0 + for i, imlt in enumerate(self.mlt): + if imlt in self.bad_mlt: + try: + j -= 1 + self.assertLessEqual( + abs(out_lat[i] - self.colat[kp][j]), 5.0) + j += 1 + except AssertionError: + j += 1 + self.assertLessEqual( + abs(out_lat[i] - self.colat[kp][j]), 5.0, + msg="failed at {:} MLT with j={:}".format( + imlt, j)) + else: + self.assertLessEqual( + abs(out_lat[i] - self.colat[kp][j]), 5.0, + msg="failed at {:} MLT with j={:}".format(imlt, j)) + j += 1 + return + + +class TestFits(unittest.TestCase): + """"Unit tests for the fitting routines.""" + + def setUp(self): + """Initialize the test case by setting some values to test against.""" + self.mlt = np.arange(0, 24, 1) + self.rvals = np.ones(shape=self.mlt.shape) + return + + def tearDown(self): + """Clean up the test environment.""" + del self.mlt, self.rvals + return + + def test_circle_fit(self): + """Test the circle fitting to a unit circle.""" + # Run the fitting routine + phi_cent, r_cent, radius, r_err = models.circle_fit( + self.mlt, self.rvals) + + # Test the output. The anglular offset can be any value with a + # radial offset of zero. It should be constrained within +/-pi + self.assertGreaterEqual(phi_cent, -np.pi) + self.assertLessEqual(phi_cent, np.pi) + self.assertAlmostEqual(r_cent, 0.0) + self.assertAlmostEqual(radius, 1.0) + self.assertAlmostEqual(r_err, 0.0) + return