diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md index 1537b560..c326d966 100644 --- a/.github/pull_request_template.md +++ b/.github/pull_request_template.md @@ -38,6 +38,7 @@ import ocbpy - [ ] Make sure you are merging into the ``develop`` (not ``main``) branch - [ ] My commits are formatted appropriately (following the SciPy/NumPy style) - [ ] My code follows the style guidelines of this project +- [ ] I assert that I have not used AI in the development of this pull request - [ ] I have performed a self-review of my own code - [ ] I have commented my code, particularly in hard-to-understand areas - [ ] I have made corresponding changes to the documentation diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 78e36ed9..c7a1634c 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -15,13 +15,13 @@ jobs: strategy: fail-fast: false matrix: - python-version: ["3.11"] + python-version: ["3.14"] name: Documentation tests steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v5 + uses: actions/setup-python@v6 with: python-version: ${{ matrix.python-version }} diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 36f417fc..328f70bd 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -13,15 +13,20 @@ jobs: fail-fast: false matrix: os: ["ubuntu-latest", "macos-latest", "windows-latest"] - python-version: ["3.10", "3.11", "3.12", "3.13"] - install-extras: ["base", "pysat_instruments", "dmsp_ssj"] + python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"] + install-extras: ["base", "dmsp_ssj"] + include: + # Pysat test for only pysat-passing Python version + - python-version: "3.10" + os: "ubuntu-latest" + install-extras: "pysat_instruments" name: Python ${{ matrix.python-version }} on ${{ matrix.os }} with ${{ matrix.install-extras }} runs-on: ${{ matrix.os }} steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v5 + uses: actions/setup-python@v6 with: python-version: ${{ matrix.python-version }} @@ -50,33 +55,13 @@ 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: + flag-name: run=${{ join(matrix.*, '-') }} + parallel: true + format: cobertura + debug: true finish: name: Finish Coverage Analysis @@ -84,9 +69,7 @@ jobs: runs-on: "ubuntu-latest" steps: - name: Coveralls Finished - env: - COVERALLS_REPO_TOKEN: ${{ secrets.COVERALLS_REPO_TOKEN }} - COVERALLS_PARALLEL: true - run: | - curl -sL https://coveralls.io/coveralls-linux.tar.gz | tar -xz - ./coveralls done --build-number ${{ github.run_number }} + uses: coverallsapp/github-action@v2 + with: + parallel-finished: true + carryforward: "all" diff --git a/.github/workflows/pip_rc_install.yml b/.github/workflows/pip_rc_install.yml index 1124b9ba..70aa2e53 100644 --- a/.github/workflows/pip_rc_install.yml +++ b/.github/workflows/pip_rc_install.yml @@ -13,21 +13,21 @@ jobs: fail-fast: false matrix: os: ["ubuntu-latest", "macos-latest", "windows-latest"] - python-version: ["3.11"] # Keep this version at the highest supported Python version + python-version: ["3.14"] # Keep this version at the highest supported Python version name: Python ${{ matrix.python-version }} on ${{ matrix.os }} runs-on: ${{ matrix.os }} steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 - name: Set up Python ${{ matrix.python-version }} - uses: actions/setup-python@v5 + uses: actions/setup-python@v6 with: python-version: ${{ matrix.python-version }} - name: Install standard dependencies run: pip install -r requirements.txt - - name: Install pysat RC + - name: Install ocbpy RC run: pip install --no-deps --pre -i https://test.pypi.org/simple/ --extra-index-url https://pypi.org/simple/ ocbpy - name: Check that installation imports correctly diff --git a/CONTRIBUTING.rst b/CONTRIBUTING.rst index 5d04921d..7c3caa71 100644 --- a/CONTRIBUTING.rst +++ b/CONTRIBUTING.rst @@ -88,3 +88,17 @@ To run a subset of tests from the test directory for a specific environment:: To run all the tests for a specific environment:: python -m unittest discover + +Use of 'AI' +----------- + +This project is human centered and developers do not use Artificial Inteligence +(AI) tools such as large language models in its creation or maintenance. +The purpose of this tool is to improve high-latitude science, and thus values +interactions with scientists and scientific programmers. Code generated by a +large language model or similar technology, such as Anthropic’s Claude, +GitHub/Microsoft’s Copilot, OpenAI’s ChatGPT, Facebook/Meta’s Code Llama et al., +is not compliant with the covenants and representations of OCBpy’s Contributor’s +Agreement, and is thus not acceptable as code for OCBpy. + + diff --git a/Changelog.rst b/Changelog.rst index 09b9bbb5..0c532ad6 100644 --- a/Changelog.rst +++ b/Changelog.rst @@ -3,6 +3,17 @@ Changelog Summary of all changes made since the first stable release +0.7.0 (05-29-2026) +------------------ +* ENH: Added the Gussenhoven (1983) model for the EAB +* ENH: Added the CH-Aurora-2014 model for the OCB and EAB +* ENH: Added an AI usage policy (no AI use allowed) +* DEP: Removed support for older versions of `zenodo_get`, updated download + checks to account for breaking changes made in the newer versions +* MAINT: Added support for Python 3.14 +* MAINT: Updated GitHub CI actions +* MAINT: Updated numpy usage + 0.6.0 (07-07-2025) ------------------ * ENH: Updated the AMPERE boundaries to include 2022-2024, inclusive diff --git a/README.md b/README.md index 1abf5c72..4e767f25 100644 --- a/README.md +++ b/README.md @@ -45,7 +45,7 @@ These routines may be used as a guide to write routines for other datasets. # Python versions -This module currently supports Python version 3.10 - 3.13. +This module currently supports Python version 3.10 - 3.14. # Dependencies @@ -53,7 +53,7 @@ The listed dependecies were tested with the following versions: * numpy * aacgmv2 * pysat (3.2.1+) - * zenodo_get + * zenodo_get (2.0.0+) Testing is performed using the python module, unittest. To limit dependency issues, the pysat and zenodo_get dependencies are optional. diff --git a/docs/citing.rst b/docs/citing.rst index fee63c5d..cb4db0b0 100644 --- a/docs/citing.rst +++ b/docs/citing.rst @@ -98,7 +98,7 @@ boundary method and data set are provided below. * **SSJ Auroral Boundaries (2010-2014)**: Kilcommons, L., et al. (2019). Defense Meteorology Satellite Program (DMSP) Electron Precipitation (SSJ) Auroral Boundaries, 2010-2014 (Version 1.0.0) [Data set]. Zenodo. - http://doi.org/10.5281/zenodo.3373812 + doi:10.5281/zenodo.3373812 .. _cite-starkov: @@ -112,3 +112,37 @@ 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. + + +.. _cite-xiong: + +CH-Aurora-2014 Model +-------------------- + +The CH-Aurora-2014 auroral boundary model provides either a mathematical +description of the poleward and equatorward auroral boundaries in each +hemisphere or an assimilated location based on FAC observations of the auroral +boundaries in that hemisphere. The data that may be used to adjust the model +are described in the first article below, while the model itself is described in +the second article. If you use this model please cite both it and your Newell +Couping Function data source when publishing your work. + +* Xiong, C., et al. (2014) Determining the boundaries of the auroral oval from + CHAMP field-aligned current signatures - Part 1, Ann. Geophys., 32, + pp 609-622, doi:10.5194/angeo-32-609-2014. +* Xiong, C. and H. Luhr (2014) An empirical model of the auroral oval derived + from CHAMP field-aligned current signatures - Part 2, Ann. Geophys., 32, + pp 623-631, doi:10.5194/angeo-32-623-2014 diff --git a/docs/examples/ex_dmsp.rst b/docs/examples/ex_dmsp.rst index 069037ad..8ad12287 100644 --- a/docs/examples/ex_dmsp.rst +++ b/docs/examples/ex_dmsp.rst @@ -13,13 +13,12 @@ need the `ssj_auroral_boundary `__ package, but now we preferentially support using the `zendodo_get `__ package to obtain the -boundary files from their -`archive `__. Once installed, -you can download DMSP SSJ data and obtain a boundary file for a specified time -period (or all available times) using :py:mod:`ocbpy.boundaries.dmsp_ssj_files`. -For this example, we'll use a single day. You can download the files into any -directory, but this example will put them in the same directory as the other -boundary files. +boundary files from their archive at the website: zenodo.org/record/3373812. +Once installed, you can download DMSP SSJ data and obtain a boundary file for a +specified time period (or all available times) using +:py:mod:`ocbpy.boundaries.dmsp_ssj_files`. For this example, we'll use a single +day. You can download the files into any directory, but this example will put +them in the same directory as the other boundary files. :: 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..c89fa6fb 100644 --- a/docs/ocb_models.rst +++ b/docs/ocb_models.rst @@ -27,6 +27,32 @@ 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-champ: + +CH-Aurora-2014 +-------------- + +The CH-Aurora-2014 model (see :ref:`cite-xiong`) uses a mathematical +formulation based on CHAMP field-aligned current (FAC) data and a delayed +time-history of the Newell Coupling Function. They specify both auroral +boundaries in each hemisphere. Although this model may be driven entirely by +the Newell Coupling Function, the authors recommend adjusting the model output +with CHAMP measurements of the FAC boundaries. + + .. _bound-model-module: Boundary Models Module diff --git a/ocbpy/boundaries/dmsp_ssj_files.py b/ocbpy/boundaries/dmsp_ssj_files.py index 74279af7..52ca5299 100644 --- a/ocbpy/boundaries/dmsp_ssj_files.py +++ b/ocbpy/boundaries/dmsp_ssj_files.py @@ -20,15 +20,14 @@ .. [7] Kilcommons, L., Redmon, R., & Knipp, D. (2019). Defense Meteorology Satellite Program (DMSP) Electron Precipitation (SSJ) Auroral Boundaries, - 2010-2014 (1.0.0) [Data set]. Zenodo. https://doi.org/10.5281/zenodo.3373812 + 2010-2014 (1.0.0) [Data set]. Zenodo. doi:10.5281/zenodo.3373812 """ import datetime as dt -from io import StringIO +import hashlib import numpy as np import os -import sys import zipfile import aacgmv2 @@ -108,50 +107,42 @@ def fetch_ssj_boundary_files(stime=None, etime=None, out_dir=None, if not os.path.isdir(out_dir): raise ValueError("can't find the output directory") - # Download the zenodo archive, capturing the output - zenodo_io = StringIO() - sys.stdout = zenodo_io - sys.stderr = zenodo_io + # Download the zenodo archive + zenodo_get.download(doi=doi, output_dir=out_dir) - # TODO(#151): remove the old (second) way of calling zenodo_get - if hasattr(zenodo_get, "download"): - zenodo_get.download(doi=doi, output_dir=out_dir) - zenodo_checksum = None - else: - zenodo_get.zenodo_get([doi, '-o', out_dir]) - zenodo_checksum = os.path.join(out_dir, 'md5sums.txt') + # Get the checksum file (does not download) + zenodo_get.download(doi=doi, output_dir=out_dir, md5=True) + zenodo_checksum = os.path.join(out_dir, "md5sums.txt") - # Parse the output and retrieve files from the zip archive - sys.stdout = sys.__stdout__ - sys.stderr = sys.__stderr__ - zen_msg = zenodo_io.getvalue() - zen_split = zen_msg.split() + # Verify the checksum and retrieve the zip archive name + if os.path.isfile(zenodo_checksum): + with open(zenodo_checksum, 'r') as zcheck: + csum, zip_name = zcheck.read().split() - if zen_msg.find('Checksum is correct') < 0 and zen_msg.find( - 'already downloaded correctly') < 0: - raise IOError('Bad checksum: {:s}'.format(zen_msg)) + # Set the archive name + archive_name = os.path.join(out_dir, zip_name) - # Remove the checksum file if the download problem wasn't found there - if zenodo_checksum is not None: - os.remove(zenodo_checksum) + if not os.path.isfile(archive_name): + raise IOError( + 'error downloading archive to output dir: {:}'.format( + archive_name)) - # Get the archive name from the output - try: - link_ind = zen_split.index('Link:') + 1 + # Get the archive checkshum + with open(os.path.join(out_dir, archive_name), 'rb') as afile: + asum = hashlib.md5(afile.read()).hexdigest() - # If the archive is already available, message may differ - zip_name = os.path.split(zen_split[link_ind]) - while zip_name[-1].find(".zip") <= 0: - zip_name = os.path.split(zip_name[0]) + zen_msg = "Checksum is correct" if asum == csum else "Corrupted data" + else: + zen_msg = "Checksum file not created" - # Set the archive name - archive_name = os.path.join(out_dir, zip_name[-1]) - except (ValueError, IndexError): - raise IOError('unable to identify zenodo archive: {:}'.format(zen_msg)) + # Remove the checksum file + if zenodo_checksum is not None: + os.remove(zenodo_checksum) - if not os.path.isfile(archive_name): - raise IOError('error downloading archive to output dir: {:}'.format( - archive_name)) + # Evaluate checksum output + if zen_msg.find('Checksum is correct') < 0 and zen_msg.find( + 'already downloaded correctly') < 0: + raise IOError('Bad checksum: {:s}'.format(zen_msg)) # Access the zip archive with zipfile.ZipFile(archive_name, 'r') as zref: diff --git a/ocbpy/boundaries/models.py b/ocbpy/boundaries/models.py index 4c36d7fa..0dec4e70 100644 --- a/ocbpy/boundaries/models.py +++ b/ocbpy/boundaries/models.py @@ -9,14 +9,23 @@ ---------- .. [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 +.. [11] Xiong and Luhr (2014) An empirical model of the auroral oval derived + from CHAMP field-aligned current signatures - Part 2, Ann. Geophys., 32, + pp 623-631, doi:10.5194/angeo-32-623-2014 """ 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. + """Calculate the location of the Starkov auroral boundaries. Parameters ---------- @@ -121,3 +130,450 @@ 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 + + +def ch_aurora_2014_boundary(mlt, em=0, bnd='ocb', hemi=1, obs_colat=None, + obs_mlt=None): + """Calculate the location of the CH-Aurora-2014 auroral boundaries. + + Parameters + ---------- + mlt : float or array-like + Magnetic local time in hours + em : float or int + The time-integrated Newell Coupling Function, as described in Equation 2 + of [11]_ (default=0) + bnd : str + Boundary to calculate, expects one of 'ocb' or 'eab' (default='ocb') + hemi : int + Sign denoting the hemisphere, 1 for North and -1 for South (default=1) + obs_colat : array-like + CHAMP, or other, boundary observation co-latitudes (default=None) + obs_mlt : array-like + MLT of the observed boundary locations (default=None) + + Returns + ------- + bnd_lat : float or array-like + Location of the boundary in degrees away from the pole in apex + geomagnetic coordinates for the specified magnetic local times. + + References + ---------- + [11]_ + + """ + # Calculate each parameter for the desired IMF conditions; all are in + # degrees except for phi0, which is in radians. + semix, semiy, x0, y0, phi0 = ch_aurora_2014_coefficient_values(em, bnd, + hemi) + + # If observations are supplied, ensure their MLT are included in the calc + calc_mlt = np.array(mlt) + iobs = 0 + if obs_mlt is not None: + obs_mlt_array = np.array(obs_mlt) + + if calc_mlt.ndim == 0: + iobs = 1 + if obs_mlt_array.ndim == 0: + calc_mlt = np.array([mlt, obs_mlt]) + else: + calc_mlt = list(obs_mlt) + calc_mlt.insert(0, mlt) + calc_mlt = np.array(calc_mlt) + else: + calc_mlt = list(mlt) + iobs = len(mlt) + if obs_mlt_array.ndim == 0: + calc_mlt.append(obs_mlt) + else: + calc_mlt.extend(list(obs_mlt)) + calc_mlt = np.array(calc_mlt) + + # Calculate the angular MLT + ang_lt = ocb_time.hr2rad(calc_mlt) + + # First round of calculations with phi == ang_lt + bnd_lat = ch_aurora_2014_radius(ang_lt, semix, semiy, x0, y0, phi0) + + # Now calculate the angular time from the radius + xprime = bnd_lat * np.cos(ang_lt + phi0) + x0 + yprime = bnd_lat * np.sin(ang_lt + phi0) + y0 + ang_prime = np.arctan2(yprime, xprime) + + # Recalculate the radius using the new local times + bnd_lat = ch_aurora_2014_radius(ang_prime, semix, semiy, x0, y0, phi0) + + # If desired, modify the boundary locations based on observations + if iobs > 0: + # Get the model values for the observation times + mod_colat = bnd_lat[iobs:] + del_lat = np.nanmean(obs_colat - mod_colat) + + # Insert the adjustment + bnd_lat = ch_aurora_2014_radius(ang_lt, semix, semiy, x0, y0, phi0, + del_lat) + + # Now calculate the angular time from the radius + xprime = bnd_lat * np.cos(ang_lt + phi0) + x0 + yprime = bnd_lat * np.sin(ang_lt + phi0) + y0 + ang_prime = np.arctan2(yprime, xprime) + + # Recalculate the radius using the new local times + bnd_lat = ch_aurora_2014_radius(ang_prime, semix, semiy, x0, y0, phi0, + del_lat) + + # Downselect the output to include only the desired MLT + bnd_lat = bnd_lat[:iobs] + + # Ensure the output co-latitude is a float if the input MLT is a float + if calc_mlt.ndim == 0: + if bnd_lat.ndim == 0: + bnd_lat = float(bnd_lat) + else: + bnd_lat = float(bnd_lat[0]) + + return bnd_lat + + +def ch_aurora_2014_coefficient_values(em, bnd, hemi): + """Retrive the CH-Aurora-2014 auroral boundary coefficients. + + Parameters + ---------- + em : float or int + The time-integrated Newell Coupling Function, as described in Equation 2 + of [11]_ (default=0) + bnd : str + Boundary to calculate, expects one of 'ocb' or 'eab' + hemi : int + Sign denoting the hemisphere, 1 for North and -1 for South + + Returns + ------- + semix : float + Semi-axis along the x-plane + semiy : float + Semi-axis along the y-plane + x0 : float + x-coordinate of the ellipse center + y0 : float + y-coordinate of the ellipse center + phi0 : float + Orientation angle of the ellipse + + References + ---------- + Tables 2 and 3 in [11]_ + + """ + # Set the ellipse parameters from Tables 2 and 3 from Xiong and Luhr. + # The list includes the p0, p1, and p2 parameters in that order, all + # provided in degrees + semix_param = {'eab': {1: [18.861, 0.95470, -0.023836], + -1: [18.559, 0.81597, -0.013209]}, + 'ocb': {1: [12.813, 0.26173, -5.9729e-4], + -1: [13.251, 0.28870, -3.0559e-4]}} + semiy_param = {'eab': {1: [20.562, 1.1504, -0.029566], + -1: [19.549, 0.10752, -0.024605]}, + 'ocb': {1: [9.5486, 0.96759, -0.029556], + -1: [11.605, 0.82006, -0.024073]}} + x0_param = {'eab': {1: [4.1263, 0.027827, 0.0], + -1: [3.6946, 0.044667, 0.0]}, + 'ocb': {1: [4.5175, -0.025310, 0.0], + -1: [4.2526, -0.021674, 0.0]}} + y0_param = {'eab': {1: [-0.32637, -0.028855, 1.6569e-3], + -1: [-0.60436, -0.011623, 6.5985e-4]}, + 'ocb': {1: [-0.39316, -0.15732, 5.613e-3], + -1: [-1.1330, -0.024479, 7.0729e-4]}} + phi0_param = {'eab': {1: [-3.1555, 0.31147, 0.0], + -1: [-8.8836, -0.10934, 0.0]}, + 'ocb': {1: [-8.5358, 1.2831, 0.0], + -1: [3.7050, -1.5508, 0.0]}} + + # Calculate each parameter for the desired IMF conditions + semix = semix_param[bnd][hemi][0] + semix_param[bnd][hemi][ + 1] * em + semix_param[bnd][hemi][2] * em * em + semiy = semiy_param[bnd][hemi][0] + semiy_param[bnd][hemi][ + 1] * em + semiy_param[bnd][hemi][2] * em * em + x0 = x0_param[bnd][hemi][0] + x0_param[bnd][hemi][1] * em + x0_param[bnd][ + hemi][2] * em * em + y0 = y0_param[bnd][hemi][0] + y0_param[bnd][hemi][1] * em + y0_param[bnd][ + hemi][2] * em * em + phi0 = np.radians(phi0_param[bnd][hemi][0] + phi0_param[bnd][hemi][1] * em + + phi0_param[bnd][hemi][2] * em * em) + + return semix, semiy, x0, y0, phi0 + + +def ch_aurora_2014_radius(ang_lt, semix, semiy, x0, y0, phi0, del_rad=0.0): + """Calculate the CH-Aurora-2014 ellipse radii at the desired times. + + Parameters + ---------- + ang_lt : float or array-like + Angular measure of MLT in radians + semix : float + Semi-axis along the x-plane in degrees + semiy : float + Semi-axis along the y-plane in degrees + x0 : float + x-coordinate of the ellipse center in degrees + y0 : float + y-coordinate of the ellipse center in degrees + phi0 : float + Orientation angle of the ellipse in radians + del_rad : float + Assimilative adjustment to the radius in degrees + + Returns + ------- + rad : float or array-like + Radius of the ellipse at the desired MLT in degrees + + References + ---------- + Equations 3 and 4 in [11]_ + + """ + # Calculate the radius from the centre of the ellipse + r0 = del_rad + semix * semiy / np.sqrt((semix * np.sin(ang_lt + phi0))**2 + + (semiy * np.cos(ang_lt + phi0))**2) + + # Adjust to provide the radius from the magnetic pole + rad = np.sqrt((r0 * np.cos(ang_lt + phi0) + x0)**2 + + (r0 * np.sin(ang_lt + phi0) + y0)**2) + + return rad diff --git a/ocbpy/ocb_scaling.py b/ocbpy/ocb_scaling.py index 0378157a..965a7f32 100644 --- a/ocbpy/ocb_scaling.py +++ b/ocbpy/ocb_scaling.py @@ -473,7 +473,7 @@ def vect_mag(self, vect_mag): self._vect_mag = vect_sqrt else: if np.any(np.greater(abs(vect_mag - vect_sqrt), 1.0e-3, - where=~np.isnan(vect_mag))): + where=~np.isnan(vect_mag), out=None)): ocbpy.logger.warning("".join([ "inconsistent vector components with a maximum difference ", "of {:} > 1.0e-3".format(abs(vect_mag - vect_sqrt).max())])) @@ -1221,7 +1221,7 @@ def archav(hav): alpha = np.full(shape=hav.shape, fill_value=np.nan) # If the number is positive, calculate the angle - norm_mask = (np.greater_equal(hav, 1.0e-16, where=~np.isnan(hav)) + norm_mask = (np.greater_equal(hav, 1.0e-16, where=~np.isnan(hav), out=None) & ~np.isnan(hav)) if np.any(norm_mask): if hav.shape == (): @@ -1231,7 +1231,7 @@ def archav(hav): # The number is small enough that machine precision may have changed # the sign, but it's a single-precission zero - small_mask = (np.less(abs(hav), 1.0e-16, where=~np.isnan(hav)) + small_mask = (np.less(abs(hav), 1.0e-16, where=~np.isnan(hav), out=None) & ~np.isnan(hav)) if np.any(small_mask): if hav.shape == (): diff --git a/ocbpy/ocb_time.py b/ocbpy/ocb_time.py index db915b27..9d0b46b7 100644 --- a/ocbpy/ocb_time.py +++ b/ocbpy/ocb_time.py @@ -434,19 +434,19 @@ def fix_range(values, min_val, max_val, val_range=None): # Fix the values, allowing for deviations that are multiples of the # value range. Also propagate NaNs - ibad = (np.greater_equal(fixed_vals, max_val, where=~np.isnan(fixed_vals)) - & ~np.isnan(fixed_vals)) + ibad = np.greater_equal(fixed_vals, max_val, where=~np.isnan(fixed_vals), + out=None) & ~np.isnan(fixed_vals) while np.any(ibad): fixed_vals[ibad] -= val_range ibad = (np.greater_equal(fixed_vals, max_val, - where=~np.isnan(fixed_vals)) + where=~np.isnan(fixed_vals), out=None) & ~np.isnan(fixed_vals)) - ibad = (np.less(fixed_vals, min_val, where=~np.isnan(fixed_vals)) + ibad = (np.less(fixed_vals, min_val, where=~np.isnan(fixed_vals), out=None) & ~np.isnan(fixed_vals)) while np.any(ibad): fixed_vals[ibad] += val_range - ibad = (np.less(fixed_vals, min_val, where=~np.isnan(fixed_vals)) - & ~np.isnan(fixed_vals)) + ibad = np.less(fixed_vals, min_val, where=~np.isnan(fixed_vals), + out=None) & ~np.isnan(fixed_vals) return fixed_vals diff --git a/ocbpy/tests/test_boundary_dual.py b/ocbpy/tests/test_boundary_dual.py index 5982e5ae..b64c8da3 100644 --- a/ocbpy/tests/test_boundary_dual.py +++ b/ocbpy/tests/test_boundary_dual.py @@ -586,7 +586,7 @@ def eval_coords(self, hemisphere=1, tol=1.0e-7, ind=None, revert=False, len(numpy.isnan(rlat))) if not numpy.isnan(rlat).all(): self.assertTrue( - numpy.less(abs(self.out[0] - rlat), tol, + numpy.less(abs(self.out[0] - rlat), tol, out=None, where=~numpy.isnan(rlat)).all(), msg="unequal {:s}: {:} != {:}".format(lat_str, self.out[0], rlat)) @@ -594,7 +594,7 @@ def eval_coords(self, hemisphere=1, tol=1.0e-7, ind=None, revert=False, len(numpy.isnan(rmlt))) if not numpy.isnan(rmlt).all(): self.assertTrue( - numpy.less(abs(self.out[1] - rmlt), tol, + numpy.less(abs(self.out[1] - rmlt), tol, out=None, where=~numpy.isnan(rmlt)).all(), msg="unequal {:s}: {:} != {:}".format(mlt_str, self.out[1], rmlt)) @@ -605,7 +605,8 @@ def eval_coords(self, hemisphere=1, tol=1.0e-7, ind=None, revert=False, if not numpy.isnan(self.out[2]).all(): self.assertTrue( numpy.less(abs(self.out[2] - self.olat[hemisphere]), - tol, where=~numpy.isnan(self.out[2])).all(), + tol, where=~numpy.isnan(self.out[2]), + out=None).all(), msg="unequal OCB latitude: {:} != {:}".format( self.out[2], self.olat[hemisphere])) @@ -615,7 +616,8 @@ def eval_coords(self, hemisphere=1, tol=1.0e-7, ind=None, revert=False, else: self.assertTrue( numpy.less(abs(self.out[3] - self.rcorr), tol, - where=~numpy.isnan(self.out[3])).all(), + where=~numpy.isnan(self.out[3]), + out=None).all(), msg="unequal radial correction: {:} != {:}".format( self.out[3], self.rcorr)) else: @@ -1119,7 +1121,7 @@ def test_get_current_aacgm_boundary_set(self): # Ensure the expected values and fill values are returned self.assertTrue( - numpy.less(abs(abound - self.bounds[i]), 1e-7, + numpy.less(abs(abound - self.bounds[i]), 1e-7, out=None, where=~numpy.isnan(abound)).all(), msg="unexpected boundary: {:} != {:}".format( abound, self.bounds[i])) diff --git a/ocbpy/tests/test_boundary_ocb.py b/ocbpy/tests/test_boundary_ocb.py index 467ed8c1..993b6da3 100644 --- a/ocbpy/tests/test_boundary_ocb.py +++ b/ocbpy/tests/test_boundary_ocb.py @@ -590,11 +590,11 @@ def test_normal_coord_north_array(self): self.out = self.ocb.normal_coord(self.lat, self.mlt) self.assertTrue(numpy.all(numpy.less(abs(self.out[0] - self.ocb_lat), - 1.0e-7, + 1.0e-7, out=None, where=~numpy.isnan(self.out[0])) | numpy.isnan(self.out[0]))) self.assertTrue(numpy.all(numpy.less(abs(self.out[1] - self.ocb_mlt), - 1.0e-7, + 1.0e-7, out=None, where=(~numpy.isnan(self.out[1]))) | numpy.isnan(self.out[1]))) self.assertTrue(numpy.where(numpy.isnan(self.out[0])) @@ -679,11 +679,11 @@ def test_revert_coord_north_array(self): self.r_corr) self.assertTrue(numpy.all(numpy.less(abs(self.out[0] - self.lat), - 1.0e-7, + 1.0e-7, out=None, where=~numpy.isnan(self.out[0])) | (numpy.isnan(self.out[0])))) self.assertTrue(numpy.all(numpy.less(abs(self.out[1] - self.mlt), - 1.0e-7, + 1.0e-7, out=None, where=(~numpy.isnan(self.out[1]) & (self.lat < 90.0))) | numpy.isnan(self.out[0]) @@ -984,10 +984,10 @@ def test_normal_coord_south_array(self): self.out = self.ocb.normal_coord(self.lat, self.mlt) self.assertTrue(numpy.all(numpy.less(abs(self.out[0] - self.ocb_lat), - 1.0e-7, + 1.0e-7, out=None, where=~numpy.isnan(self.out[0])))) self.assertTrue(numpy.all(numpy.less(abs(self.out[1] - self.ocb_mlt), - 1.0e-7, + 1.0e-7, out=None, where=~numpy.isnan(self.out[1])))) self.assertTrue(numpy.all(numpy.where(numpy.isnan(self.out[0]))[0] == numpy.where(numpy.isnan(self.ocb_lat))[0])) @@ -1041,11 +1041,11 @@ def test_revert_coord_south_array(self): self.r_corr) self.assertTrue(numpy.all(numpy.less(abs(self.out[0] - self.lat), - 1.0e-7, + 1.0e-7, out=None, where=~numpy.isnan(self.out[0])) | numpy.isnan(self.out[0]))) self.assertTrue(numpy.all(numpy.less(abs(self.out[1] - self.mlt), - 1.0e-7, + 1.0e-7, out=None, where=(~numpy.isnan(self.out[1]) & (self.lat > -90.0))) | numpy.isnan(self.out[0]) diff --git a/ocbpy/tests/test_models.py b/ocbpy/tests/test_models.py index a437d792..a0433a3e 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,347 @@ 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 + + +class TestCHAMPModel(unittest.TestCase): + """"Unit tests for the CH-Aurora-2014 routines.""" + + def setUp(self): + """Initialize the test case by setting some values to test against.""" + self.mlt = np.arange(0, 24, 1) + self.coeff_out = {'semix': {'ocb': {1: 12.813, -1: 13.251}, + 'eab': {1: 18.861, -1: 18.559}}, + 'semiy': {'ocb': {1: 9.5486, -1: 11.605}, + 'eab': {1: 20.562, -1: 19.549}}, + 'x0': {'ocb': {1: 4.5175, -1: 4.2526}, + 'eab': {1: 4.1263, -1: 3.6946}}, + 'y0': {'ocb': {1: -0.39316, -1: -1.1330}, + 'eab': {1: -0.32637, -1: -0.60436}}, + 'phi0': {'ocb': {1: -0.1489778, -1: 0.0646644}, + 'eab': {1: -0.055074, -1: -0.155048}}} + self.em = [0, 10.5] + self.iobs = [0, 12] + self.obs_mlt = self.mlt[self.iobs] + self.obs_colat = {'ocb': np.array([18.0, 8.0]), + 'eab': np.array([25.0, 16.0])} + self.max_lat = {'ocb': {1: [17.31668454, 20.1346482], + -1: [17.5881323, 20.741303775]}, + 'eab': {1: [23.052935, 31.313414], + -1: [22.354591, 29.69814346]}} + return + + def tearDown(self): + """Clean up the test environment.""" + del self.mlt, self.coeff_out, self.em, self.max_lat, self.obs_mlt + del self.obs_colat, self.iobs + return + + def test_coeff_construction(self): + """Test coefficient calculation for an Em of 0.""" + + for coeff in self.coeff_out.keys(): + for bnd in self.coeff_out[coeff].keys(): + for hemi in [1, -1]: + with self.subTest(coeff=coeff, bnd=bnd, hemi=hemi): + # Calculate the coefficient value + out = models.ch_aurora_2014_coefficient_values( + self.em[0], bnd, hemi) + + # Compare the output + for ic, coeff in enumerate(['semix', 'semiy', 'x0', + 'y0']): + self.assertEqual( + out[ic], self.coeff_out[coeff][bnd][hemi], + msg="{:s} does not match".format(coeff)) + + # Phi0 has been converted to radians, so equality + # will not be exact. Use significance from paper + self.assertAlmostEqual( + out[-1], self.coeff_out['phi0'][bnd][hemi], + places=5, msg="phi0 does not match") + return + + def test_coeff_bad_hemi(self): + """Test a KeyError is raised for an unknown hemisphere.""" + hemi = "north" + with self.assertRaisesRegex(KeyError, hemi): + models.ch_aurora_2014_coefficient_values( + self.em[0], list(self.max_lat.keys())[0], hemi) + return + + def test_coeff_bad_bnd(self): + """Test a KeyError is raised for an unknown boundary name.""" + bound = "not a boundary" + with self.assertRaisesRegex(KeyError, bound): + models.ch_aurora_2014_coefficient_values(self.em[0], bound, 1) + return + + def test_bound_loc_array(self): + """Test the expected boundary location across an MLT array.""" + # Get the boundary keys + bnds = list(self.max_lat.keys()) + + # Cycle through low and high Em values + for ie, in_em in enumerate(self.em): + for hemi in self.max_lat[bnds[0]].keys(): + with self.subTest(em=in_em, hemi=hemi): + lats = { + bnd: models.ch_aurora_2014_boundary( + self.mlt, em=in_em, bnd=bnd, hemi=hemi) + for bnd in bnds} + + # Test the output latitude shape and values + for bnd in bnds: + self.assertTupleEqual(self.mlt.shape, lats[bnd].shape) + self.assertAlmostEqual(max(lats[bnd]), + self.max_lat[bnd][hemi][ie], + places=5) + + # Test that the OCB is greater than zero and the EAB is + # greater than the OCB + self.assertGreaterEqual(min(lats['ocb']), 0) + self.assertTrue(np.all(lats['eab'] > lats['ocb'])) + + return + + def test_bound_loc_float(self): + """Test the expected boundary location across an MLT value.""" + # Cycle through low and high Em values + for ie, in_em in enumerate(self.em): + for bnd in self.max_lat.keys(): + for hemi in self.max_lat[bnd].keys(): + with self.subTest(em=in_em, bnd=bnd, hemi=hemi): + lat = models.ch_aurora_2014_boundary( + self.mlt[0], em=in_em, bnd=bnd, hemi=hemi) + + # Test the output latitude shape and values + self.assertTrue(isinstance(lat, float)) + self.assertGreaterEqual(lat, 0) + self.assertLessEqual(lat, self.max_lat[bnd][hemi][ie]) + + return + + def test_bound_assim(self): + """Test the expected assimilated boundary location.""" + + # Cycle through low and high Em values + for ie, in_em in enumerate(self.em): + for bnd in self.max_lat.keys(): + for hemi in self.max_lat[bnd].keys(): + with self.subTest(em=in_em, bnd=bnd, hemi=hemi): + mlat = models.ch_aurora_2014_boundary( + self.mlt, em=in_em, bnd=bnd, hemi=hemi) + alat = models.ch_aurora_2014_boundary( + self.mlt, em=in_em, bnd=bnd, hemi=hemi, + obs_mlt=self.obs_mlt, obs_colat=self.obs_colat[bnd]) + + # Test the output latitude shape + self.assertTupleEqual(self.mlt.shape, alat.shape) + + # Model and assimilation should differ + self.assertGreater(abs(alat - mlat).min(), 1.0e-2) + + # Assimilated output should be closer to the provided + # data points than the original model in at least one + # of the assimilated locations + self.assertTrue( + np.any(abs(alat[self.iobs] - self.obs_colat[bnd]) + < abs(mlat[self.iobs] + - self.obs_colat[bnd]))) + return diff --git a/ocbpy/tests/test_ocb_scaling.py b/ocbpy/tests/test_ocb_scaling.py index f2936023..e326ff97 100644 --- a/ocbpy/tests/test_ocb_scaling.py +++ b/ocbpy/tests/test_ocb_scaling.py @@ -540,10 +540,10 @@ def test_update_loc_coords_float(self): # Evaluate the output self.assertAlmostEqual( - float(self.vdata.lat), mag_out[coord][0], places=4, + self.vdata.lat.item(), mag_out[coord][0], places=4, msg="unexpected magnetic latitude") self.assertAlmostEqual( - float(self.vdata.lt), mag_out[coord][1], places=4, + self.vdata.lt.item(), mag_out[coord][1], places=4, msg="unexpected MLT") self.assertRegex(self.vdata.loc_coord, "magnetic") @@ -552,9 +552,9 @@ def test_update_loc_coords_float(self): coord=coord) # Evaluate the output; note the loss of precision - self.assertAlmostEqual(float(self.vdata.lat), 75.0, places=1, + self.assertAlmostEqual(self.vdata.lat.item(), 75.0, places=1, msg="unexpected geographic latitude") - self.assertAlmostEqual(float(self.vdata.lt), 22.0, places=1, + self.assertAlmostEqual(self.vdata.lt.item(), 22.0, places=1, msg="unexpected SLT") self.assertRegex(self.vdata.loc_coord, coord) return @@ -722,15 +722,15 @@ def test_update_vect_and_loc_coords_float(self): # Evaluate the output self.assertAlmostEqual( - float(self.vdata.vect_n), + self.vdata.vect_n.item(), mag_out[coord][loc_coord]['vect_n'], places=4, msg="unexpected north component") self.assertAlmostEqual( - float(self.vdata.vect_e), + self.vdata.vect_e.item(), mag_out[coord][loc_coord]['vect_e'], places=4, msg="unexpected east component") self.assertAlmostEqual( - float(self.vdata.vect_z), + self.vdata.vect_z.item(), mag_out[coord][loc_coord]['vect_z'], places=4, msg="unexpected vertical component") self.assertRegex(self.vdata.vect_coord, "magnetic") @@ -738,10 +738,10 @@ def test_update_vect_and_loc_coords_float(self): if loc_coord in mag_out.keys(): self.assertAlmostEqual( - float(self.vdata.lat), mag_out[loc_coord]['lat'], + self.vdata.lat.item(), mag_out[loc_coord]['lat'], places=4, msg="unexpected magnetic latitude") self.assertAlmostEqual( - float(self.vdata.lt), mag_out[loc_coord]['lt'], + self.vdata.lt.item(), mag_out[loc_coord]['lt'], places=4, msg="unexpected MLT") return diff --git a/ocbpy/vectors.py b/ocbpy/vectors.py index d176eb79..83427c4d 100644 --- a/ocbpy/vectors.py +++ b/ocbpy/vectors.py @@ -98,10 +98,10 @@ def calc_vec_pole_angle(data_lt, data_lat, pole_lt, pole_lat): else: flat_mask = (((del_long == 0) | (abs(del_long) == np.pi)) & np.greater(abs(data_lat), abs(pole_lat), - where=~np.isnan(del_long))) + where=~np.isnan(del_long), out=None)) zero_mask = (((del_long == 0) | (abs(del_long) == np.pi)) & np.less_equal(abs(data_lat), abs(pole_lat), - where=~np.isnan(del_long))) + where=~np.isnan(del_long), out=None)) pole_angle[flat_mask] = 180.0 pole_angle[zero_mask] = 0.0 @@ -176,39 +176,43 @@ def define_pole_quadrants(data_lt, pole_lt, pole_angle): pole_quad = np.zeros(shape=np.asarray(del_lt).shape) # Determine which differences need to be - neg_mask = np.less(del_lt, 0.0, where=~np.isnan(del_lt)) & ~np.isnan(del_lt) + neg_mask = np.less(del_lt, 0.0, where=~np.isnan(del_lt), + out=None) & ~np.isnan(del_lt) while np.any(neg_mask): if len(del_lt.shape) == 0: del_lt += 24.0 - neg_mask = np.less(del_lt, 0.0) # Has one finite value + neg_mask = np.less(del_lt, 0.0, out=None) # Has one finite value else: del_lt[neg_mask] += 24.0 - neg_mask = np.less(del_lt, 0.0, + neg_mask = np.less(del_lt, 0.0, out=None, where=~np.isnan(del_lt)) & ~np.isnan(del_lt) - large_mask = np.greater_equal(abs(del_lt), 24.0, + large_mask = np.greater_equal(abs(del_lt), 24.0, out=None, where=~np.isnan(del_lt)) & ~np.isnan(del_lt) while np.any(large_mask): if len(del_lt.shape) == 0: del_lt -= 24.0 * np.sign(del_lt) - large_mask = np.greater_equal(abs(del_lt), 24.0) # One finite value + # One finite value + large_mask = np.greater_equal(abs(del_lt), 24.0, out=None) else: del_lt[large_mask] -= 24.0 * np.sign(del_lt[large_mask]) - large_mask = np.greater_equal(abs(del_lt), 24.0, + large_mask = np.greater_equal(abs(del_lt), 24.0, out=None, where=~np.isnan(del_lt)) & ~np.isnan( del_lt) # Find the quadrant in which the OCB pole lies nan_mask = ~np.isnan(pole_angle) & ~np.isnan(del_lt) - quad1_mask = np.less(pole_angle, 90.0, where=nan_mask) & np.less( - del_lt, 12.0, where=nan_mask) & nan_mask - quad2_mask = np.less(pole_angle, 90.0, where=nan_mask) & np.greater_equal( - del_lt, 12.0, where=nan_mask) & nan_mask + quad1_mask = np.less(pole_angle, 90.0, where=nan_mask, out=None) & np.less( + del_lt, 12.0, where=nan_mask, out=None) & nan_mask + quad2_mask = np.less( + pole_angle, 90.0, where=nan_mask, out=None) & np.greater_equal( + del_lt, 12.0, where=nan_mask, out=None) & nan_mask quad3_mask = np.greater_equal( - pole_angle, 90.0, where=nan_mask) & np.greater_equal( - del_lt, 12.0, where=nan_mask) & nan_mask - quad4_mask = np.greater_equal(pole_angle, 90.0, where=nan_mask) & np.less( - del_lt, 12.0, where=nan_mask) & nan_mask + pole_angle, 90.0, where=nan_mask, out=None) & np.greater_equal( + del_lt, 12.0, where=nan_mask, out=None) & nan_mask + quad4_mask = np.greater_equal( + pole_angle, 90.0, where=nan_mask, out=None) & np.less( + del_lt, 12.0, where=nan_mask, out=None) & nan_mask if len(pole_quad.shape) == 0: if np.all(quad1_mask): @@ -255,14 +259,16 @@ def define_vect_quadrants(vect_n, vect_e): # Get the masks for non-fill values in each quadrant nan_mask = ~np.isnan(vect_n) & ~np.isnan(vect_e) quad1_mask = np.greater_equal( - vect_n, 0.0, where=nan_mask) & np.greater_equal( - vect_e, 0.0, where=nan_mask) & nan_mask - quad2_mask = np.greater_equal(vect_n, 0.0, where=nan_mask) & np.less( - vect_e, 0.0, where=nan_mask) & nan_mask - quad3_mask = np.less(vect_n, 0.0, where=nan_mask) & np.less( - vect_e, 0.0, where=nan_mask) & nan_mask - quad4_mask = np.less(vect_n, 0.0, where=nan_mask) & np.greater_equal( - vect_e, 0.0, where=nan_mask) & nan_mask + vect_n, 0.0, where=nan_mask, out=None) & np.greater_equal( + vect_e, 0.0, where=nan_mask, out=None) & nan_mask + quad2_mask = np.greater_equal( + vect_n, 0.0, where=nan_mask, out=None) & np.less( + vect_e, 0.0, where=nan_mask, out=None) & nan_mask + quad3_mask = np.less(vect_n, 0.0, where=nan_mask, out=None) & np.less( + vect_e, 0.0, where=nan_mask, out=None) & nan_mask + quad4_mask = np.less( + vect_n, 0.0, where=nan_mask, out=None) & np.greater_equal( + vect_e, 0.0, where=nan_mask, out=None) & nan_mask # Initialize the output vect_quad = np.zeros(shape=nan_mask.shape) @@ -466,13 +472,14 @@ def calc_dest_vec_sign(pole_quad, vect_quad, base_naz_angle, pole_angle, pmask = (quads[1][1] | quads[2][2] | quads[3][3] | quads[4][4] | ((quads[1][4] | quads[2][3]) & np.less_equal( - base_naz_angle, pole_plus, where=nan_mask)) + base_naz_angle, pole_plus, where=nan_mask, out=None)) | ((quads[1][2] | quads[2][1]) & np.less_equal( - base_naz_angle, minus_pole, where=nan_mask)) + base_naz_angle, minus_pole, where=nan_mask, out=None)) | ((quads[3][4] | quads[4][3]) & np.greater_equal( - base_naz_angle, 180.0 - pole_minus, where=nan_mask)) + base_naz_angle, 180.0 - pole_minus, where=nan_mask, + out=None)) | ((quads[3][2] | quads[4][1]) & np.greater_equal( - base_naz_angle, pole_minus, where=nan_mask))) + base_naz_angle, pole_minus, where=nan_mask, out=None))) if np.any(pmask): if len(vsigns["north"].shape) == 0: @@ -491,13 +498,13 @@ def calc_dest_vec_sign(pole_quad, vect_quad, base_naz_angle, pole_angle, pmask = (quads[1][4] | quads[2][1] | quads[3][2] | quads[4][3] | ((quads[1][1] | quads[4][4]) & np.greater_equal( - base_naz_angle, pole_angle, where=nan_mask)) + base_naz_angle, pole_angle, where=nan_mask, out=None)) | ((quads[3][1] | quads[2][4]) & np.less_equal( - base_naz_angle, minus_pole, where=nan_mask)) + base_naz_angle, minus_pole, where=nan_mask, out=None)) | ((quads[4][2] | quads[1][3]) & np.greater_equal( - base_naz_angle, minus_pole, where=nan_mask)) + base_naz_angle, minus_pole, where=nan_mask, out=None)) | ((quads[2][2] | quads[3][3]) & np.less_equal( - base_naz_angle, pole_angle, where=nan_mask))) + base_naz_angle, pole_angle, where=nan_mask, out=None))) if np.any(pmask): if len(vsigns["east"].shape) == 0: @@ -618,7 +625,8 @@ def adjust_vector(vect_lt, vect_lat, vect_n, vect_e, vect_z, vect_quad, # Determine if the measurement is on or between the poles. This does # not affect the vertical direction sign_mask = (pole_angle == 0.0) & np.greater_equal( - vect_lat, pole_lat, where=~np.isnan(vect_lat)) & ~np.isnan(vect_lat) + vect_lat, pole_lat, where=~np.isnan(vect_lat), + out=None) & ~np.isnan(vect_lat) if np.any(sign_mask): if len(out_shape) == 0: diff --git a/pyproject.toml b/pyproject.toml index e843883a..f654358f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "ocbpy" -version = "0.6.0" +version = "0.7.0" license = {file = "LICENSE"} description = 'Convert between magnetic/geodetic and adaptive, polar boundary coordinates' maintainers = [ @@ -12,7 +12,7 @@ maintainers = [ ] requires-python = ">=3.10" dependencies = [ - "aacgmv2", + "aacgmv2>=2.7.1", "numpy", ] readme = "README.md" @@ -46,6 +46,7 @@ classifiers = [ "Programming Language :: Python :: 3.11", "Programming Language :: Python :: 3.12", "Programming Language :: Python :: 3.13", + "Programming Language :: Python :: 3.14", 'Operating System :: Unix', 'Operating System :: POSIX', 'Operating System :: POSIX :: Linux', @@ -55,14 +56,14 @@ classifiers = [ [project.optional-dependencies] pysat_instruments = [ "pysat>=3.2.1" ] -dmsp_ssj = [ "zenodo-get" ] +dmsp_ssj = [ "zenodo-get>=2.0.0" ] doc = [ "numpydoc", "pyproject-parser", "pysat>=3.2.1", "sphinx>=1.3", "sphinx-rtd-theme", - "zenodo-get", + "zenodo-get>=2.0.0", ] test = [ "coverage[toml]",