Working prototype for the implementation of reference path dynamic time-warping

measprocess/__init__.py

@@ -2,6 +2,11 @@ from .plotting import plot_locations_osm, plot_street_nodes
 
 from .simulation import TCP_on_lines, ShadowFading
 
+from .dtw import (
+    warp_geodataframe,
+    raw_rp_dtw
+)
+
 from .rtr.process import (
     process,
     interpolate_signal_loc

measprocess/dtw.py

+import geopandas as gpd
+from typing import List, Union
+import pandas as pd
+import numpy as np
+from shapely.geometry import LineString, Point
+from tqdm import tqdm
+
+def raw_rp_dtw(rp, y, dist_func=None, return_path_cost_matrix=False, return_moves_matrices=False, verbose=False, band_width=np.inf, export_path_cost_matrix=None, tqdm_text="Warping..."):
+    '''
+    The raw reference path dtw implementation from Vaclav
+    '''
+
+    if dist_func is None:
+        dist_func = lambda a,b: np.sqrt( np.sum((a-b)**2) )  # Euclidean distance
+
+    if rp.ndim == 1:
+        rp = rp.reshape(-1,1)
+        y = y.reshape(-1,1)
+
+    N = rp.shape[0]
+    M = y.shape[0]
+
+    path_cost_matrix = np.full((N,M), np.inf)
+    best_moves_x = np.full((N,M), np.nan, np.int8)
+    best_moves_y = np.full((N,M), np.nan, np.int8)
+
+    path_cost_matrix[0,0] = dist_func(rp[0,:],y[0,:])
+
+    # Calculate the distance
+    #if verbose: print('Warping...')
+    for ii in tqdm(range(N), desc=tqdm_text, disable=(not verbose)):
+        #if verbose & (np.mod(ii,1000)==0): print('{}/{}'.format(ii,N-1))
+        #for jj in range(M):
+        #for jj in np.arange(start=np.min([np.max([0,ii-band_width]),M-band_width]), stop=np.min([ii+1,M]), step=1):
+        jj_min = np.max([0, ii*(M-1)/(N-1)-np.floor(band_width/2)])
+        jj_max = np.min([M, ii*(M-1)/(N-1)+np.ceil(band_width/2)])
+        for jj in np.arange(start=jj_min, stop=jj_max, step=1, dtype=int):
+            if (ii == 0) and (jj == 0): continue
+            if jj > ii: continue
+
+            # Calculate cost of each allowed action
+            cost_move_right = np.inf
+            cost_move_up    = np.inf
+            cost_move_diag  = np.inf
+            if (ii > 0):
+                cost_move_right = path_cost_matrix[ii-1, jj  ]  # Take next sample from the first signal, keep the second fixed
+
+            # This move is never allowed (we always take next sample from the first signal => the first signal stays unchanged)
+            #if (jj > 0):
+            #    cost_move_up    = path_cost_matrix[ii  , jj-1]  # Take next sample from the second signal, keep the first fixed
+
+            if (ii > 0) and (jj > 0):
+                cost_move_diag  = path_cost_matrix[ii-1, jj-1]  # Take next sample from both signals
+
+            cost_moves = np.array([cost_move_right, cost_move_up, cost_move_diag])
+
+            # Select the best cost (i.e. the best previous position)
+            action_idx = cost_moves.argmin()
+            best_move_cost = cost_moves[action_idx]
+            if action_idx == 0:
+                    best_moves_x[ii,jj] = -1  # We came here from left
+                    best_moves_y[ii,jj] = 0
+            elif action_idx == 1:             # Never happens.
+                    best_moves_x[ii,jj] = 0   # We came here from below
+                    best_moves_y[ii,jj] = -1
+            elif action_idx == 2:
+                    best_moves_x[ii,jj] = -1  # We came here diagonally
+                    best_moves_y[ii,jj] = -1
+
+            path_cost_matrix[ii,jj] = dist_func(rp[ii,:], y[jj,:]) + best_move_cost
+
+    # Calculate the indexes: Go back through the distance matrix, following the best (cheapest) moves
+    #x_idx = np.full(M+N, np.nan)  # Preallocation
+    y_idx = np.full(M+N, np.nan)
+
+    xidx_temp = N-1
+    yidx_temp = M-1
+    ii = 0
+
+    while (xidx_temp > 0) or (yidx_temp > 0):
+        if ii > 0:
+            # Get the previous indexes on the path by a lookup from the "moves" matrixes
+            xidx_prev = xidx_temp + best_moves_x[xidx_temp, yidx_temp]
+            yidx_prev = yidx_temp + best_moves_y[xidx_temp, yidx_temp]
+            xidx_temp = xidx_prev
+            yidx_temp = yidx_prev
+
+        # Store the current x and y indexes
+        #x_idx[ii] = xidx_temp
+        y_idx[ii] = yidx_temp
+
+        ii += 1
+
+    #ix = np.flipud( x_idx[~np.isnan(x_idx)] )
+    iy = np.flipud( y_idx[~np.isnan(y_idx)] ).astype(np.int)
+
+    if export_path_cost_matrix is not None:
+        from scipy.io import savemat
+        savemat(export_path_cost_matrix, {'M':path_cost_matrix, 'idx':iy})
+
+    if return_path_cost_matrix and return_moves_matrices:
+        return iy, path_cost_matrix, best_moves_x, best_moves_y
+    if return_path_cost_matrix and (not return_moves_matrices):
+        return iy, path_cost_matrix
+    if (not return_path_cost_matrix) and return_moves_matrices:
+        return iy, best_moves_x, best_moves_y
+    if (not return_path_cost_matrix) and (not return_moves_matrices):
+        return iy
+
+def warp_geodataframe(
+        reference_path : gpd.GeoSeries,
+        measurement_series : List[gpd.GeoDataFrame],
+        resampling_distance=1, used_crs="EPSG:31287",
+        band_width=np.inf, verbose=True
+    ) -> List[gpd.GeoDataFrame]:
+    '''
+    Reindex a list of measurements to a reference path by applying dynamic timewarping
+
+    :param reference_path: geopandas series with EPSG:4326 crs either consisting of multiple points or a single linestring
+    :param measurement_series: geopandas dataframe with geometry specified by a sequence of points in EPSG:4326
+    :param resampling_distance: distance in meters used for resampling of the reference_path. If None then no resampling will be conducted
+    :param used_crs: the coordinate system used for the distance calculation
+    :param band_width: band_width parameter for the the dynamic time_warping algorithm
+    :param verbose: print a progress animation during warping
+
+    :return:
+        A list of geodadaframes reindexed at the reference path positions
+        the previous measurement sample index is stored under original_index
+    '''
+
+    #Implement the inclusion of a single linestring resampling needs to be set
+    if any(series.crs != "EPSG:4326" for series in [reference_path] + measurement_series):
+        raise ValueError("Pass GeoSeries in EPSG:4326!")
+
+    #Check if reference path is linestring
+    linestring_input = False
+    if any(type(el) == LineString for el in reference_path):
+        linestring_input = True
+        if not len(reference_path) == 1:
+            raise ValueError("Only a single LineString can be used for the reference path")
+        elif resampling_distance is None:
+            raise ValueError("Resampling distance cannot be None for LineString Reference Path")
+
+    #Convert to correct CRS
+    reference_path = reference_path.to_crs(used_crs)
+    measurement_series = [series.to_crs(used_crs) for series in measurement_series]
+
+    #Resample the reference path if a resampling distance is passed
+    if resampling_distance is not None:
+        if linestring_input:
+            reference_path_ls = reference_path.iloc[0]
+        else:
+            reference_path_ls = LineString(reference_path)
+
+        reference_path = gpd.GeoSeries(
+            reference_path_ls.interpolate(dist) for dist in np.arange(0, reference_path_ls.length, resampling_distance)
+        )
+
+    #Basically we should simply summarize the indices here. then append them again.
+    try:
+        index_array = [
+            raw_rp_dtw(
+                np.column_stack((reference_path.x, reference_path.y)),
+                np.column_stack((series.geometry.x, series.geometry.y)),
+                verbose=verbose,
+                tqdm_text=f"Warping {i+1}/{len(measurement_series)}",
+                band_width=band_width
+            ) for i, series in enumerate(measurement_series)
+        ]
+    except IndexError:
+        raise IndexError("Encountered Index Error - Choose smaller resampling distance or remove low speed measurement sequences")
+
+    warped_measurement_series = []
+    for i, indices in enumerate(index_array):
+        #Reindex the warped series
+
+        #Save the previous index under original index
+        old_index_name = measurement_series[i].index.name
+        warped_series = measurement_series[i].loc[indices].reset_index()
+        warped_series = warped_series.rename(columns={old_index_name: 'original_index'})
+
+        #Assign the reference path to the warped series
+        warped_series.geometry = reference_path
+        warped_series = warped_series.to_crs("EPSG:4326")
+
+        warped_measurement_series.append(warped_series)
+
+    return warped_measurement_series
+
+    '''
+    # --> What will we return?
+    # --> we should remove values where we have low velocity.
+    # --> what does the bandwidth parameter actually mean?
+    # --> Out of bounds error might mean that the distance was not correct? resampling should be reset basically
+    # --> Should be run on each of the subexamples --> !!maybe get vaclavs simple example from somewhere!!
+    # --> Would it make sense to pass the reference path as a linestring directly?
+    '''

setup.py

@@ -12,7 +12,7 @@ README = (HERE / "README.md").read_text()
 # This call to setup() does all the work
 setup(
     name="measprocess",
-    version="0.5.15",
+    version="0.6",
     description="Collection of measurement processing tools",
     long_description=README,
     long_description_content_type="text/markdown",

tests/dynamic_time_warping_test.py

+import unittest
+import geopandas as gpd
+from context import measprocess as mpc
+import pandas as pd
+import numpy as np
+
+from shapely.geometry import LineString, Point
+
+class TestDTW(unittest.TestCase):
+    def setUp(self):
+        '''
+        Note, that the Dummy Example does not use coordinates but x, y positions
+        Still we transform it here to be consistent with the default use-case
+        '''
+
+        #Prepare Reference Path
+        rp = pd.read_csv("tests/example_files/dtw_example/reference_path.csv", index_col=0)
+        rp = gpd.GeoSeries(
+            gpd.points_from_xy(rp.x_m, rp.y_m)
+        ).set_crs("EPSG:31287").to_crs("EPSG:4326")
+        self._rp = rp.drop(columns=["x_m, y_m"])
+
+        #Prepare Measurements
+        measurement_files = (
+            pd.read_csv(f"tests/example_files/dtw_example/measurement_{i}.csv", index_col=0)
+            for i in range(3)
+        )
+        self._measurements = [
+            gpd.GeoDataFrame(
+                meas,
+                geometry=gpd.points_from_xy(meas.x_m, meas.y_m)
+            ).set_crs("EPSG:31287").to_crs("EPSG:4326").drop(columns=["x_m", "y_m"]) for meas in measurement_files
+        ]
+
+    def test_epsg_exception(self):
+        rp = self._rp
+        rp.crs = None
+
+        measurements = self._measurements
+        for meas in measurements:
+            meas.crs = None
+
+        with self.assertRaises(ValueError) as context:
+            mpc.dtw.warp_geodataframe(
+                rp, measurements
+            )
+
+    def test_original_index(self):
+        #Check the no measurement was lost
+        measurement = self._measurements[0]
+        warped_measurement = mpc.dtw.warp_geodataframe(
+            self._rp, [measurement],
+            resampling_distance=1
+        )[0]
+
+        self.assertTrue(
+            all(i in warped_measurement.original_index for i in measurement.index)
+        )
+
+    def test_basic_example_resampling(self):
+        measurement_list = mpc.dtw.warp_geodataframe(
+            self._rp, self._measurements,
+            resampling_distance=1
+        )
+
+        self.assertTrue(
+            ((measurement_list[0].geometry == measurement_list[1].geometry) & (measurement_list[1].geometry == measurement_list[2].geometry)).all()
+        )
+
+    def test_basic_example(self):
+        measurement_list = mpc.dtw.warp_geodataframe(
+            self._rp, self._measurements,
+            resampling_distance=None
+        )
+
+        self.assertTrue(
+            ((measurement_list[0].geometry == measurement_list[1].geometry) & (measurement_list[1].geometry == measurement_list[2].geometry)).all()
+        )
+
+    def test_linestring_input(self):
+        rp = gpd.GeoSeries(LineString(self._rp)).set_crs("EPSG:4326")
+
+        measurement_list = mpc.dtw.warp_geodataframe(
+            rp, self._measurements,
+            resampling_distance=1
+        )
+
+        self.assertTrue(
+            ((measurement_list[0].geometry == measurement_list[1].geometry) & (measurement_list[1].geometry == measurement_list[2].geometry)).all()
+        )
+
+    def test_linestring_input_warning(self):
+        rp = gpd.GeoSeries(LineString(self._rp)).set_crs("EPSG:4326")
+
+        measurement_list = mpc.dtw.warp_geodataframe(
+            rp, self._measurements,
+            resampling_distance=1
+        )
+
+        self.assertTrue(
+            ((measurement_list[0].geometry == measurement_list[1].geometry) & (measurement_list[1].geometry == measurement_list[2].geometry)).all()
+        )
+
+    def test_index_warning(self):
+        rp = gpd.GeoSeries(LineString(self._rp)).set_crs("EPSG:4326")
+
+        with self.assertRaises(IndexError) as context:
+            measurement_list = mpc.dtw.warp_geodataframe(
+                rp, self._measurements,
+                resampling_distance=5
+            )

tests/example_files/dtw_example/reference_path.csv

+,x_m,y_m
+0,10,89
+1,11,89
+2,12,89
+3,13,89
+4,14,89
+5,15,89
+6,16,89
+7,17,89
+8,18,89
+9,19,89
+10,20,89
+11,21,89
+12,22,89
+13,23,89
+14,24,89
+15,25,89
+16,26,89
+17,26,88
+18,26,87
+19,26,86
+20,26,85
+21,26,84
+22,26,83
+23,26,82
+24,26,81
+25,26,80
+26,26,79
+27,26,78
+28,26,77
+29,26,76
+30,26,75
+31,26,74
+32,26,73
+33,26,72
+34,26,71
+35,26,70
+36,26,69
+37,26,68
+38,26,67
+39,26,66
+40,26,65
+41,26,64
+42,26,63
+43,26,62
+44,26,61
+45,27,61
+46,28,61
+47,29,61
+48,30,61
+49,31,61
+50,32,61
+51,33,61
+52,34,61
+53,35,61
+54,36,61
+55,37,61
+56,38,61
+57,39,61
+58,40,61
+59,41,61
+60,42,61
+61,43,61
+62,44,61
+63,45,61
+64,46,61
+65,47,61
+66,48,61
+67,49,61
+68,50,61
+69,51,61
+70,52,61
+71,53,61
+72,54,61
+73,55,61
+74,56,61
+75,57,61
+76,58,61
+77,59,61
+78,60,61
+79,61,61
+80,62,61
+81,63,61
+82,64,61
+83,65,61
+84,66,61
+85,67,61
+86,68,61
+87,69,61
+88,70,61
+89,71,61
+90,72,61
+91,73,61
+92,74,61
+93,75,61
+94,76,61
+95,77,61
+96,78,61
+97,78,62
+98,78,63
+99,78,64
+100,78,65
+101,78,66
+102,78,67
+103,78,68
+104,78,69
+105,78,70
+106,78,71
+107,78,72
+108,78,73
+109,78,74
+110,78,75
+111,78,76
+112,78,77
+113,78,78
+114,78,79
+115,78,80
+116,78,81
+117,78,82
+118,78,83
+119,78,84
+120,78,85
+121,78,86
+122,78,87
+123,78,88
+124,78,89
+125,78,90
+126,78,91
+127,78,92
+128,78,93
+129,78,94
+130,78,95
+131,78,96
+132,78,97
+133,78,98
+134,78,99
+135,78,100
+136,78,101
+137,78,102
+138,78,103
+139,78,104
+140,78,105
+141,78,106
+142,78,107
+143,78,108
+144,78,109
+145,78,110
+146,78,111
+147,78,112
+148,78,113
+149,78,114
+150,78,115
+151,78,116
+152,78,117
+153,78,118
+154,78,119
+155,78,120
+156,78,121
+157,78,122
+158,78,123
+159,78,124
+160,78,125
+161,78,126
+162,79,126
+163,80,126
+164,81,126
+165,82,126
+166,83,126
+167,84,126
+168,85,126
+169,86,126
+170,87,126
+171,88,126
+172,89,126
+173,90,126
+174,91,126
+175,92,126
+176,93,126
+177,94,126
+178,95,126
+179,96,126
+180,97,126
+181,98,126
+182,99,126
+183,100,126
+184,101,126
+185,102,126
+186,103,126
+187,104,126
+188,105,126
+189,106,126
+190,107,126
+191,108,126
+192,109,126
+193,110,126
+194,111,126
+195,112,126
+196,113,126
+197,114,126
+198,115,126
+199,116,126
+200,117,126
+201,118,126
+202,119,126
+203,120,126
+204,121,126
+205,122,126
+206,123,126
+207,124,126
+208,125,126
+209,126,126
+210,127,126
+211,128,126
+212,129,126
+213,130,126
+214,131,126
+215,132,126
+216,133,126
+217,134,126
+218,135,126
+219,136,126
+220,137,126
+221,138,126
+222,139,126
+223,139,125
+224,139,124
+225,139,123
+226,139,122
+227,139,121
+228,139,120
+229,139,119
+230,139,118
+231,139,117
+232,139,116
+233,139,115
+234,139,114
+235,139,113
+236,139,112
+237,139,111
+238,139,110
+239,139,109
+240,139,108
+241,139,107
+242,139,106
+243,139,105
+244,139,104
+245,139,103
+246,139,102
+247,139,101
+248,139,100
+249,139,99
+250,139,98
+251,139,97
+252,139,96
+253,139,95
+254,139,94
+255,139,93
+256,139,92
+257,139,91
+258,139,90
+259,139,89
+260,139,88
+261,139,87
+262,139,86
+263,139,85
+264,139,84
+265,139,83
+266,139,82
+267,139,81
+268,139,80
+269,139,79
+270,139,78
+271,139,77
+272,139,76
+273,139,75
+274,139,74
+275,139,73
+276,139,72
+277,139,71
+278,139,70
+279,139,69
+280,139,68
+281,139,67
+282,139,66
+283,139,65
+284,139,64
+285,139,63
+286,139,62
+287,139,61
+288,139,60
+289,139,59
+290,139,58
+291,139,57
+292,139,56
+293,139,55
+294,139,54
+295,139,53
+296,139,52
+297,139,51
+298,139,50
+299,139,49
+300,139,48
+301,139,47
+302,139,46
+303,139,45
+304,139,44
+305,138,44
+306,137,44
+307,136,44
+308,135,44
+309,134,44
+310,133,44
+311,132,44
+312,131,44
+313,130,44
+314,129,44
+315,128,44
+316,127,44
+317,126,44
+318,125,44
+319,124,44
+320,123,44
+321,122,44
+322,121,44
+323,120,44
+324,119,44
+325,118,44
+326,117,44
+327,116,44
+328,115,44
+329,114,44
+330,113,44
+331,112,44
+332,111,44
+333,110,44
+334,109,44
+335,108,44
+336,107,44
+337,106,44
+338,105,44
+339,104,44
+340,103,44
+341,102,44
+342,101,44
+343,101,45
+344,101,46
+345,101,47
+346,101,48
+347,101,49
+348,101,50
+349,101,51
+350,101,52
+351,101,53
+352,101,54
+353,101,55
+354,101,56
+355,101,57
+356,101,58
+357,101,59