added tracking of aews with filtering based on avg speed over time frame

enstools/feature/identification/african_easterly_waves/configuration.py

@@ -39,7 +39,7 @@ levels = [70000]  # 700 hPa
 # time of interest, if None all
 # june-oct is AEW season
 start_date = '2008-08-01T00:00' #  # '2008-08-01T00:00'
-end_date = '2008-08-15T00:00' # None # '2008-08-03T00:00'
+end_date = '2008-08-24T00:00' # None # '2008-08-03T00:00'
 
 # Algorithm parameters
 # max u wind (m/s) (0 = only keep west-propagating). Belanger: 2.5; Berry: 0.0
@@ -63,6 +63,11 @@ duration_threshold = timedelta(days=2)
 # max time allowed between detected waves
 max_cmp_delta = timedelta(hours=12)
 
+# duration of min avg speed
+avg_speed_t = timedelta(hours=48)
+avg_speed_min_m_per_s = 3
+avg_speed_min_deg_per_h = avg_speed_min_m_per_s * 3.6 * 0.00914
+
 # speed range of AEWs
 # at 10°N we have in longitude direction 0.00914 degrees/km (360/(40,000*cos(10deg)))
 speed_range_m_per_s = [3.0, 12.0] # [5,10], but be more gentle with polygons.

enstools/feature/tracking/african_easterly_waves/tracking.py

@@ -2,10 +2,14 @@ from ..object_compare_tracking import ObjectComparisonTracking
 
 import enstools.feature.identification.african_easterly_waves.configuration as cfg
 from datetime import datetime, timedelta
+from collections import defaultdict
+import statistics
 
 from shapely.geometry import Polygon, LineString
 from enstools.feature.util.data_utils import pb_str_to_datetime
+from enstools.feature.util.graph import DataGraph
 
+zero_dt = timedelta(seconds=0)
 
 class AEWTracking(ObjectComparisonTracking):
     """
@@ -49,7 +53,7 @@ class AEWTracking(ObjectComparisonTracking):
             p2_pt_list.append((lineseg2[i].lon, lineseg2[i].lat))
         actual_linestring = LineString(p2_pt_list)
 
-        if not predicted_polygon.intersects(actual_linestring): # linesting poly not where expected?
+        if not predicted_polygon.intersects(actual_linestring):  # linesting poly not where expected?
             return False
 
         # if bb1_lon_center - bb2_lon_center  > time_diff_h / 2: # TODO weird
@@ -61,39 +65,37 @@ class AEWTracking(ObjectComparisonTracking):
         print("postprocess @ aew")
         return
 
+    # get dict of key=endnode, values=nodes on way) for endnodes at end if time_delta
+    @staticmethod
+    def get_nodes_after(time_delta, edges, start_index):
 
+        this_node = edges[start_index].parent
+        if time_delta <= zero_dt or len(edges[start_index].childs) == 0:
+            dd = defaultdict(list) # TODO set instead of list!
+            dd[start_index] = [start_index]
+            return dd
 
-    def get_downstream_nodes_after(self, time_delta, edges, start_index):
+        obj_node_list = [e.parent for e in edges]
+        childs_ = defaultdict(list)
 
-        if time_delta < 0:
-            ??
+        for child in edges[start_index].childs:
 
-        node_indices = [start_index]
-        obj_node_list = [e.parent for e in edges]  # this up
-
-        co = edges[start_index]
-        node, connected_nodes = co.parent, co.childs
-
-
-
-        for c_node in connected_nodes:
-            dt = pb_str_to_datetime(c_node.time) - pb_str_to_datetime(node.time)
+            dt = pb_str_to_datetime(child.time) - pb_str_to_datetime(this_node.time)
             remaining_time_delta = time_delta - dt
 
             # get index of connected node in graph
-            c_node_idx = obj_node_list.index(c_node)
-
-            # call recursively on this connected node
-            end_way_node_list = self.get_downstream_nodes_after(remaining_time_delta, edges, c_node_idx)
-            node_indices.extend(c_node_downstream_indices)
-
-
-        return [final node, nodes on the way]
+            c_node_idx = obj_node_list.index(child)
+            after_nodes = AEWTracking.get_nodes_after(remaining_time_delta, edges, c_node_idx)
+            # add own node to it (this_node)
+            for endn_, pathn_ in after_nodes.items():
+                childs_[endn_].extend(pathn_ + [start_index])
 
+        return childs_
 
+    def adjust_track(self, tracking_graph):
 
-    def adjust_track(self, track): # TODO track should be graph?
         # keep track if persists longer than duration_threshold
+        track = tracking_graph.graph
 
         nodes = [edge.parent for edge in track.edges]
         min_time = pb_str_to_datetime(nodes[0].time)
@@ -106,59 +108,75 @@ class AEWTracking(ObjectComparisonTracking):
         times = [pb_str_to_datetime(node.time) for node in nodes]
         center_lons = [(node.object.properties.bb.max.lon + node.object.properties.bb.min.lon) / 2.0 for node in nodes]
         center_lats = [(node.object.properties.bb.max.lat + node.object.properties.bb.min.lat) / 2.0 for node in nodes]
+        extent_lons = [(node.object.properties.bb.max.lon - node.object.properties.bb.min.lon) for node in nodes]
+        extent_lats = [(node.object.properties.bb.max.lat - node.object.properties.bb.min.lat) for node in nodes]
 
         # TODO need general statement that nodes are part of exactly one track.
 
         # init nodes as False = should be discarded
         keep_node = [False] * len(nodes)
 
+        # TODO throw out also nodes which are wrongly orientated on their own?! no, once is fine...
+
         # iterate over all nodes: if fast enough from this to this+dt, keep
         for node_idx, node in enumerate(nodes):
-            node_time_plus_delta = pb_str_to_datetime(node.time) + self.config.delta48
 
-            # TODO v this should instantly return pairs (endnote, all nodes on path) for later
             # downstream_nodes_with_path = self.get_downstream_nodes_after(self.config.slow_duration_window, track.edges, node_idx)
             # TODO if none? maybe skipped?
-            ds_nodes = AEWTracking.get_downstream_node_indices(nodes, node_idx, until_time=node_time_plus_delta)
-
-            end_nodes = ends_of(ds_nodes)
-            for end_node in end_nodes:
-
-                ds_node_lon_center = (ds.object.properties.bb.max.lon + ds.object.properties.bb.min.lon) / 2.0
+            # get list of (end_node, nodes_on_way)
+            ds_nodes = AEWTracking.get_nodes_after(self.config.avg_speed_t, track.edges, node_idx)
 
-                deg_per_hr_avg = (ds_node_lon_center - center_lons[node_idx]) / self.config.slow_duration_window.hours
-                if abs(deg_per_hr_avg) > abs(self.config.avg_speed_thr):
+            for end_node_idx, way_node_idxs in ds_nodes.items():
+                way_node_idxs = list(set(way_node_idxs))
+                end_node = nodes[end_node_idx]
 
-                if speed A->B is not fine:
+                # if dt < 48h -> didnt reach 2 days.
+                if pb_str_to_datetime(end_node.time) - pb_str_to_datetime(node.time) < self.config.avg_speed_t:
+                    # could not track wave that long, dont set to True
                     continue
 
+                # if nodes along track are predominantly horizontal orientated discard them.
+                lon_lat_ratio = [extent_lons[i] / extent_lats[i] for i in way_node_idxs]
+                if statistics.mean(lon_lat_ratio) > 1:
+                    continue # TODO doc this
 
-                us_nodes = get_upstream_nodes_at_time(t)
-                nodes_in_path_ab = cut(ds_nodes, us_nodes)
+                # check avg speed in that 48h
+                end_node_lon_center = (end_node.object.properties.bb.max.lon + end_node.object.properties.bb.min.lon) / 2.0
 
-                keep_node[nodes_in_path_ab] = True
+                avg_speed_window_h = self.config.avg_speed_t.total_seconds() / 3600
+                deg_per_hr_avg = (end_node_lon_center - center_lons[node_idx]) / avg_speed_window_h
+                # speed check. fast enough -> set to True
+                if abs(deg_per_hr_avg) > abs(self.config.avg_speed_min_deg_per_h):
+                    for way_node_idx in way_node_idxs:  # TODO indices of way_nodes
+                        keep_node[way_node_idx] = True
 
-            """
-            # for each +48h node: check if avg speed fulfilled. if yes, mark this path.
-            for ds, path_nodes in downstream_nodes_with_path:
-                ds_node_lon_center = (ds.object.properties.bb.max.lon + ds.object.properties.bb.min.lon) / 2.0
-
-                deg_per_hr_avg = (ds_node_lon_center - center_lons[node_idx]) / self.config.slow_duration_window.hours
-                if abs(deg_per_hr_avg) > abs(self.config.avg_speed_thr):
-                    # fast enough, keep all nodes along this path.
-                    for path_node in path_nodes:
-                        keep_node[path_node] = True
-            """
-
-        kept_nodes = [node for x, node in enumerate(nodes) if keep_node[x]]
+                # us_nodes = get_upstream_nodes_at_time(t)
+                # nodes_in_path_ab = cut(ds_nodes, us_nodes)
 
+        # clean up -> degenerate nodes (child without connection TODO cleanup)
+        # then generate_tracks()
+        kept_edges = [edge for x, edge in enumerate(track.edges) if keep_node[x]]
         # regenerate tracks out of leftover edges. this also splits a track etc, and removes short ones
-        new_tracks = generate_tracks(kept_nodes)
 
-        return new_tracks
+        # create return object
+        reduced_tg = tracking_graph.tr_tech.pb_reference.TrackableSet()
+        reduced_tg.CopyFrom(tracking_graph.set_desc)
+        del reduced_tg.graph.edges[:]
+        del reduced_tg.tracks[:]
+        reduced_tg.graph.edges.extend(kept_edges)
+
+        dg = DataGraph(reduced_tg, self, fix=True)
+        # if dg nodes same as input to this method, nothing changed, dont regenerate tracks.
+        # otherwise stuck in infinite loop
+        if len(dg.graph.edges) == len(tracking_graph.graph.edges):
+            # TODO here graph.tracks stays empty!!!
+            print("Nothing changed.")
+            return [dg]
 
+        print("Regenerate.")
+        tracks = dg.generate_tracks(apply_filter=True)
 
         # TODO
         #   redo clim, plot which WTs not considered
 
-        return track
+        return tracks

enstools/feature/tracking/tracking.py

@@ -120,7 +120,7 @@ class TrackingTechnique(ABC):
 
         pass
 
-    def adjust_track(self, track):
+    def adjust_track(self, track): # TODO is DataGraph, update
         """
         Override this method to filter, change and split certain tracks.
         When generate_tracks() on a Graph is called, each track is checked via this method.

enstools/feature/util/graph.py

@@ -8,7 +8,7 @@ class DataGraph:
     It contains methods to iterate a graph and to generate the tracks.
     """
 
-    def __init__(self, set_desc, tracking_tech):
+    def __init__(self, set_desc, tracking_tech, fix=False):
         self.set_desc = set_desc
 
         # sort the edges in the graph by time.
@@ -18,6 +18,9 @@ class DataGraph:
 
         self.tr_tech = tracking_tech
 
+        if fix:
+            self.fix_graph()
+
     def get_earliest_nodes(self):
         """
         Get the nodes which are earliest in the graph.
@@ -64,6 +67,24 @@ class DataGraph:
 
         return nodes
 
+    def fix_graph(self):
+        """
+        Fixes graph by removing orphaned edges. (Edges where endnode is not a valid node)
+
+        Returns
+        -------
+        """
+        nodes = [e.parent for e in self.graph.edges]
+
+        print("Fixing graph")
+        print("")
+        for edge in self.graph.edges:
+            for child in edge.childs[:]:
+                if child not in nodes:
+                    print("Remove edge")
+                    edge.childs.remove(child)
+
+
     def get_parents(self, node):
         """
         Get the parents of a node (previous timestep)
@@ -161,9 +182,16 @@ class DataGraph:
             track = self.tr_tech.pb_reference.ObjectGraph()
             track.edges.extend(cur_object_nodes)
 
+            # make a graph out of it
+            track_set_ = self.tr_tech.pb_reference.TrackableSet()
+            track_set_.CopyFrom(self.set_desc)
+            track_set_.graph.CopyFrom(track)
+
+            track_graph = DataGraph(track_set_, self.tr_tech)
+
             if apply_filter:
                 # let user change track, reject it, or split it
-                current_tracks = self.tr_tech.adjust_track(track)
+                current_tracks = self.tr_tech.adjust_track(track_graph)
                 # if user just returns the track itself, make it a list
                 if isinstance(current_tracks, self.tr_tech.pb_reference.ObjectGraph):
                     current_tracks = [current_tracks]
@@ -174,9 +202,10 @@ class DataGraph:
                 current_tracks = [track]
 
             # did user reject this track?
-            if current_tracks is None or current_tracks is []:
+            if current_tracks is None or current_tracks == []:
                 continue
 
+            # TODO current_tracks not used????
             # add to desc
             self.set_desc.tracks.append(track)
 

requirements.txt

 # for base functionality
 protobuf
 netcdf4
-
+metpy