diff --git a/Changelog.rst b/Changelog.rst index 13368e1c..c2e80bfd 100644 --- a/Changelog.rst +++ b/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) ------------------ diff --git a/README.md b/README.md index 25d023eb..1abf5c72 100644 --- a/README.md +++ b/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 diff --git a/docs/citing.rst b/docs/citing.rst index 9b0ff895..fee63c5d 100644 --- a/docs/citing.rst +++ b/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. diff --git a/docs/example.rst b/docs/example.rst index 26b19ef3..2bdd7897 100644 --- a/docs/example.rst +++ b/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 diff --git a/docs/examples/ex_dmsp.rst b/docs/examples/ex_dmsp.rst index 441b3a0b..069037ad 100644 --- a/docs/examples/ex_dmsp.rst +++ b/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 diff --git a/docs/examples/ex_init.rst b/docs/examples/ex_init.rst index 2ad530f5..e1da5665 100644 --- a/docs/examples/ex_init.rst +++ b/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,6 +43,7 @@ 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 @@ -48,6 +51,8 @@ the default IMAGE FUV file and will take a few minutes to load. 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'}) diff --git a/docs/examples/ex_model_boundaries.rst b/docs/examples/ex_model_boundaries.rst new file mode 100644 index 00000000..a0731001 --- /dev/null +++ b/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 diff --git a/docs/examples/ex_pysat_eab.rst b/docs/examples/ex_pysat_eab.rst index e2e40f9d..19d95f27 100644 --- a/docs/examples/ex_pysat_eab.rst +++ b/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 diff --git a/docs/examples/ex_vector.rst b/docs/examples/ex_vector.rst index 4323955e..60257033 100644 --- a/docs/examples/ex_vector.rst +++ b/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 diff --git a/docs/figures/example_model_starkov_ocb.png b/docs/figures/example_model_starkov_ocb.png new file mode 100644 index 00000000..e699d3ba Binary files /dev/null and b/docs/figures/example_model_starkov_ocb.png differ diff --git a/docs/figures/example_model_starkov_ocb_lat_vs_al.png b/docs/figures/example_model_starkov_ocb_lat_vs_al.png new file mode 100644 index 00000000..6af3a018 Binary files /dev/null and b/docs/figures/example_model_starkov_ocb_lat_vs_al.png differ diff --git a/docs/index.rst b/docs/index.rst index 2108b3e2..6df64006 100644 --- a/docs/index.rst +++ b/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 diff --git a/docs/ocb_models.rst b/docs/ocb_models.rst new file mode 100644 index 00000000..fc4b3652 --- /dev/null +++ b/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: diff --git a/ocbpy/_boundary.py b/ocbpy/_boundary.py index 931e66c8..ad0f6210 100644 --- a/ocbpy/_boundary.py +++ b/ocbpy/_boundary.py @@ -43,33 +43,39 @@ class OCBoundary(object): If NoneType, no file is loaded. If 'default', `ocbpy.boundaries.files.get_default_file` is called. (default='default') instrument : str - Instrument providing the OCBoundaries. Requires 'image', 'ampere', or - 'dmsp-ssj' if a file is provided. If using filename='default', also - accepts 'amp', 'si12', 'si13', 'wic', and ''. (default='') + Instrument providing the OCBoundaries. Requires 'image', 'ampere', or + 'dmsp-ssj' if a file is provided through `filename`. If using + filename='default', also accepts 'amp', 'si12', 'si13', 'wic', 'model', + and ''. (default='') hemisphere : int Integer (+/- 1) denoting northern/southern hemisphere (default=1) boundary_lat : float Typical OCBoundary latitude in AACGM coordinates. Hemisphere will give this boundary the desired sign. (default=74.0) - stime : dt.datetime or NoneType + stime : dt.datetime, array-like, or NoneType First time to load data or beginning of file. If specifying time, be sure to start before the time of the data to allow the best match within - the allowable time tolerance to be found. (default=None) + the allowable time tolerance to be found. If running a model (e.g., + `instrument='model'`, provide an array of datetime values + corresponding to the required driving index(es). (default=None) etime : dt.datetime or NoneType Last time to load data or ending of file. If specifying time, be sure to end after the last data point you wish to match to, to ensure the best match within the allowable time tolerance is made. (default=None) rfunc : numpy.ndarray, function, or NoneType - OCB radius correction function. If None, will use the instrument - default. Function must have AACGM MLT (in hours) as argument input. - To allow the boundary shape to change with univeral time, each temporal - instance may have a different function (array input). If a single - function is provided, will recast as an array that specifies this - function for all times. (default=None) + OCB radius correction or model function (if `instrument` is 'model'). + If None, will use the instrument default. Function must have AACGM MLT + (in hours) as argument input. To allow the boundary shape to change + with univeral time, each temporal instance may have a different function + (array input). If a single function is provided, will recast as an array + that specifies this function for all times. (default=None) rfunc_kwargs : numpy.ndarray, dict, or NoneType - Optional keyword arguements for `rfunc`. If None is specified, - uses function defaults. If dict is specified, recasts as an array - of this dict for all times. Array must be an array of dicts. + Optional keyword arguments for `rfunc`. If None is specified, + uses function defaults. If dict is specified and `instrument` is + 'model', any array-like values are expected to be the same length as + `dtime` and will be recast as an array of dicts with single values. + Otherwise, dict inputs are recasts as an array of this dict for all + times. Array must be an array of dicts the same length as `dtime`. (default=None) Attributes @@ -104,7 +110,7 @@ class OCBoundary(object): inst_defaults Get the instrument-specific boundary file loading information. load - Load the data from the specified boundary file. + Load the data from the specified boundary file or model. get_next_good_ocb_ind Cycle to the the next quality boundary record. normal_coord @@ -136,7 +142,7 @@ def __init__(self, filename="default", instrument='', hemisphere=1, self.instrument = instrument.lower() # If a filename wanted and not provided, get one - if filename is None: + if filename is None or self.instrument == "model": self.filename = None elif not hasattr(filename, "lower"): logger.warning("filename is not a string [{:}]".format( @@ -174,9 +180,6 @@ def __init__(self, filename="default", instrument='', hemisphere=1, self.min_fom = -np.inf self.max_fom = np.inf - # Get the instrument defaults - hlines, ocb_cols, datetime_fmt = self.inst_defaults() - # Set the boundary latitude, if supplied self.boundary_lat = 74.0 if boundary_lat is None else boundary_lat @@ -186,11 +189,17 @@ def __init__(self, filename="default", instrument='', hemisphere=1, # If possible, load the data. Any boundary correction is applied here. if self.filename is not None: + # Get the Instrument defaults + hlines, ocb_cols, datetime_fmt = self.inst_defaults() + + # Choose the approriate loading method if len(ocb_cols) > 0: self.load(hlines=hlines, ocb_cols=ocb_cols, datetime_fmt=datetime_fmt, stime=stime, etime=etime) else: self.load(stime=stime, etime=etime) + elif self.instrument == "model": + self.load(stime=stime) return @@ -242,10 +251,14 @@ def __str__(self): class_name = repr(self.__class__).split("'")[1].split(".")[-1] - if self.filename is None: + if self.filename is None and self.instrument != "model": out = "No {:s} file specified\n".format(class_name) else: - out = "{:s} file: {:s}\n".format(class_name, self.filename) + if self.filename is not None: + out = "{:s} file: {:s}\n".format(class_name, self.filename) + else: + out = "{:s}\n".format(class_name) + out = "{:s}Source instrument: ".format(out) out = "{:s}{:s}\n".format(out, self.instrument.upper()) out = "{:s}Boundary reference latitude: ".format(out) @@ -258,11 +271,14 @@ def __str__(self): self.dtime[0]) out = "{:s} to {:}\n\n".format(out, self.dtime[-1]) + imid = -1 if self.records == 1: irep = [0] else: irep = np.unique( np.arange(0, self.records, 1)[[0, 1, -2, -1]]) + if self.records > 4: + imid = 1 head = "YYYY-MM-DD HH:MM:SS Phi_Centre R_Centre R" out = "{:s}{:s}\n{:-<77s}\n".format(out, head, "") @@ -271,66 +287,111 @@ def __str__(self): self.phi_cent[i]) out = "{:s} {:.2f} {:.2f}\n".format(out, self.r_cent[i], self.r[i]) + if i == imid: + out = "".join([out, "...\n"]) # Determine which scaling functions are used if self.rfunc is not None: - out = "{:s}\nUses scaling function(s):\n".format(out) + out = "".join([out, "\nUses ", "boundary" + if self.instrument == "model" else "scaling", + " function(s):\n"]) fnames = list(set([".".join([ff.__module__, ff.__name__]) for ff in self.rfunc])) for ff in fnames: - kw = list(set([self.rfunc_kwargs[i].__str__() - for i, rf in enumerate(self.rfunc) - if rf.__name__ == ff.split(".")[-1]])) - - for kk in kw: - out = "{:s}{:s}(**{:s})\n".format(out, ff, kk) + kw = [self.rfunc_kwargs[i].__str__() + for i, rf in enumerate(self.rfunc) + if rf.__name__ == ff.split(".")[-1]] + + ikw = list(irep) + kmid = imid + if len(kw) < len(ikw): + ikw = [kind for kind in range(len(kw))] + if len(ikw) <= 1: + kmid = -1 + + for k in ikw: + out = "{:s}{:s}(**{:s})\n".format(out, ff, kw[k]) + if k == kmid: + out = "".join([out, "...\n"]) return out - def inst_defaults(self): - """Get the instrument-specific OCB file loading information. - - Returns - ------- - hlines : int - Number of header lines - ocb_cols : str - String containing the names for each data column - datetime_fmt : str - String containing the datetime format + def _set_default_rfunc(self): + """Set the default instrument OCB boundary function. Notes ----- - Updates the min_fom attribute for AMPERE and DMSP-SSJ + Assign a function for each time in case we have a data set with a + correction that changes with UT """ - if self.instrument == "image": - hlines = 0 - ocb_cols = "year soy num_sectors phi_cent r_cent r a r_err fom" - datetime_fmt = "" - self.max_fom = 5.0 # From Chisham et al. (2022) + if self.instrument in ["image", "dmsp-ssj"]: + self.rfunc = np.full(shape=self.records, + fill_value=ocbcor.circular) elif self.instrument == "ampere": - hlines = 0 - ocb_cols = "date time r x y fom" - datetime_fmt = "%Y%m%d %H:%M" - self.min_fom = 0.15 # From Milan et al. (2015) - elif self.instrument == "dmsp-ssj": - hlines = 1 - ocb_cols = "sc date time r x y fom x_1 y_1 x_2 y_2" - datetime_fmt = "%Y-%m-%d %H:%M:%S" - self.min_fom = 3.0 # From Burrell et al. (2019) + self.rfunc = np.full(shape=self.records, + fill_value=ocbcor.elliptical) else: - hlines = 0 - ocb_cols = "" - datetime_fmt = "" + raise ValueError("unknown instrument: {:}".format(self.instrument)) - return hlines, ocb_cols, datetime_fmt + return - def load(self, hlines=0, - ocb_cols="year soy num_sectors phi_cent r_cent r a r_err fom", - datetime_fmt="", stime=None, etime=None): + def _set_rfunc(self): + """Adjust the formatting of the boundary function and kwargs.""" + # Set the boundary function + if self.rfunc is None: + self._set_default_rfunc() + elif isinstance(self.rfunc, types.FunctionType): + self.rfunc = np.full(shape=self.records, fill_value=self.rfunc) + elif hasattr(self.rfunc, "shape"): + if self.rfunc.shape != self.dtime.shape: + raise ValueError("Misshaped correction function array") + else: + raise ValueError("Unknown input type for correction function") + + # Set the boundary function keyword inputs + if self.rfunc_kwargs is None: + self.rfunc_kwargs = np.full(shape=self.records, fill_value={}) + elif isinstance(self.rfunc_kwargs, dict): + if self.instrument == "model": + # Inputs may be an array that needs to be parsed + reshape = False + for key in self.rfunc_kwargs.keys(): + val = self.rfunc_kwargs[key] + if not hasattr(val, 'lower') and len(val) == self.records: + if reshape is False: + reshape = [key] + else: + reshape.append(key) + else: + reshape = False + + if reshape is False: + self.rfunc_kwargs = np.full(shape=self.records, + fill_value=self.rfunc_kwargs) + else: + kwarg_list = list() + for i in range(self.records): + kwarg_list.append({ + key: self.rfunc_kwargs[key][i] + if key in reshape else self.rfunc_kwargs[key] + for key in self.rfunc_kwargs.keys()}) + + self.rfunc_kwargs = np.asarray(kwarg_list) + elif hasattr(self.rfunc_kwargs, "shape"): + if self.rfunc_kwargs.shape != self.dtime.shape: + raise ValueError("Misshaped correction function keyword array") + else: + raise ValueError("Unknown input type for correction keywords") + + return + + def _load_file( + self, hlines=0, + ocb_cols="year soy num_sectors phi_cent r_cent r a r_err fom", + datetime_fmt="", stime=None, etime=None): """Load the data from the specified boundary file. Parameters @@ -354,7 +415,7 @@ def load(self, hlines=0, (default=None) """ - + # Get the column formatting cols = ocb_cols.split() dflag = -1 ldtype = [(k, float) if k != "num_sectors" else (k, int) for k in cols] @@ -424,28 +485,8 @@ def load(self, hlines=0, self.records = len(dt_list) self.dtime = np.array(dt_list) - # Set the boundary function - if self.rfunc is None: - self._set_default_rfunc() - elif isinstance(self.rfunc, types.FunctionType): - self.rfunc = np.full(shape=self.records, fill_value=self.rfunc) - elif hasattr(self.rfunc, "shape"): - if self.rfunc.shape != self.dtime.shape: - raise ValueError("Misshaped correction function array") - else: - raise ValueError("Unknown input type for correction function") - - # Set the boundary function keyword inputs - if self.rfunc_kwargs is None: - self.rfunc_kwargs = np.full(shape=self.records, fill_value={}) - elif isinstance(self.rfunc_kwargs, dict): - self.rfunc_kwargs = np.full(shape=self.records, - fill_value=self.rfunc_kwargs) - elif hasattr(self.rfunc_kwargs, "shape"): - if self.rfunc_kwargs.shape != self.dtime.shape: - raise ValueError("Misshaped correction function keyword array") - else: - raise ValueError("Unknown input type for correction keywords") + # Set the boundary function and function keyword arguments + self._set_rfunc() # Load the attributes saved in odata for nn in oname: @@ -453,6 +494,105 @@ def load(self, hlines=0, return + def _load_model(self, rec_times): + """Load the data from the specified boundary file. + + Parameters + ---------- + rec_times : array-like + Times for which the model will be run. + + """ + # Set the times and record numbers + self.dtime = np.asarray(rec_times) + self.records = len(rec_times) + + # Set the boundary function and function keyword arguments + self._set_rfunc() + + # Update the "observed" location to be a point at the pole + self.phi_cent = np.zeros(shape=self.records) + self.r_cent = np.zeros(shape=self.records) + self.r = np.zeros(shape=self.records) + self.fom = np.zeros(shape=self.records) + + return + + def inst_defaults(self): + """Get the instrument-specific OCB file loading information. + + Returns + ------- + hlines : int + Number of header lines + ocb_cols : str + String containing the names for each data column + datetime_fmt : str + String containing the datetime format + + Notes + ----- + Updates the min_fom attribute for AMPERE and DMSP-SSJ + + """ + + if self.instrument == "image": + hlines = 0 + ocb_cols = "year soy num_sectors phi_cent r_cent r a r_err fom" + datetime_fmt = "" + self.max_fom = 5.0 # From Chisham et al. (2022) + elif self.instrument == "ampere": + hlines = 0 + ocb_cols = "date time r x y fom" + datetime_fmt = "%Y%m%d %H:%M" + self.min_fom = 0.15 # From Milan et al. (2015) + elif self.instrument == "dmsp-ssj": + hlines = 1 + ocb_cols = "sc date time r x y fom x_1 y_1 x_2 y_2" + datetime_fmt = "%Y-%m-%d %H:%M:%S" + self.min_fom = 3.0 # From Burrell et al. (2019) + else: + hlines = 0 + ocb_cols = "" + datetime_fmt = "" + + return hlines, ocb_cols, datetime_fmt + + def load(self, hlines=0, + ocb_cols="year soy num_sectors phi_cent r_cent r a r_err fom", + datetime_fmt="", stime=None, etime=None): + """Load the data from the specified boundary file. + + Parameters + ---------- + ocb_cols : str + String specifying format of OCB file. All but the first two + columns must be included in the string, additional data values will + be ignored. If 'year soy' aren't used, expects + 'date time' in 'YYYY-MM-DD HH:MM:SS' format. + (default='year soy num_sectors phi_cent r_cent r a r_err r_merit') + hlines : int + Number of header lines preceeding data in the OCB file (default=0) + datetime_fmt : str + A string used to read in 'date time' data. Not used if 'year soy' + is specified. (default='') + stime : dt.datetime, array-like, or NoneType + Time to start loading data, array of times for model calculation, + or None to start at beginning of an instrument file. (default=None) + etime : datetime or NoneType + Time to stop loading data or None to end at the end of the file. + (default=None) + + """ + + if self.instrument == "model": + self._load_model(stime) + else: + self._load_file(hlines=hlines, ocb_cols=ocb_cols, + datetime_fmt=datetime_fmt, stime=stime, etime=etime) + + return + def get_next_good_ocb_ind(self, min_merit=None, max_merit=None, **kwargs): """Cycle to the the next quality OCB record. @@ -851,27 +991,6 @@ def get_aacgm_boundary_lat(self, aacgm_mlt, rec_ind=None, return - def _set_default_rfunc(self): - """Set the default instrument OCB boundary function. - - Notes - ----- - Assign a function for each time in case we have a data set with a - correction that changes with UT - - """ - - if self.instrument in ["image", "dmsp-ssj"]: - self.rfunc = np.full(shape=self.records, - fill_value=ocbcor.circular) - elif self.instrument == "ampere": - self.rfunc = np.full(shape=self.records, - fill_value=ocbcor.elliptical) - else: - raise ValueError("unknown instrument: {:}".format(self.instrument)) - - return - def to_dict(self, xarray_style=False, sel_inds=None): """Convert the boundary data into a pair of dictionaries. @@ -971,31 +1090,36 @@ class EABoundary(OCBoundary): instrument : str Instrument providing the EABoundaries. Requires 'image' or 'dmsp-ssj' if a file is provided. If using filename='default', also accepts - 'si12', 'si13', 'wic', and ''. (default='') + 'si12', 'si13', 'wic', 'model', and ''. (default='') hemisphere : int Integer (+/- 1) denoting northern/southern hemisphere (default=1) boundary_lat : float Typical EABoundary latitude in AACGM coordinates. Hemisphere will give this boundary the desired sign. (default=64.0) - stime : dt.datetime or NoneType + stime : dt.datetime, array-like, or NoneType First time to load data or beginning of file. If specifying time, be sure to start before the time of the data to allow the best match within - the allowable time tolerance to be found. (default=None) + the allowable time tolerance to be found. If running a model (e.g., + `instrument='model'`, provide an array of datetime values + corresponding to the required driving index(es). (default=None) etime : dt.datetime or NoneType Last time to load data or ending of file. If specifying time, be sure to end after the last data point you wish to match to, to ensure the best match within the allowable time tolerance is made. (default=None) rfunc : numpy.ndarray, function, or NoneType - EAB radius correction function. If None, will use the instrument - default. Function must have AACGM MLT (in hours) as argument input. - To allow the boundary shape to change with univeral time, each temporal - instance may have a different function (array input). If a single - function is provided, will recast as an array that specifies this - function for all times. (default=None) + EAB radius correction or model function (if `instrument` is 'model'). + If None, will use the instrument default. Function must have AACGM MLT + (in hours) as argument input. To allow the boundary shape to change + with univeral time, each temporal instance may have a different + function (array input). If a single function is provided, will recast + as an array that specifies this function for all times. (default=None) rfunc_kwargs : numpy.ndarray, dict, or NoneType - Optional keyword arguements for `rfunc`. If None is specified, - uses function defaults. If dict is specified, recasts as an array - of this dict for all times. Array must be an array of dicts. + Optional keyword arguments for `rfunc`. If None is specified, + uses function defaults. If dict is specified and `instrument` is + 'model', any array-like values are expected to be the same length as + `dtime` and will be recast as an array of dicts with single values. + Otherwise, dict inputs are recasts as an array of this dict for all + times. Array must be an array of dicts the same length as `dtime`. (default=None) See Also @@ -1021,7 +1145,7 @@ def __init__(self, filename="default", instrument='', hemisphere=1, self.rfunc = None if(hasattr(filename, "lower") and hasattr(instrument, "lower") - and filename.lower() == "default"): + and filename.lower() == "default" and instrument != "model"): filename, instrument = get_default_file(stime, etime, hemisphere, instrument, bound='eab') @@ -1076,11 +1200,11 @@ class DualBoundary(object): eab_instrument : str Instrument providing the EABoundaries. Requires 'image', 'ampere', or 'dmsp-ssj' if a file is provided. If using filename='default', also - accepts 'si12', 'si13', 'wic', and ''. (default='') + accepts 'si12', 'si13', 'wic', 'model', and ''. (default='') ocb_instrument : str Instrument providing the OCBoundaries. Requires 'image', 'ampere', or 'dmsp-ssj' if a file is provided. If using filename='default', also - accepts 'si12', 'si13', 'wic', and ''. (default='') + accepts 'si12', 'si13', 'wic', 'model', and ''. (default='') hemisphere : int Integer (+/- 1) denoting northern/southern hemisphere (default=1) eab_lat : float @@ -1089,37 +1213,45 @@ class DualBoundary(object): ocb_lat : float Typical OCBoundary latitude in AACGM coordinates. Hemisphere will give this boundary the desired sign. (default=74.0) - stime : dt.datetime or NoneType + stime : dt.datetime, array-like, or NoneType First time to load data or beginning of file. If specifying time, be sure to start before the time of the data to allow the best match within - the allowable time tolerance to be found. (default=None) + the allowable time tolerance to be found. If running a model (e.g., + `instrument='model'`, provide an array of datetime values + corresponding to the required driving index(es). (default=None) etime : dt.datetime or NoneType Last time to load data or ending of file. If specifying time, be sure to end after the last data point you wish to match to, to ensure the best match within the allowable time tolerance is made. (default=None) eab_rfunc : numpy.ndarray, function, or NoneType - EAB radius correction function. If None, will use the instrument - default. Function must have AACGM MLT (in hours) as argument input. - To allow the boundary shape to change with univeral time, each temporal - instance may have a different function (array input). If a single - function is provided, will recast as an array that specifies this - function for all times. (default=None) + EAB radius correction or model function (if `instrument` is 'model'). + If None, will use the instrument default. Function must have AACGM MLT + (in hours) as argument input. To allow the boundary shape to change + with univeral time, each temporal instance may have a different + function (array input). If a single function is provided, will recast + as an array that specifies this function for all times. (default=None) eab_rfunc_kwargs : numpy.ndarray, dict, or NoneType - Optional keyword arguements for `eab_rfunc`. If None is specified, - uses function defaults. If dict is specified, recasts as an array - of this dict for all times. Array must be an array of dicts. + Optional keyword arguments for `eab_rfunc`. If None is specified, + uses function defaults. If dict is specified and `instrument` is + 'model', any array-like values are expected to be the same length as + `dtime` and will be recast as an array of dicts with single values. + Otherwise, dict inputs are recasts as an array of this dict for all + times. Array must be an array of dicts the same length as `dtime`. (default=None) ocb_rfunc : numpy.ndarray, function, or NoneType - OCB radius correction function. If None, will use the instrument - default. Function must have AACGM MLT (in hours) as argument input. - To allow the boundary shape to change with univeral time, each temporal - instance may have a different function (array input). If a single - function is provided, will recast as an array that specifies this - function for all times. (default=None) + OCB radius correction or model function (if `instrument` is 'model'). + If None, will use the instrument default. Function must have AACGM MLT + (in hours) as argument input. To allow the boundary shape to change + with univeral time, each temporal instance may have a different function + (array input). If a single function is provided, will recast as an array + that specifies this function for all times. (default=None) ocb_rfunc_kwargs : numpy.ndarray, dict, or NoneType Optional keyword arguements for `ocb_rfunc`. If None is specified, - uses function defaults. If dict is specified, recasts as an array - of this dict for all times. Array must be an array of dicts. + uses function defaults. If dict is specified and `instrument` is + 'model', any array-like values are expected to be the same length as + `dtime` and will be recast as an array of dicts with single values. + Otherwise, dict inputs are recasts as an array of this dict for all + times. Array must be an array of dicts the same length as `dtime`. (default=None) eab : ocbpy.EABoundary or NoneType Equatorward auroral boundary data object or None to initialize here diff --git a/ocbpy/boundaries/__init__.py b/ocbpy/boundaries/__init__.py index 25bb7b05..8894afcf 100644 --- a/ocbpy/boundaries/__init__.py +++ b/ocbpy/boundaries/__init__.py @@ -10,6 +10,7 @@ from ocbpy import logger from ocbpy.boundaries import files # noqa F401 +from ocbpy.boundaries import models # noqa F401 try: from ocbpy.boundaries import dmsp_ssj_files # noqa F401 diff --git a/ocbpy/boundaries/models.py b/ocbpy/boundaries/models.py new file mode 100644 index 00000000..4c36d7fa --- /dev/null +++ b/ocbpy/boundaries/models.py @@ -0,0 +1,123 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# DOI: 10.5281/zenodo.1179230 +# Full license can be found in License.md +# ---------------------------------------------------------------------------- +"""Functions that provide boundary locations through a mathematical model. + +References +---------- +.. [8] Starkov, G. V. (1994) Mathematical model of the auroral boundaries, + Geomagnetism and Aeronomy, English Translation, 34(3), 331-336. + +""" + +import numpy as np + + +def starkov_auroral_boundary(mlt, al=-1, bnd='ocb'): + """Calculate the location of the auroral boundaries. + + Parameters + ---------- + mlt : float or array-like + Magnetic local time in hours + al : float or int + AL geomagnetic index, Auroral Electrojet Lower envelope in nT + (default=-1) + bnd : str + Boundary to calculate, expects one of 'ocb', 'eab', or 'diffuse' + (default='ocb') + + 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. + + References + ---------- + [8]_ + + """ + + # Calculate the AL dependence of the coefficients + A0 = starkov_coefficient_values(al, "A0", bnd) + A1 = starkov_coefficient_values(al, "A1", bnd) + alpha1 = starkov_coefficient_values(al, "alpha1", bnd) + A2 = starkov_coefficient_values(al, "A2", bnd) + alpha2 = starkov_coefficient_values(al, "alpha2", bnd) + A3 = starkov_coefficient_values(al, "A3", bnd) + alpha3 = starkov_coefficient_values(al, "alpha3", bnd) + + # Calculate the angular inputs in radians + rad1 = np.radians(15.0 * (mlt + alpha1)) + rad2 = np.radians(15.0 * (2.0 * mlt + alpha2)) + rad3 = np.radians(15.0 * (3.0 * mlt + alpha3)) + + # Calculate the boundary location in degrees latitude + bnd_lat = A0 + A1 * np.cos(rad1) + A2 * np.cos(rad2) + A3 * np.cos(rad3) + + # Ensure all co-latitudes are positive or zero + if np.asarray(bnd_lat).shape == (): + if bnd_lat < 0: + bnd_lat = 0.0 + else: + bnd_lat[bnd_lat < 0] = 0.0 + + return bnd_lat + + +def starkov_coefficient_values(al, coeff_name, bnd): + """Calculate the Starkov auroral model coefficient values. + + Parameters + ---------- + al : float or int + AL geomagnetic index, Auroral Electrojet Lower envelope in nT + coeff_name : str + Coefficient name, expects one of 'A0', 'A1', 'A2', 'A3', 'alpha1', + 'alpha2', or 'alpha3' + bnd : str + Boundary to calculate, expects one of 'ocb', 'eab', or 'diffuse' + + Returns + ------- + coeff : float + Coefficient value in hours (alpha) or degrees latitude (A) + + References + ---------- + [8]_ + + """ + # Define the model coefficients for each type of boundary + coeff_terms = {'A0': {'ocb': [-.07, 24.54, -12.53, 2.15], + 'eab': [1.16, 23.21, -10.97, 2.03], + 'diffuse': [3.44, 29.77, -16.38, 3.35]}, + 'A1': {'ocb': [-10.06, 19.83, -9.33, 1.24], + 'eab': [-9.59, 17.78, -7.20, 0.96], + 'diffuse': [-2.41, 7.89, -4.32, 0.87]}, + 'alpha1': {'ocb': [-6.61, 10.17, -5.80, 1.19], + 'eab': [-2.22, 1.50, -0.58, 0.08], + 'diffuse': [-1.68, -2.48, 1.58, -0.28]}, + 'A2': {'ocb': [-4.44, 7.47, -3.01, 0.25], + 'eab': [-12.07, 17.49, -7.96, 1.15], + 'diffuse': [-0.74, 3.94, -3.09, 0.72]}, + 'alpha2': {'ocb': [6.37, -1.10, 0.34, -0.38], + 'eab': [-23.98, 42.79, -26.96, 5.56], + 'diffuse': [8.69, -20.73, 13.03, -2.14]}, + 'A3': {'ocb': [-3.77, 7.90, -4.73, 0.91], + 'eab': [-6.56, 11.44, -6.73, 1.31], + 'diffuse': [-2.12, 3.24, -1.67, 0.31]}, + 'alpha3': {'ocb': [-4.48, 10.16, -5.87, 0.98], + 'eab': [-20.07, 36.67, -24.20, 5.11], + 'diffuse': [8.61, -5.34, -1.36, 0.76]}} + + # Calculate the desired coefficient + log_al = np.log10(abs(al)) + coeff = coeff_terms[coeff_name][bnd][0] + coeff_terms[coeff_name][bnd][ + 1] * log_al + coeff_terms[coeff_name][bnd][2] * ( + log_al**2) + coeff_terms[coeff_name][bnd][3] * (log_al**3) + + return coeff diff --git a/ocbpy/tests/test_boundary_dual.py b/ocbpy/tests/test_boundary_dual.py index 6706e5d3..5982e5ae 100644 --- a/ocbpy/tests/test_boundary_dual.py +++ b/ocbpy/tests/test_boundary_dual.py @@ -8,12 +8,14 @@ # ----------------------------------------------------------------------------- """Tests the boundary DualBoundary class.""" +import datetime import numpy from os import path import sys import unittest import ocbpy +from ocbpy.boundaries import models import ocbpy.tests.class_common as cc import ocbpy.tests.test_boundary_ocb as test_ocb @@ -79,14 +81,16 @@ def setUp(self): "r_err", "fom"], "ampere": ["date", "time", "x", "y", "fom"], "dmsp-ssj": ["date", "time", "sc", "x", "y", "fom", - "x_1", "x_2", "y_1", "y_2"]} + "x_1", "x_2", "y_1", "y_2"], + "model": ["fom", "rfunc_kwargs", "rfunc"]} self.not_attrs = {"image": ["date", "time", "x", "y", "x_1", "x_2", "y_1", "y_2", "sc"], "ampere": ["year", "soy", "x_1", "y_1", "x_2", "y_2", "sc", "num_sectors", "a", "r_err"], "dmsp-ssj": ["year", "soy", "num_sectors", "a", - "r_err"]} + "r_err"], + "model": ["date", "time", "year", "soy", "x", "y"]} self.inst_init = [{"eab_instrument": "dmsp-ssj", "hemisphere": 1, "eab_filename": path.join(cc.test_dir, "dmsp-ssj_north_out.eab"), @@ -104,7 +108,17 @@ def setUp(self): "dmsp-ssj_south_out.eab"), "ocb_instrument": "ampere", "ocb_filename": path.join(cc.test_dir, - "test_south_ocb")}] + "test_south_ocb")}, + {"ocb_instrument": "model", "hemisphere": 1, + "eab_instrument": "model", + "ocb_rfunc": models.starkov_auroral_boundary, + "eab_rfunc": models.starkov_auroral_boundary, + "stime": [datetime.datetime(2001, 1, 1, i) + for i in range(3)], + "ocb_rfunc_kwargs": {"al": [-50, -100, -300], + "bnd": "ocb"}, + "eab_rfunc_kwargs": {"al": [-50, -100, -300], + "bnd": "diffuse"}}] return diff --git a/ocbpy/tests/test_boundary_eab.py b/ocbpy/tests/test_boundary_eab.py index 40aa1eb6..e0284917 100644 --- a/ocbpy/tests/test_boundary_eab.py +++ b/ocbpy/tests/test_boundary_eab.py @@ -8,10 +8,12 @@ # ----------------------------------------------------------------------------- """Tests the boundary EABoundary class.""" +import datetime import numpy from os import path import ocbpy +from ocbpy.boundaries import models import ocbpy.tests.class_common as cc import ocbpy.tests.test_boundary_ocb as test_ocb @@ -43,11 +45,13 @@ def setUp(self): self.inst_attrs = {"image": ["year", "soy", "num_sectors", "a", "r_err", "fom"], "dmsp-ssj": ["date", "time", "sc", "x", "y", "fom", - "x_1", "x_2", "y_1", "y_2"]} + "x_1", "x_2", "y_1", "y_2"], + "model": ["fom", "rfunc_kwargs", "rfunc"]} self.not_attrs = {"image": ["date", "time", "x", "y", "x_1", "x_2", "y_1", "y_2", "sc"], "dmsp-ssj": ["year", "soy", "num_sectors", "a", - "r_err"]} + "r_err"], + "model": ["date", "time", "year", "soy", "x", "y"]} self.inst_init = [{"instrument": "image", "hemisphere": 1, "filename": path.join(cc.test_dir, "test_north_ocb")}, @@ -56,7 +60,13 @@ def setUp(self): "dmsp-ssj_north_out.eab")}, {"instrument": "dmsp-ssj", "hemisphere": -1, "filename": path.join(cc.test_dir, - "dmsp-ssj_south_out.eab")}] + "dmsp-ssj_south_out.eab")}, + {"instrument": "model", "hemisphere": 1, + "rfunc": models.starkov_auroral_boundary, + "stime": [datetime.datetime(2001, 1, 1, i) + for i in range(3)], + "rfunc_kwargs": {"al": [-50, -100, -300], + "bnd": "eab"}}] return diff --git a/ocbpy/tests/test_boundary_ocb.py b/ocbpy/tests/test_boundary_ocb.py index 3915146f..467ed8c1 100644 --- a/ocbpy/tests/test_boundary_ocb.py +++ b/ocbpy/tests/test_boundary_ocb.py @@ -14,6 +14,7 @@ import unittest import ocbpy +from ocbpy.boundaries import models import ocbpy.tests.class_common as cc @@ -126,14 +127,16 @@ def setUp(self): "r_err", "fom"], "ampere": ["date", "time", "x", "y", "fom"], "dmsp-ssj": ["date", "time", "sc", "x", "y", "fom", - "x_1", "x_2", "y_1", "y_2"]} + "x_1", "x_2", "y_1", "y_2"], + "model": ["fom", "rfunc_kwargs", "rfunc"]} self.not_attrs = {"image": ["date", "time", "x", "y", "x_1", "x_2", "y_1", "y_2", "sc"], "ampere": ["year", "soy", "x_1", "y_1", "x_2", "y_2", "sc", "num_sectors", "a", "r_err"], "dmsp-ssj": ["year", "soy", "num_sectors", "a", - "r_err"]} + "r_err"], + "model": ["date", "time", "year", "soy", "x", "y"]} self.inst_init = [{"instrument": "image", "hemisphere": 1, "filename": path.join(cc.test_dir, "test_north_ocb")}, @@ -145,7 +148,13 @@ def setUp(self): "dmsp-ssj_south_out.ocb")}, {"instrument": "ampere", "hemisphere": -1, "filename": path.join(cc.test_dir, - "test_south_ocb")}] + "test_south_ocb")}, + {"instrument": "model", "hemisphere": 1, + "rfunc": models.starkov_auroral_boundary, + "stime": [datetime.datetime(2001, 1, 1, i) + for i in range(3)], + "rfunc_kwargs": {"al": [-50, -100, -300], + "bnd": "ocb"}}] self.ocb = None return diff --git a/ocbpy/tests/test_models.py b/ocbpy/tests/test_models.py new file mode 100644 index 00000000..a437d792 --- /dev/null +++ b/ocbpy/tests/test_models.py @@ -0,0 +1,99 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# DOI: 10.5281/zenodo.1179230 +# Full license can be found in License.md +# ----------------------------------------------------------------------------- +"""Tests the boundaries.models functions.""" + +import numpy as np +import unittest + +from ocbpy.boundaries import models + + +class TestStarkovModel(unittest.TestCase): + """"Unit tests for the Starkov 1994 routines.""" + + def setUp(self): + """Initialize the test case by copying over necessary files.""" + self.mlt = np.arange(0, 24, 1) + self.coeff_out = { + 'A0': {'ocb': -.07, 'eab': 1.16, 'diffuse': 3.44}, + 'A1': {'ocb': -10.06, 'eab': -9.59, 'diffuse': -2.41}, + 'alpha1': {'ocb': -6.61, 'eab': -2.22, 'diffuse': -1.68}, + 'A2': {'ocb': -4.44, 'eab': -12.07, 'diffuse': -0.74}, + 'alpha2': {'ocb': 6.37, 'eab': -23.98, 'diffuse': 8.69}, + 'A3': {'ocb': -3.77, 'eab': -6.56, 'diffuse': -2.12}, + 'alpha3': {'ocb': -4.48, 'eab': -20.07, 'diffuse': 8.61}} + self.al = [-1, -500] + self.max_lat = {'ocb': [11.66543814, 19.2853977], + 'eab': [24.1150303, 29.14429888], + 'diffuse': [7.00627994, 34.72168093]} + return + + def tearDown(self): + """Clean up the test environment.""" + del self.mlt, self.coeff_out, self.al, self.max_lat + return + + def test_coeff_construction(self): + """Test coefficient calculation for an AL of -1.""" + + for coeff in self.coeff_out.keys(): + for bnd in self.coeff_out[coeff].keys(): + with self.subTest(coeff=coeff, bnd=bnd): + # Calculate the coefficient value + out = models.starkov_coefficient_values( + self.al[0], coeff, bnd) + + # Compare the output + self.assertEqual(out, self.coeff_out[coeff][bnd]) + return + + def test_coeff_bad_coeff(self): + """Test a KeyError is raised for an unknown coeffcient name.""" + coeff = "not a coefficient" + with self.assertRaisesRegex(KeyError, coeff): + models.starkov_coefficient_values(self.al[0], coeff, + list(self.max_lat.keys())[0]) + 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.starkov_coefficient_values( + self.al[0], list(self.coeff_out.keys())[0], bound) + return + + def test_bound_loc_array(self): + """Test the expected boundary location across an MLT array.""" + # Cycle through low and high AL values + for ia, in_al in enumerate(self.al): + for bnd in self.max_lat.keys(): + with self.subTest(al=in_al, bnd=bnd): + lat = models.starkov_auroral_boundary( + self.mlt, al=in_al, bnd=bnd) + + # Test the output latitude shape and values + self.assertTupleEqual(self.mlt.shape, lat.shape) + self.assertGreaterEqual(min(lat), 0) + self.assertAlmostEqual(max(lat), self.max_lat[bnd][ia]) + + return + + def test_bound_loc_float(self): + """Test the expected boundary location across an MLT value.""" + # Cycle through low and high AL values + for ia, in_al in enumerate(self.al): + for bnd in self.max_lat.keys(): + with self.subTest(al=in_al, bnd=bnd): + lat = models.starkov_auroral_boundary( + self.mlt[0], al=in_al, bnd=bnd) + + # Test the output latitude shape and values + self.assertTrue(isinstance(lat, float)) + self.assertGreaterEqual(lat, 0) + self.assertLessEqual(lat, self.max_lat[bnd][ia]) + + return