Add graph-based linejoin

fronts.py

@@ -3,6 +3,7 @@
 import scipy.signal
 import xarray as xr
 import scipy.spatial.distance as sp_dist
+import scipy.spatial
 import geopy.distance as gp_dist
 
 def wetbulb(ta,hus,plev,steps=100,ta_units=None):
@@ -199,6 +200,64 @@ def linejoin(inpts,searchdist=1.5,minlength=250,lonex=0):
 			filt_lines.append(line)
 	lines=filt_lines
 	return lines
+
+
+def linejoin_graph(inpts, searchdist: float=1.5, minlength: float=250, lonex: float=0):
+    """
+    Turns a list of lat-lon points into a list of joined lines
+
+    Args:
+        inpts - the list of points (list of lat-lon points)
+        searchdist - degree radius around each point that other points within are
+                     deemed to be part of the same line
+        minlength - minimum end-to-end length of the lines (in km)
+        lonex - minimum end-to-end longitudinal extent
+
+    Returns:
+        List with each member a tuple of (lats, lons) for each line found
+    """
+
+    tree = scipy.spatial.KDTree(inpts)
+
+    assert tree.m == 2, f"Expected 2d input points, found {tree.m}d"
+
+    distances = tree.sparse_distance_matrix(tree, searchdist)
+
+    # Cull to only the minimal connections
+    span = scipy.sparse.csgraph.minimum_spanning_tree(distances)
+
+    # Get the components
+    ncomp, comp_labels = scipy.sparse.csgraph.connected_components(span)
+
+    lines = []
+
+    for c in range(ncomp):
+        comp_indices = np.arange(comp_labels.size)[comp_labels == c]
+        if len(comp_indices) <= 1:
+            continue
+
+        start = comp_indices[1]
+
+        # This ordering may not start at the start of the span, but it will end at the end of it
+        order = scipy.sparse.csgraph.depth_first_order(span, start, return_predecessors=False, directed=False)
+
+        # Search the other way to get the correct ordering for the full path
+        order = scipy.sparse.csgraph.depth_first_order(span, order[-1], return_predecessors=False, directed=False)
+
+        # Add the coordinates from this path to the output list
+        lines.append((inpts[order,0], inpts[order,1]))
+
+    # Filter from the previous version
+    filt_lines=[]
+    for line in lines:
+        ln_dist=gp_dist.distance((line[0][0],line[1][0]),(line[0][-1],line[1][-1])).km
+        lon_extent=max(line[1])-min(line[1])
+        if ln_dist>minlength and lon_extent>lonex:
+            filt_lines.append(line)
+    lines=filt_lines
+
+    return lines
+
 
 def smoother(data,numsmooth=9,smooth_kernel=np.ones((3,3))/9):
 	# smooths an input 2-d xarray using a given kernel and number of passes