ENH: Improve calc_measured_loc performance and support

CHANGELOG.md

@@ -2,6 +2,12 @@
 All notable changes to this project will be documented in this file.
 This project adheres to [Semantic Versioning](https://semver.org/).
 
+## [0.1.0] - 2021-XX-XX
+- Documentation
+   - Added examples for JRO methods
+- Testing
+   - Added unit tests for JRO methods
+
 ## [0.0.4] - 2021-06-11
 - Made changes to structure to comply with updates in pysat 3.0.0
 - Deprecations

docs/examples.rst

@@ -6,3 +6,4 @@ tools
 
 .. toctree::
    examples/ex_init.rst
+   examples/ex_jro_isr_beam.rst

docs/examples/ex_jro_isr_beam.rst

@@ -0,0 +1,86 @@
+.. _ex-jro-beam-loc:
+
+Calculate JRO ISR Beam Location
+===============================
+
+For measurements made using a single beam, it may be more appropriate to account
+for the changes in beam direction with altitude when determining the measurement
+location instead of using the location of the radar.  The method
+:py:meth:`instruments.methods.jro.calc_measurement_loc` (see :ref:`meth-jro`)
+uses the beam azimuth and elevation measurements to determine the geodetic
+latitude and longitude.
+
+This method is designed to be used with the JRO ISR data, and so assumes the
+azimuths and elevations have data variable names with the format ``'eldir#'``
+and ``'azdir#'`` (where **#** is the beam number), or ``'elm'`` and ``'azm'``.
+It will modify the :py:attr:`pysat.Instrument.data` object by adding latitude
+(``'gdlat#'``) and longitude (``'gdlon#'``) variables for every beam that has
+appropriately labeled azimuth and elevatiton data.  If the azimuth and elevation
+angle variables don't specify the beam number, **#** will be set to ``'_bm'``.
+
+The easiest way to use :py:meth:`instruments.methods.jro.calc_measurement_loc`
+is to attach it to the JRO ISR :py:class:`pysat.Instrument` as a `custom
+function <https://pysat.readthedocs.io/en/latest/tutorial/tutorial_custom.html>`_
+before loading data.
+
+.. code::
+
+   import datetime as dt
+   import pysat
+   import pysatMadrigal as pysat_mad
+
+   jro_obl = pysat.Instrument(inst_module=pysat_mad.instruments.jro_isr,
+                              tag='oblique_long')
+   jro_obl.custom_attach(pysat_mad.instruments.methods.jro.calc_measurement_loc)
+
+
+If necessary, download the desired data before loading it.  The geographic
+beam locations will be present alongside the azimuths and elevations.
+
+.. code::
+
+   ftime = dt.datetime(2010, 4, 12)
+
+   if not ftime in jro_obl.files.files.index:
+       jro_obl.download(start=ftime)
+
+   jro_obl.load(date=ftime)
+   'gdlat_bm' in jro_obl.variables and 'gdlon_bm' in jro_obl.variables
+
+
+The result of the above command should be ``True``.  To better visualize the
+beam location calculation, let us plot the locations of the beam range gates
+and the radar location.
+
+.. code::
+
+   import matplotlib.pyplot as plt
+
+   # Initialize the figure and axes
+   fig = plt.figure()
+   ax_alt = fig.add_subplot(211)
+   ax_geo = fig.add_subplot(212)
+
+   # Plot the altitude data
+   ax_alt.plot(jro_obl['gdlatr'], 0.52, 'X', color='green')
+   ax_alt.plot(jro_obl['gdlat_bm'], jro_obl['gdalt'], 'P', color='springgreen')
+
+   # Plot the lat/lon data
+   ax_geo.plot(jro_obl['gdlatr'], jro_obl['gdlonr'], 'X', color='green')
+   ax_geo.plot(jro_obl['gdlat_bm'], jro_obl['gdlon_bm'], 'P',
+               color='springgreen')
+
+   # Format the figure
+   ax_geo.set_xlabel('Latitude ($\circ$)')
+   ax_geo.set_ylabel('Longitude ($\circ$)')
+   ax_alt.set_ylabel('Altitude (km)')
+   ax_alt.legend(['ISR Location', 'ISR Beam'], fontsize='medium')
+   fig.suptitle('{:s} {:s} {:s} data at {:s}\n`pysatMadrigal.instruments.method.jro.calc_measurement_loc` results'.format(
+       jro_obl.platform, jro_obl.name, jro_obl.tag,
+       jro_obl.index[0].strftime('%d %b %Y')), fontsize='medium')
+
+
+.. image:: ../figures/ex_jro_isr_beam.png
+    :width: 800px
+    :align: center
+    :alt: Beam location in altiutde (top) and longitude (bottom) as a function of latitude along with the radar location.

docs/methods.rst

@@ -1,10 +1,14 @@
+.. _methods:
+
 Methods
 =======
 
 Several methods exist to help combine multiple data sets and convert between
 equivalent indices.
 
 
+.. _meth-dmsp:
+
 DMSP
 ----
 
@@ -15,6 +19,9 @@ common custom routines alongside reference and acknowledgement information.
 .. automodule:: pysatMadrigal.instruments.methods.dmsp
    :members:
 
+
+.. _meth-gnss:
+
 GNSS
 ----
 
@@ -25,6 +32,9 @@ reference and acknowledgement information.
 .. automodule:: pysatMadrigal.instruments.methods.gnss
    :members:
 
+
+.. _meth-jro:
+
 JRO
 ---
 
@@ -35,6 +45,9 @@ routines alongside reference and acknowledgement information.
 .. automodule:: pysatMadrigal.instruments.methods.jro
    :members:
 
+
+.. _meth-gen:
+
 General
 -------
 Supports the Madrigal data access.

pysatMadrigal/instruments/methods/jro.py

@@ -52,42 +52,79 @@ def calc_measurement_loc(inst):
     inst : pysat.Instrument
         JRO ISR Instrument object
 
+    Raises
+    ------
+    ValueError
+        If no appropriate azimuth and elevation angles are found, if no range
+        variable is found, or if multiple range variables are found.
+
+
     Note
     ----
     Adds 'gdlat#', 'gdlon#' to the instrument, for all directions that
     have azimuth and elevation keys that match the format 'eldir#' and 'azdir#'
+    or 'azm' and 'elm' (in this case # will be replaced with '_bm')
 
     """
-
-    az_keys = [kk[5:] for kk in list(inst.data.keys())
-               if kk.find('azdir') == 0]
-    el_keys = [kk[5:] for kk in list(inst.data.keys())
-               if kk.find('eldir') == 0]
     good_dir = list()
+    good_pre = list()
+
+    # Assume the 'dir#' format is used
+    az_keys = [kk[5:] for kk in inst.variables if kk.find('azdir') == 0]
+    el_keys = [kk[5:] for kk in inst.variables if kk.find('eldir') == 0]
 
     for i, kk in enumerate(az_keys):
         if kk in el_keys:
             try:
-                good_dir.append(int(kk))
+                good_dir.append("{:d}".format(int(kk)))
+                good_pre.append('dir{:s}'.format(kk))
             except ValueError:
-                logger.warning("unknown direction number [{:}]".format(kk))
+                logger.warning("Unknown direction number [{:}]".format(kk))
 
-    # Calculate the geodetic latitude and longitude for each direction
+    # Assume the 'm' format is used
+    if 'azm' in inst.variables and 'elm' in inst.variables:
+        good_dir.append('_bm')
+        good_pre.append('m')
+
+    # Test the success of finding the azimuths and elevations
     if len(good_dir) == 0:
         raise ValueError("No matching azimuth and elevation data included")
 
-    for dd in good_dir:
+    # Set common meta data variables. Includes determining the longitude range,
+    # which is only possible because JRO is in the western hemisphere.
+    lon_min = 0.0 if inst['gdlonr'] > 0.0 else -180.0
+    lon_max = 360.0 + lon_min
+    notes = 'Calculated using {:s}'.format(__name__)
+
+    # Get the range key
+    range_data = None
+    for rkey in ['range', 'rgate']:
+        if rkey in inst.variables:
+            if range_data is None:
+                if rkey == 'rgate':
+                    range_data = inst['gdalt']
+                else:
+                    range_data = inst[rkey]
+            else:
+                raise ValueError('Multiple range variables found')
+
+    if range_data is None:
+        raise ValueError('No range variable found')
+
+    # Calculate the geodetic latitude and longitude for each direction
+    for i, dd in enumerate(good_dir):
         # Format the direction location keys
-        az_key = 'azdir{:d}'.format(dd)
-        el_key = 'eldir{:d}'.format(dd)
-        lat_key = 'gdlat{:d}'.format(dd)
-        lon_key = 'gdlon{:d}'.format(dd)
+        az_key = 'az{:s}'.format(good_pre[i])
+        el_key = 'el{:s}'.format(good_pre[i])
+        lat_key = 'gdlat{:s}'.format(dd)
+        lon_key = 'gdlon{:s}'.format(dd)
+
         # JRO is located 520 m above sea level (jro.igp.gob.pe./english/)
         # Also, altitude has already been calculated
         gdaltr = np.ones(shape=inst['gdlonr'].shape) * 0.52
         gdlat, gdlon, _ = coords.local_horizontal_to_global_geo(inst[az_key],
                                                                 inst[el_key],
-                                                                inst['range'],
+                                                                range_data,
                                                                 inst['gdlatr'],
                                                                 inst['gdlonr'],
                                                                 gdaltr,
@@ -98,16 +135,20 @@ def calc_measurement_loc(inst):
         inst.data = inst.data.assign({lat_key: gdlat, lon_key: gdlon})
 
         # Add metadata for the new data values
-        bm_label = "Beam {:d} ".format(dd)
+        bm_label = "Beam" if dd[0] == "_" else "Beam {:s} ".format(dd)
         inst.meta[lat_key] = {inst.meta.labels.units: 'degrees',
                               inst.meta.labels.name: bm_label + 'latitude',
+                              inst.meta.labels.notes: notes,
                               inst.meta.labels.desc: bm_label + 'latitude',
                               inst.meta.labels.min_val: -90.0,
                               inst.meta.labels.max_val: 90.0,
                               inst.meta.labels.fill_val: np.nan}
         inst.meta[lon_key] = {inst.meta.labels.units: 'degrees',
+                              inst.meta.labels.notes: notes,
                               inst.meta.labels.name: bm_label + 'longitude',
                               inst.meta.labels.desc: bm_label + 'longitude',
+                              inst.meta.labels.min_val: lon_min,
+                              inst.meta.labels.max_val: lon_max,
                               inst.meta.labels.fill_val: np.nan}
 
     return