Create function `node_to_prop`

We extract part of function `calculate_radii_rossby` into new function
`node_to_prop`. So we now separate reading data from the shapefiles
and computing the average radius and Rossby number. Thus, we prepare
for rewriting the algorithm in the main part of the script. The
quantity of reading is still the same, for now.

cost_function.py

@@ -28,9 +28,11 @@ import argparse
 
 Omega = 2 * math.pi / 86164.
 
-def calculate_radii_rossby(list_eddies, e_overestim, handlers, array_d_init):
-    """Compute average on list_eddies of Rossby number and radius of
-    maximum speed contour.
+def calculate_radii_rossby(properties):
+    """Compute average on some instantaneous eddies of Rossby number and
+    radius of maximum speed contour. The required properties for each
+    eddy are position, radius and speed. "properties" is a list of
+    dictionaries. Each dictionary in the list contains the three properties.
 
     """
 
@@ -38,26 +40,17 @@ def calculate_radii_rossby(list_eddies, e_overestim, handlers, array_d_init):
     avg_Rossby = 0
     n_valid_Rossby = 0
 
-    for n in list_eddies:
-        date_index, eddy_index = report_graph.node_to_date_eddy(n, e_overestim)
-        i_SHPC = get_SHPC(array_d_init, date_index)
-        ishape = util_eddies.comp_ishape(handlers[i_SHPC], date_index,
-                                         eddy_index)
+    for prop in properties:
+        f = 2 * Omega * math.sin(math.radians(prop["pos"][1])) # in s-1
+        radius = prop["radius"] * 1000 # in m
 
-        # Now that we have the location in the shapefiles, we need to
-        # get the radius and the rossby number:
-        shapeRec = handlers[i_SHPC]["readers"]["extremum"].shapeRecord(ishape)
-        f = 2 * Omega * math.sin(math.radians(shapeRec.shape.points[0][1]))
-        radius = handlers[i_SHPC]["readers"]["max_speed_contour"]\
-            .record(ishape).r_eq_area * 1000 # in m
-
-        if abs(shapeRec.record.speed) < 100:
-            avg_Rossby += shapeRec.record.speed / (f * radius)
+        if abs(prop["speed"]) < 100:
+            avg_Rossby += prop["speed"] / (f * radius)
             n_valid_Rossby += 1
 
         avg_rad += radius # in m
 
-    avg_rad /= len(list_eddies)
+    avg_rad /= len(properties)
     if n_valid_Rossby != 0: avg_Rossby /= n_valid_Rossby
     return avg_rad, avg_Rossby
 
@@ -66,6 +59,30 @@ def get_SHPC(array_d_ini, date_index):
     assert i_SHPC >= 0
     return i_SHPC
 
+def node_to_prop(node_list, e_overestim, array_d_init, handlers):
+    """node_list is a list of node identification numbers for
+    instantaneous eddies. This function returns some properties of the
+    eddies, read from shapefiles: position of extremum, radius of
+    outermost contour or maximum speed contour, and speed. The three
+    properties are in a dictionary, for each eddy.
+
+    """
+
+    properties = []
+
+    for n in node_list:
+        date_index, eddy_index = report_graph.node_to_date_eddy(n, e_overestim)
+        i_SHPC = get_SHPC(array_d_init, date_index)
+        ishape = util_eddies.comp_ishape(handlers[i_SHPC], date_index,
+                                         eddy_index)
+        shapeRec = handlers[i_SHPC]["readers"]["extremum"].shapeRecord(ishape)
+        prop = {"pos": shapeRec.shape.points[0], "speed": shapeRec.record.speed}
+        prop["radius"] = handlers[i_SHPC]["readers"]["max_speed_contour"]\
+            .record(ishape).r_eq_area
+        properties.append(prop)
+
+    return properties
+
 t0 = time.perf_counter()
 timings = open("timings.txt", "w")
 parser = argparse.ArgumentParser()
@@ -156,20 +173,23 @@ for n in g.vertices():
             # to average
 
             # First 7 days calculation
+            properties = node_to_prop(segment[:n_days_avg], e_overestim,
+                                      array_d_init, handlers)
             g.vp.first_av_rad[n], g.vp.first_av_ros[n] \
-                = calculate_radii_rossby(segment[:n_days_avg], e_overestim,
-                                         handlers, array_d_init)
+                = calculate_radii_rossby(properties)
 
             # Last 7 days calculation:
+            properties = node_to_prop(segment[- n_days_avg:], e_overestim,
+                                      array_d_init, handlers)
             g.vp.last_av_rad[n], g.vp.last_av_ros[n] \
-                = calculate_radii_rossby(segment[- n_days_avg:], e_overestim,
-                                         handlers, array_d_init)
+                = calculate_radii_rossby(properties)
         else:
             # The number of eddies in a segment is lower than the number
             # of days over which to average. The values will be the same
             # except for the positions.
-            avg_rad, avg_Rossby = calculate_radii_rossby(segment, e_overestim,
-                                                         handlers, array_d_init)
+            properties = node_to_prop(segment, e_overestim, array_d_init,
+                                      handlers)
+            avg_rad, avg_Rossby = calculate_radii_rossby(properties)
 
             # Average and assign the rossbies:
             g.vp.first_av_ros[n] = avg_Rossby