Merge pull request #134 from aburrell/geo_coords

Add geographic coordinate options

.github/workflows/main.yml

@@ -13,10 +13,12 @@ jobs:
       fail-fast: false
       matrix:
         os: ["ubuntu-latest", "macos-latest", "windows-latest"]
-        python-version: ["3.7", "3.8", "3.9", "3.10", "3.11"]
+        python-version: ["3.9", "3.10", "3.11", "3.12"]
         install-extras: [0, 1, 2]  # TODO(#129): update to replace extra flag
         exclude:
-          - os: "windows-latest"
+          - os: "windows-latest"  # TODO(#135): remove exclusion
+            install-extras: 2
+          - os: "macos-latest"  # TODO(#135): remove exclusion
             install-extras: 2
 
     name: Python ${{ matrix.python-version }} on ${{ matrix.os }} with extras ${{ matrix.install-extras }}

Changelog.rst

@@ -7,6 +7,7 @@ Summary of all changes made since the first stable release
 ------------------
 * DEP: Deprecated functions that depend on ssj_auroral_boundary_package
 * DEP: Removed classes and kwargs deprecated in v0.3.0
+* MAINT: Cycled supported Python versions
 * MAINT: Added a pyproject.toml file and removed setup.py
 * MAINT: Updated numpy logic to address DeprecationWarnings
 * MAINT: Updated GitHub Action yamls to use pyproject.toml and cycle versions
@@ -15,12 +16,17 @@ Summary of all changes made since the first stable release
 * ENH: Added AMPERE EABs, using the Heppner-Maynard boundary as a valid EAB
 * ENH: Changed default directory ID to use `pathlib`
 * ENH: Allow data padding in `pysat_instrument` functions
+* ENH: Created separate vector transformation functions to support multiple
+       coordinate systems
+* ENH: Updated VectorData and pysat_instruments to allow geographic inputs
 * BUG: Fixed a typo in the documentation's pysat example
 * BUG: Added an error catch for badly formatted SuperMAG file reading
+* BUG: Fixed the flat/zero masking for vector pole angles with array input
 * TST: Added a new CI test for the Test PyPi Release Candidate
 * TST: Reduced duplication by creating a common test class and test variable
 * TST: Added a ReadTheDocs yaml
 * DOC: Improved docstring style compliance and correctness
+* DOC: Updated pysat example to use geodetic inputs and improve the EAB plot
 
 0.3.0 (10-21-2022)
 ------------------

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.7 - 3.10.
+This module currently supports Python version 3.9 - 3.12.
 
 # Dependencies
 

docs/examples/ex_pysat_eab.rst

@@ -88,96 +88,81 @@ using the following commands.
 prodedures while loading the desired data.  The
 :py:mod:`ocbpy.instrument.pysat_instrument` module contains functions that may
 be applied using the :py:mod:`pysat` `custom interface
-<https://pysat.readthedocs.io/en/latest/tutorial/tutorial_custom.html>`_.
-However, before this can be done the magnetic locations need to be calculated.
-This can be done by writing an appropriate function that takes the
-:py:class:`pysat.Instrument` object as input and updates it within the function.
-
+<https://pysat.readthedocs.io/en/latest/tutorial/tutorial_custom.html>`_.  The
+EAB conversion can handle magnetic, geodetic, or geocentric coordinates for
+scalar or vector data types using the :py:var:`loc_coord` and
+:py:var:`vect_coord` keyword arguments, respectively. We do need to calculate
+local time before calculating the EAB coordinates, though.
 
 ::
 
-   import aacgmv2
-   import numpy as np
-
-   def add_mag_coords(inst, lat='gdlat', lon='glon', alt='gdalt'):
-       """Add AACGMV2 magnetic coordinates.
+   def add_slt(inst, lon='glon', slt='slt'):
+       """Add solar local time.
 
        Parameters
        ----------
        inst : pysat.Instrument
            Data object
-       lat : str
-           Geodetic latitude key (default='gdlat')
        lon : str
            Geographic longitude key (default='glon')
-       alt : str
-           Geodetic altitude key (default='gdalt')
+       slt : str
+           Key for new solar local tim data (default='slt')
+
        """
        # Initalize the data arrays
-       mlat = np.full(shape=tuple(tec.data.dims[kk]
-                                  for kk in ['time', lat, lon]),
-                      fill_value=np.nan)
-       mlt = np.full(shape=mlat.shape, fill_value=np.nan)
-
-       # Cycle through all times, calculating magnetic locations
-       for i, dtime in enumerate(inst.index):
-           for j, gdlat in enumerate(inst[lat].values):
-               height = inst[alt][i, j].values
-               if not np.isnan(height).all():
-                   mlat[i, j], mlon, r = aacgmv2.convert_latlon_arr(
-                       gdlat, inst[lon].values, height, dtime)
-                   mlt[i, j] = aacgmv2.convert_mlt(mlon, dtime)
-
-       # Assign the magnetic data to the input Instrument
-       inst.data = inst.data.assign({"mlat": (("time", lat, lon), mlat),
-                                     "mlt": (("time", lat, lon), mlt)})
+       lt = [ocbpy.ocb_time.glon2slt(inst[lon], dtime) for dtime in inst.index]
+
+       # Assign the SLT data to the input Instrument
+       inst.data = inst.data.assign({slt: (("time", lon), lt)})
        return
 
     
-Assign this function and the ocbpy function in the desired order of operations.
+Assign this function and the :py:mod:`ocbpy` function in the desired order of
+operations. When calculating magnetic coordinates, it is important to specify
+the height, which can be a single value or specified for each observation.
 
 ::
 
-   
-   tec.custom_attach(add_mag_coords)
+   tec.custom_attach(add_slt, kwargs={'lon': 'glon'})
    tec.custom_attach(ocbpy.instruments.pysat_instruments.add_ocb_to_data,
-                     kwargs={'ocb': eab, 'mlat_name': 'mlat',
-                     'mlt_name': 'mlt', 'max_sdiff': 150})
+                     kwargs={'ocb': eab, 'mlat_name': 'gdlat',
+                     'mlt_name': 'slt', 'height_name': 'gdalt',
+                     'loc_coord': 'geodetic', 'max_sdiff': 150})
    tec.load(date=eab.dtime[eab.rec_ind])
    print(tec.variables)
 
-   ['time', 'gdlat', 'glon', 'dtec', 'gdalt', 'tec', 'mlat', 'mlt', 'mlat_ocb',
-    'mlt_ocb', 'r_corr_ocb']
+   ['time', 'gdlat', 'glon', 'dtec', 'gdalt', 'tec', 'gdlat_eab',
+    'slt_eab', 'r_corr_eab']
 
 Now we have the EAB coordinates for each location in the Northern Hemisphere
 where a good EAB was found within 2.5 minutes of the data record.  This time
 difference was chosen because the VTEC data has a 5 minute resolution.
 
-Now, let's plot the average of the VTEC poleward of the EAB. To do this we will
-first need to calculate these averages.
+Now, let's plot the average of the VTEC equatorward of the EAB. To do this we
+will first need to calculate these averages.
 
 ::
 
    del_lat = 2.0
    del_mlt = 2.0
-   ave_lat = np.arange(eab.boundary_lat, 90, del_lat)
+   ave_lat = np.arange(45, eab.boundary_lat, del_lat)
    ave_mlt = np.arange(0, 24, del_mlt)
    ave_tec = np.full(shape=tec['tec'].shape, fill_value=np.nan)
 
    for lat in ave_lat:
        for mlt in ave_mlt:
            # We are not overlapping bins, so don't need to worry about MLT
-	   # rollover from 0-24
+           # rollover from 0-24
            sel_tec = tec['tec'].where(
-               (tec['mlat_ocb'] > lat) & (tec['mlat_ocb'] <= lat + del_lat)
-               & (tec['mlt_ocb'] >= mlt) & (tec['mlt_ocb'] < mlt + del_mlt))
+               (tec['gdlat_eab'] > lat) & (tec['gdlat_eab'] <= lat + del_lat)
+               & (tec['slt_eab'] >= mlt) & (tec['slt_eab'] < mlt + del_mlt))
            inds = np.where(~np.isnan(sel_tec.values))
            if len(inds[0]) > 0:
                ave_tec[inds] = np.nanmean(sel_tec.values)
 
 Now let us plot these averages at the EAB location of each measurement.  This
 will provide us with knowledge of the coverage as well as knowledge of the
-average behaviour.
+average behaviour of the sub-auroral VTEC on this day.
 
 ::
 
@@ -200,11 +185,12 @@ average behaviour.
    ax.plot(lon, lat, "m-", linewidth=2, label="EAB")
 
    # Plot the VTEC data
-   tec_lon = tec['mlt_ocb'].values * np.pi / 12.0
-   tec_lat = 90.0 - tec['mlat_ocb'].values
+   tec_lon = (tec['slt_eab'].values * np.pi / 12.0)
+   tec_lat = (90.0 - tec['gdlat_eab'].values)
    tec_max = np.ceil(np.nanmax(ave_tec))
-   con = ax.scatter(tec_lon, tec_lat, c=ave_tec, marker="s",
-                    cmap=mpl.cm.get_cmap("viridis"), s=5, vmin=0, vmax=tec_max)
+   con = ax.scatter(tec_lon, tec_lat, c=ave_tec.transpose([0, 2, 1]),
+                    marker="s", cmap=mpl.colormaps["viridis"], s=5,
+                    vmin=0, vmax=tec_max)
 
    # Add a colourbar and labels
    tticks = np.linspace(0, tec_max, 6, endpoint=True)

ocbpy/__init__.py

@@ -6,11 +6,7 @@
 """Auroral oval and polar cap normalised location calculation tools."""
 
 import logging
-
-try:
-    from importlib import metadata
-except ImportError:
-    import importlib_metadata as metadata
+from importlib import metadata
 
 # Define a logger object to allow easier log handling
 logging.raiseExceptions = False
@@ -25,6 +21,7 @@
 from ocbpy import ocb_correction  # noqa F401
 from ocbpy import ocb_scaling  # noqa F401
 from ocbpy import ocb_time  # noqa F401
+from ocbpy import vectors  # noqa F401
 
 from ocbpy._boundary import DualBoundary  # noqa F401
 from ocbpy._boundary import EABoundary  # noqa F401

ocbpy/_boundary.py

@@ -861,7 +861,7 @@ def _set_default_rfunc(self):
             self.rfunc = np.full(shape=self.records,
                                  fill_value=ocbcor.elliptical)
         else:
-            raise ValueError("unknown instrument")
+            raise ValueError("unknown instrument: {:}".format(self.instrument))
 
         return
 

ocbpy/boundaries/dmsp_ssj_files.py

@@ -25,18 +25,30 @@
 import warnings
 import zipfile
 
+import aacgmv2
+
 import ocbpy
 
-err = ''.join(['unable to load the DMSP SSJ module; ssj_auroral_boundary ',
-               'is available at: ',
-               'https://github.com/lkilcommons/ssj_auroral_boundary'])
 try:
     from spacepy import pycdf
-    import aacgmv2
     import ssj_auroral_boundary as ssj
+    err = ''
+except ImportError as ierr:
+    ssj = None
+    err = ''.join(['unable to load the DMSP SSJ module; ssj_auroral_boundary ',
+                   'is available at: ',
+                   'https://github.com/lkilcommons/ssj_auroral_boundary\n',
+                   str(ierr)])
+
+try:
     import zenodo_get
 except ImportError as ierr:
-    raise ImportError("{:s}\n{:}".format(err, ierr))
+    zenodo_get = None
+    err = ''.join([err, '\nunable to load `zenodo_get` module; avalable ',
+                   'from PyPi.\n', str(ierr)])
+
+if ssj is None and zenodo_get is None:
+    raise ImportError(err)
 
 
 def fetch_ssj_files(stime, etime, out_dir=None, sat_nums=None):
@@ -77,6 +89,11 @@ def fetch_ssj_files(stime, etime, out_dir=None, sat_nums=None):
         ssj_auroral_boundaries package is no longer supported; use
         `fetch_ssj_boundary_files`
 
+    Raises
+    ------
+    ImportError
+        If called and ssj_auroral_boundaries is not available
+
     See Also
     ---------
     requests.exceptions.ProxyError
@@ -88,6 +105,10 @@ def fetch_ssj_files(stime, etime, out_dir=None, sat_nums=None):
                            " will be removed in version 0.4.1+."]),
                   DeprecationWarning, stacklevel=2)
 
+    if ssj is None:
+        raise ImportError(
+            'depends on uninstalled package ssj_auroral_boundaries')
+
     # Get and test the output directory
     if out_dir is None:
         out_dir = ocbpy.boundaries.files.get_boundary_directory()
@@ -135,8 +156,8 @@ def fetch_ssj_files(stime, etime, out_dir=None, sat_nums=None):
                 try:
                     ssj.files.download_cdf_from_noaa(remote, local)
                     out_files.append(local)
-                except RuntimeError as err:
-                    ocbpy.logger.info(err)
+                except RuntimeError as rerr:
+                    ocbpy.logger.info(rerr)
 
             # Cycle by one day
             ctime += dt.timedelta(days=1)
@@ -175,6 +196,8 @@ def create_ssj_boundary_files(cdf_files, out_dir=None,
     ------
     ValueError
         If incorrect input is provided
+    ImportError
+        If called and ssj_auroral_boundaries is not available
 
     Warnings
     --------
@@ -189,6 +212,10 @@ def create_ssj_boundary_files(cdf_files, out_dir=None,
                            " will be removed in version 0.4.1+."]),
                   DeprecationWarning, stacklevel=2)
 
+    if ssj is None:
+        raise ImportError(
+            'depends on uninstalled package ssj_auroral_boundaries')
+
     # Test the directory inputs
     if out_dir is None:
         out_dir = ocbpy.boundaries.files.get_boundary_directory()
@@ -225,8 +252,8 @@ def create_ssj_boundary_files(cdf_files, out_dir=None,
                                                  make_plot=make_plots,
                                                  csvvars=out_cols)
                     out_files.append(absd.csv.csvfn)
-            except pycdf.CDFError as err:
-                ocbpy.logger.warning("{:}".format(err))
+            except pycdf.CDFError as cerr:
+                ocbpy.logger.warning("{:}".format(cerr))
             except Warning as war:
                 ocbpy.logger.warning("{:}".format(war))
         else:
@@ -268,13 +295,17 @@ def fetch_ssj_boundary_files(stime=None, etime=None, out_dir=None,
         If an unknown satellite ID is provided.
     IOError
         If unable to donwload the target archive and identify the zip file
+    ImportError
+        If called and zenodo_get is not available
 
     Notes
     -----
     If a file already exists, the routine will add the file to the output list
     without downloading it again.
 
     """
+    if zenodo_get is None:
+        raise ImportError('depends on uninstalled package zenodo_get')
 
     # Test the requested satellite outputs. SSJ5 was carried on F16 onwards.
     # F19 was short lived, F20 was not launched. Ref: