Added a boundary models

Changelog.rst

@@ -7,6 +7,8 @@ Summary of all changes made since the first stable release
 ------------------
 * ENH: Updated the AMPERE boundaries to include 2022-2024, inclusive
 * ENH: Added `to_dict` method to the boundary class objects
+* ENH: Added Starkov (1994) auroral model
+* ENH: Adapted boundary classes to accept model boundaries
 
 0.5.0 (01-28-2025)
 ------------------

README.md

@@ -117,6 +117,7 @@ YYYY-MM-DD HH:MM:SS Phi_Centre R_Centre R
 -----------------------------------------------------------------------------
 2000-05-04 03:03:20 4.64 2.70 21.00
 2000-05-04 03:07:15 147.24 2.63 7.09
+...
 2002-10-31 20:03:16 207.11 5.94 22.86
 2002-10-31 20:05:16 335.47 6.76 11.97
 

docs/citing.rst

@@ -99,3 +99,16 @@ boundary method and data set are provided below.
   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
+
+
+.. _cite-starkov:
+
+Starkov Model
+-------------
+
+The Starkov mathematical auroral boundary model is described in the following
+article. If you use this model please cite both it and your AL data source
+when publishing your work.
+
+* Starkov, G. V. (1994) Mathematical model of the auroral boundaries,
+  Geomagnetism and Aeronomy, English Translation, 34(3), 331-336.

docs/example.rst

@@ -15,3 +15,4 @@ convert between AACGM and OCB coordinates, and more.
    examples/ex_vector.rst
    examples/ex_pysat_eab.rst
    examples/ex_save_boundaries.rst
+   examples/ex_model_boundaries.rst

docs/examples/ex_dmsp.rst

@@ -71,6 +71,7 @@ instrument, and hemisphere or merely the instrument and hemisphere.
    -----------------------------------------------------------------------------
    2010-12-31 00:27:23 356.72 14.02 12.13
    2010-12-31 12:27:56 324.82 0.86 14.73
+   ...
    2010-12-31 18:49:58 233.68 6.12 14.10
    2010-12-31 22:11:38 318.60 4.64 12.34
 
@@ -92,6 +93,7 @@ instrument, and hemisphere or merely the instrument and hemisphere.
    -----------------------------------------------------------------------------
    2010-12-31 01:19:13 191.07 10.69 8.59
    2010-12-31 06:27:18 195.29 13.52 6.77
+   ...
    2010-12-31 21:21:32 259.27 2.73 10.45
    2010-12-31 23:02:48 234.73 3.94 10.79
 

docs/examples/ex_init.rst

@@ -8,6 +8,8 @@ for most functions and examples, is the Open-Closed field line Boundary class,
 :py:class:`~ocbpy._boundary.OCBoundary`.
 
 
+.. _exinit-ocb:
+
 Initialise an OCBoundary object
 -------------------------------
 Start a python or iPython session, and begin by importing ocbpy, numpy,
@@ -41,13 +43,16 @@ the default IMAGE FUV file and will take a few minutes to load.
    -----------------------------------------------------------------------------
    2000-05-04 03:03:20 4.64 2.70 21.00
    2000-05-04 03:07:15 147.24 2.63 7.09
+   ...
    2002-10-31 20:03:16 207.11 5.94 22.86
    2002-10-31 20:05:16 335.47 6.76 11.97
 
    Uses scaling function(s):
    ocbpy.ocb_correction.circular(**{})
 
 
+.. _exinit-other:
+
 Other Boundary classes
 ----------------------
 
@@ -83,6 +88,7 @@ classes.
    -----------------------------------------------------------------------------
    2000-05-05 11:35:27 111.80 2.34 25.12
    2000-05-05 11:37:23 296.23 1.42 26.57
+   ...
    2000-05-07 21:50:20 220.84 7.50 16.89
    2000-05-07 23:32:14 141.27 10.33 18.32
 
@@ -99,9 +105,63 @@ classes.
    -----------------------------------------------------------------------------
    2000-05-05 11:19:52 218.54 9.44 11.48
    2000-05-05 11:35:27 304.51 8.69 15.69
+   ...
    2000-05-07 23:24:14 199.84 10.91 12.69
    2000-05-07 23:32:14 141.53 9.24 13.03
 
 
    Uses scaling function(s):
    ocbpy.ocb_correction.circular(**{})
+
+
+.. _exinit-model:
+
+Initialising with a Model
+-------------------------
+
+Any of these boundary classes may be initialized with a modelled specification
+of the appropriate boundary or boundaries, instead of measurements. This
+requires two major changes. First, a boundary model function that specifies the
+location of the desired boundary in degrees latitude away from the pole must
+be provided through the :py:attr:`rfunc` keyword. This function should have
+magnetic local time (MLT) as the input argument; any subsequent inputs should be
+keyword arguments.  Second, the user must provide an array or list of input
+times for boundary calculations through the :py:attr:`stime` keyword argument.
+Model inputs other than MLT are then provided through the :py:attr:`rfunc_kwarg`
+keyword.
+
+
+::
+
+
+   # Set the hourly AL values for an active period
+   stime = [dt.datetime(2000, 1, 4, 4 + i) for i in range(5)]
+   al_list = [-36, -87, -104, -166, -326]
+
+   # Initalize the class for the desired period
+   starkov = ocbpy.OCBoundary(
+       instrument='model', stime=stime,
+       rfunc=ocbpy.boundaries.models.starkov_auroral_boundary,
+       rfunc_kwargs={'al': al_list, 'bnd': 'ocb'})
+   print(starkov)
+
+   OCBoundary
+   Source instrument: MODEL
+   Boundary reference latitude: 74.0 degrees
+
+   5 records from 2000-01-04 04:00:00 to 2000-01-04 08:00:00
+
+   YYYY-MM-DD HH:MM:SS Phi_Centre R_Centre R
+   -----------------------------------------------------------------------------
+   2000-01-04 04:00:00 0.00 0.00 0.00
+   2000-01-04 05:00:00 0.00 0.00 0.00
+   ...
+   2000-01-04 07:00:00 0.00 0.00 0.00
+   2000-01-04 08:00:00 0.00 0.00 0.00
+
+   Uses boundary function(s):
+   ocbpy.boundaries.models.starkov_auroral_boundary(**{'al': -36, 'bnd': 'ocb'})
+   ocbpy.boundaries.models.starkov_auroral_boundary(**{'al': -87, 'bnd': 'ocb'})
+   ...
+   ocbpy.boundaries.models.starkov_auroral_boundary(**{'al': -166, 'bnd': 'ocb'})
+   ocbpy.boundaries.models.starkov_auroral_boundary(**{'al': -326, 'bnd': 'ocb'})

docs/examples/ex_model_boundaries.rst

@@ -0,0 +1,102 @@
+.. _exmodel:
+
+Using Model Boundaries
+======================
+
+The boundary models supplied in this package can be used in various ways. The
+boundary locations themselves can be calculated for a range of MLTs and
+upstream conditions and extracted for various uses. They can also be used to
+grid data in adaptive coordinates, as if they were a boundary data set. This may
+be useful for statistical studies where sufficient observations of the desired
+boundaries are not available.
+
+Calculate Boundary Locations
+----------------------------
+
+Following on from the model OCBoundary initialisation example
+(see :ref:`exinit-model`), the OCB locations for the desired period of time can
+be calculated for a set of MLTs using the
+:py:meth:`~ocbpy.OCBoundary.get_aacgm_boundary_lat` method.
+
+
+::
+
+   import numpy as np
+
+   # Define the desired MLT range
+   mlt = np.arange(0, 24, 0.5)
+
+   # Calculate the boundaries at all times
+   starkov.get_aacgm_boundary_lat(mlt)
+
+
+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`.
+
+::
+
+
+   # Set up fthe figure
+   fig = plt.figure()
+   ax = fig.add_subplot(111, projection='polar')
+   set_up_polar_plot(ax)
+
+   # Construct an array of AL values to color code the boundaries
+   al_array = np.full(shape=(mlt.shape[0], starkov.records),
+                      fill_value=al_list).transpose()
+
+   # Plot the boundaries
+   con = ax.scatter(np.asarray(starkov.aacgm_boundary_mlt) * np.pi / 12.0,
+                    90 - np.asarray(starkov.aacgm_boundary_lat), c=al_array,
+		    marker='.', vmin=-400, vmax=0,
+		    cmap=mpl.colormaps.get_cmap('plasma'))
+
+   # Add a colorbar
+   cb = fig.colorbar(con, ax=ax, orientation='horizontal')
+   cb.set_label('AL (nT)')
+
+
+.. image:: ../figures/example_model_starkov_ocb.png
+
+
+Calculate Boundary Locations
+----------------------------
+
+You can also use these empirical boundaries to grid data. When using models to
+do this, be aware that the accuracy of the model will impact the results of any
+statistical studies. The example below tracks the location of a location near
+the edge of the polar cap at magnetic midnight.
+
+::
+
+
+   # Ensure the record index iterator is at zero to start
+   starkov.rec_ind = 0
+
+   # Cycle through each time and find the OCB location of the magnetic pole
+   ocb_lat = list()
+   ocb_mlt = list()
+   while starkov.rec_ind < starkov.records:
+       out = starkov.normal_coord(75.0, 0.0, height=0.0)
+       ocb_lat.append(out[0])
+       ocb_mlt.append(float(out[1]))
+       starkov.rec_ind += 1
+
+   # Note that this model is consistently centred at the same location, so
+   # MLT does not change
+   print(ocb_mlt)  # Will show all 0.0
+
+   # Plot the change in OCB latitude as a function of AL
+   fig = plt.figure()
+   ax = fig.add_subplot(111)
+   ax.plot([rkwargs['al'] for rkwargs in starkov.rfunc_kwargs], ocb_lat, 'ko')
+
+   # Format the figure
+   ax.set_xlabel('AL (nT)')
+   ax.set_ylabel('OCB Latitude of 75$^\circ$ ($\circ$)')
+   fig.suptitle("Starkov (1994) Model Gridded Location")
+   fig.tight_layout()
+
+
+.. image:: ../figures/example_model_starkov_ocb_lat_vs_al.png

docs/examples/ex_pysat_eab.rst

@@ -34,6 +34,7 @@ very similar to :ref:`exinit`, the first example in this section.
     ------------------------------------------
     2000-05-03 02:41:43 251.57 5.45 14.16
     2000-05-05 11:35:27 111.80 2.34 25.12
+    ...
     2002-10-31 20:03:16 354.41 6.22 20.87
     2002-10-31 20:05:16 2.87 14.67 13.10
 

docs/examples/ex_vector.rst

@@ -69,6 +69,7 @@ data set available.
     ----------------------------------------------------------------------------
     2001-02-14 01:33:24 169.01 1.94 15.85
     2001-02-14 01:35:26 170.18 1.56 16.47
+    ...
     2001-02-14 01:53:51 239.81 1.57 15.70
     2001-02-14 01:55:54 210.02 2.09 15.89
 
@@ -239,6 +240,7 @@ Now let us transform the same data using both the EAB and OCB.
    -----------------------------------------------------------------------------
    2001-02-14 01:33:24 26.27 3.24 26.76
    2001-02-14 01:35:26 26.54 3.59 26.53
+   ...
    2001-02-14 01:53:51 9.72 4.20 26.37
    2001-02-14 01:55:54 13.46 3.06 26.78
 
@@ -255,6 +257,7 @@ Now let us transform the same data using both the EAB and OCB.
    -----------------------------------------------------------------------------
    2001-02-14 01:33:24 169.01 1.94 15.85
    2001-02-14 01:35:26 170.18 1.56 16.47
+   ...
    2001-02-14 01:53:51 239.81 1.57 15.70
    2001-02-14 01:55:54 210.02 2.09 15.89
 

docs/index.rst

@@ -15,6 +15,7 @@ process and provides examples for usage.
    overview.rst
    citing.rst
    ocb_datasets.rst
+   ocb_models.rst
    ocb_gridding.rst
    ocb_boundary_correction.rst
    supported_datasets.rst

docs/ocb_models.rst

@@ -0,0 +1,35 @@
+.. _bound-model:
+
+
+Boundary Models
+===============
+
+Another option for providing boundaries are model specifications. Several
+empirically derived mathematical models are included in
+:py:mod:`ocbpy.boundaries.models` to allow access to these formulations. These
+models typically depend on magnetic local time (MLT) and a geomagnetic or solar
+wind index. OCBpy requires that these functions provide the boundary location in
+co-latitude (degrees of magnetic latitude away from the pole), have MLT as
+the first input argument, and that all other inputs be keyword arguments.
+
+
+.. _bound-model-starkov:
+
+Starkov
+-------
+
+The Starkov 1994 model (see :ref:`cite-starkov`) uses a mathematical formulation
+based on All-Sky Imager data and the Auroral Electrojet Lower envelope index.
+They specify three boundaries: the polar edge of the auroral oval, the
+equatorward edge of the auroral oval, and the equatorward edge of the diffuse
+aurora (usually more equatorward than the discrete edge). These may be accessed
+in the code here using the boundary keyworkds: 'ocb', 'eab', and 'diffuse',
+respectively.
+
+
+.. _bound-model-module:
+
+Boundary Models Module
+----------------------
+.. automodule:: ocbpy.boundaries.models
+    :members: