Fixes for REI calculation, docs, and tests

 Changes to be committed:
	modified:   docs/index.rst
	modified:   docs/notebooks/interval_data.ipynb.rst
	modified:   docs/notebooks/receiver_efficiency_index.ipynb.rst
	modified:   py_notebooks/interval_data.ipynb
	modified:   py_notebooks/receiver_efficiency_index.ipynb
	modified:   resonate/filters.py
	modified:   resonate/receiver_efficiency.py
	modified:   tests/assertion_files/nsbs_interval.csv
	modified:   tests/interval_test.py

docs/index.rst

@@ -101,7 +101,7 @@ This residence index tool will take a compressed or uncompressed detection file
 Receiver Efficiency Index
 -------------------------
 
-*(Ellis, R., Flaherty-Walia, K., Collins, A., Bickford, J., Walters Burnsed,  Lowerre-Barbieri S. 2018. Acoustic telemetry array evolution: from species- and project-specific designs to large-scale, multispecies, cooperative networks)*
+`(Ellis, R., Flaherty-Walia, K., Collins, A., Bickford, J., Walters Burnsed,  Lowerre-Barbieri S. 2018. Acoustic telemetry array evolution: from species- and project-specific designs to large-scale, multispecies, cooperative networks) <https://doi.org/10.1016/j.fishres.2018.09.015>`_
 
 The receiver efficiency index is number between ``0`` and ``1`` indicating the amount of relative activity at each receiver compared to the entire set of receivers, regardless of positioning. The function takes a set detections and a deployment history of the receivers to create a context for the detections. Both the amount of unique tags and number of species are taken into consideration in the calculation. For the exact method, see the details in :ref:`Receiver Efficiency Index<receiver_efficiency_index_page>`.
 

docs/notebooks/interval_data.ipynb.rst

@@ -47,11 +47,11 @@ You can modify individual stations if needed by using
 
 .. code:: python
 
-    station_name = 'station'
+    station_name = 'HFX001'
     
     station_detection_radius = 500
     
-    station_det_radius.set_value(station_name, 'radius', geopy.distance.Distance( station_detection_radius/1000.0 ))
+    station_det_radius.at[station_name, 'radius'] = geopy.distance.Distance( station_detection_radius/1000.0 )
 
 Create the interval data by passing the compressed detections, the
 matrix, and the station radii.

docs/notebooks/receiver_efficiency_index.ipynb.rst

@@ -16,7 +16,7 @@ formula of:
 .. container:: large-math
 
    REI =
-   :math:`\frac{\left(\frac{T_r}{T_a} \times \frac{S_r}{S_a}\right) / \left(\frac{DD_a}{DD_r}\right)}{D_r}`
+   :math:`\frac{T_r}{T_a} \times \frac{S_r}{S_a} \times \frac{DD_r}{DD_a} \times \frac{D_a}{D_r}`
 
 .. raw:: html
 
@@ -31,6 +31,7 @@ formula of:
    receivers
 -  :math:`DD_r` = The number of unique days with detections on the
    receiver
+-  :math:`D_a` = The number of days the array was active
 -  :math:`D_r` = The number of days the receiver was active
 
 Each REI is then normalized against the sum of all considered stations.

py_notebooks/interval_data.ipynb

@@ -71,11 +71,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "station_name = 'station'\n",
+    "station_name = 'HFX001'\n",
     "\n",
     "station_detection_radius = 500\n",
     "\n",
-    "station_det_radius.set_value(station_name, 'radius', geopy.distance.Distance( station_detection_radius/1000.0 ))"
+    "station_det_radius.at[station_name, 'radius'] = geopy.distance.Distance( station_detection_radius/1000.0 )"
    ]
   },
   {
@@ -163,5 +163,5 @@
   }
  },
  "nbformat": 4,
- "nbformat_minor": 1
+ "nbformat_minor": 2
 }

py_notebooks/receiver_efficiency_index.ipynb

@@ -6,7 +6,9 @@
    "source": [
     "# Receiver Efficiency Index\n",
     "\n",
-    "The receiver efficiency index is number between ``0`` and ``1`` indicating the amount of relative activity at each receiver compared to the entire set of receivers, regardless of positioning. The function takes a set detections and a deployment history of the receivers to create a context for the detections. Both the amount of unique tags and number of species are taken into consideration in the calculation.\n",
+    "The receiver efficiency index is number between ``0`` and ``1`` indicating the amount of relative activity at each receiver compared to the entire set of receivers, regardless of positioning. \n",
+    "The function takes a set detections and a deployment history of the receivers to create a context for the detections. Both the amount of unique tags and number of species are taken into \n",
+    "consideration in the calculation.\n",
     "\n",
     "The receiver efficiency index implement is implemented based on the paper [paper place holder]. Each receiver's index is calculated on the formula of:\n",
     "\n",
@@ -15,7 +17,7 @@
     "\n",
     "<div class=\"large-math\">\n",
     "    \n",
-    "REI = $\\frac{\\left(\\frac{T_r}{T_a} \\times \\frac{S_r}{S_a}\\right) / \\left(\\frac{DD_a}{DD_r}\\right)}{D_r}$\n",
+    "REI = $\\frac{T_r}{T_a} \\times \\frac{S_r}{S_a} \\times \\frac{DD_r}{DD_a} \\times \\frac{D_a}{D_r}$\n",
     "\n",
     "</div>\n",
     "\n",
@@ -28,6 +30,7 @@
     "* $S_a$ = The number of species detected across all receivers\n",
     "* $DD_a$ = The number of unique days with detections across all receivers\n",
     "* $DD_r$ = The number of unique days with detections on the receiver\n",
+    "* $D_a$ = The number of days the array was active\n",
     "* $D_r$ = The number of days the receiver was active\n",
     "\n",
     "\n",

resonate/filters.py

@@ -2,7 +2,7 @@ from datetime import datetime, timedelta
 
 import numpy as np
 import pandas as pd
-from geopy.distance import vincenty
+from geopy.distance import geodesic
 from resonate.library.exceptions import GenericException
 
 
@@ -28,7 +28,7 @@ def get_distance_matrix(detections):
                       stn_locs.loc[cstation, 'longitude'])
             rpoint = (stn_locs.loc[rstation, 'latitude'],
                       stn_locs.loc[rstation, 'longitude'])
-            dist_mtx.loc[rstation, cstation] = vincenty(cpoint, rpoint).m
+            dist_mtx.loc[rstation, cstation] = geodesic(cpoint, rpoint).m
     dist_mtx.index.name = None
     return dist_mtx
 

resonate/receiver_efficiency.py

@@ -31,11 +31,12 @@ def REI(detections, deployments):
         # date columns
         deployments = deployments.copy(deep=True)
         detections = detections.copy(deep=True)
-        deployments['recovery_notes'] = deployments.recovery_date.str.extract(
-            '([A-Za-z\//:]+)', expand=False)
-        deployments.recovery_date = deployments.recovery_date.str.extract(
-            '(\d+-\d+-\d+)', expand=False)
-        deployments = deployments.replace('-', np.nan)
+        if deployments.recovery_date.dtype != np.dtype('<M8[ns]'):
+            deployments['recovery_notes'] = deployments.recovery_date.str.extract(
+                '([A-Za-z\//:]+)', expand=False)
+            deployments.recovery_date = deployments.recovery_date.str.extract(
+                '(\d+-\d+-\d+)', expand=False)
+            deployments = deployments.replace('-', np.nan)
         deployments.loc[deployments.recovery_date.isnull(
         ), 'recovery_date'] = deployments.last_download
 
@@ -62,6 +63,9 @@ def REI(detections, deployments):
         array_unique_species = len(detections.scientificname.unique())
         days_with_detections = len(pd.to_datetime(
             detections.datecollected).dt.date.unique())
+        array_days_active = (max(deployments.last_download.max(
+        ), deployments.recovery_date.max()) - min(deployments.deploy_date)).days
+
         station_reis = pd.DataFrame(columns=['station', 'rei'])
 
         # Loop through each station in the detections and Calculate REI for
@@ -80,11 +84,10 @@ def REI(detections, deployments):
             if name in days_active.index:
                 receiver_days_active = days_active.loc[name].days_deployed.days
                 if receiver_days_active > 0:
-                    rei = ((receiver_unique_tags / array_unique_tags) *
-                           (receiver_unique_species / array_unique_species)
-                           ) / \
-                        (days_with_detections /
-                         receiver_days_with_detections) / receiver_days_active
+                    rei = (receiver_unique_tags / array_unique_tags) * \
+                        (receiver_unique_species / array_unique_species) * \
+                        (receiver_days_with_detections / days_with_detections) * \
+                        (array_days_active / receiver_days_active)
                     station_reis = station_reis.append({
                         'station': name,
                         'rei': rei,

tests/interval_test.py

 # -*- coding: utf-8 -*-
-from resonate.filters import get_distance_matrix
-from resonate.compress import compress_detections
-from resonate.interval_data_tool import interval_data
-import pandas as pd
-import geopy
 import unittest
+
+import geopy
+import pandas as pd
 import pandas.testing as pt
 from colorama import Fore as c
-
+from resonate.compress import compress_detections
+from resonate.filters import get_distance_matrix
+from resonate.interval_data_tool import interval_data
 
 
 class IntervalTest(unittest.TestCase):
 
     def test_filter(self):
-        print( c.YELLOW+'Testing Interval...'+c.RESET)
+        print(c.YELLOW + 'Testing Interval...' + c.RESET)
         input_file = pd.read_csv('tests/assertion_files/nsbs.csv')
-        compressed = compress_detections(input_file) # compressed detections
-        matrix = get_distance_matrix(input_file) # station distance matrix
+        compressed = compress_detections(input_file)  # compressed detections
+        matrix = get_distance_matrix(input_file)  # station distance matrix
 
-        detection_radius = 400 # (in meters) applies same detection radius to all stations
-        station_det_radius = pd.DataFrame([(x, geopy.distance.Distance(detection_radius/1000.0)) for x in matrix.columns.tolist()], columns=['station','radius'])
+        # (in meters) applies same detection radius to all stations
+        detection_radius = 400
+        station_det_radius = pd.DataFrame([(x, geopy.distance.Distance(
+            detection_radius / 1000.0)) for x in matrix.columns.tolist()], columns=['station', 'radius'])
         station_det_radius.set_index('station', inplace=True)
-        station_det_radius # preview radius values
-        dfa = interval_data(compressed_df=compressed, dist_matrix_df=matrix, station_radius_df=station_det_radius)
-
+        station_det_radius  # preview radius values
+        dfa = interval_data(compressed_df=compressed,
+                            dist_matrix_df=matrix, station_radius_df=station_det_radius)
 
         dfb = pd.read_csv('tests/assertion_files/nsbs_interval.csv')
 
@@ -35,12 +37,11 @@ class IntervalTest(unittest.TestCase):
         dfb.to_arrive = pd.to_datetime(dfb.to_arrive)
         dfb.to_leave = pd.to_datetime(dfb.to_leave)
 
-
         dfa.intervaltime = pd.to_timedelta(dfa.intervaltime)
         dfb.intervaltime = pd.to_timedelta(dfb.intervaltime)
 
         pt.assert_frame_equal(dfa, dfb)
-        print( c.GREEN+'OK!\n'+c.RESET)
+        print(c.GREEN + 'OK!\n' + c.RESET)
 
 
 if __name__ == '__main__':