Recurse algorithm

ecoscope/analysis/__init__.py

@@ -1,6 +1,7 @@
 from ecoscope.analysis import UD, astronomy, seasons
 from ecoscope.analysis.ecograph import Ecograph, get_feature_gdf
 from ecoscope.analysis.percentile import get_percentile_area
+from ecoscope.analysis.recurse import get_recursions
 from ecoscope.analysis.speed import SpeedDataFrame
 
 __all__ = [
@@ -11,4 +12,5 @@
     "get_feature_gdf",
     "get_percentile_area",
     "seasons",
+    "get_recursions",
 ]

ecoscope/analysis/ecograph.py

@@ -7,7 +7,6 @@
 import pandas as pd
 import rasterio
 import sklearn.base
-from affine import Affine
 from shapely.geometry import shape
 from skimage.draw import line
 
@@ -58,19 +57,7 @@ def __init__(self, trajectory, resolution=15, radius=2, cutoff=None, tortuosity_
         eastings = np.array([geom.iloc[i].coords.xy[0] for i in range(len(geom))]).flatten()
         northings = np.array([geom.iloc[i].coords.xy[1] for i in range(len(geom))]).flatten()
 
-        self.xmin = floor(np.min(eastings)) - self.resolution
-        self.ymin = floor(np.min(northings)) - self.resolution
-        self.xmax = ceil(np.max(eastings)) + self.resolution
-        self.ymax = ceil(np.max(northings)) + self.resolution
-
-        self.xmax += self.resolution - ((self.xmax - self.xmin) % self.resolution)
-        self.ymax += self.resolution - ((self.ymax - self.ymin) % self.resolution)
-
-        self.transform = Affine(self.resolution, 0.00, self.xmin, 0.00, -self.resolution, self.ymax)
-        self.inverse_transform = ~self.transform
-
-        self.n_rows = int((self.xmax - self.xmin) // self.resolution)
-        self.n_cols = int((self.ymax - self.ymin) // self.resolution)
+        self.grid = ecoscope.base.Grid(eastings, northings, self.resolution)
 
         def compute(df):
             subject_name = df.name
@@ -143,7 +130,7 @@ def to_geotiff(self, feature, output_path, individual="all", interpolation=None,
             band_count=1,
         )
         raster_profile.raster_extent = ecoscope.io.raster.RasterExtent(
-            x_min=self.xmin, x_max=self.xmax, y_min=self.ymin, y_max=self.ymax
+            x_min=self.grid.xmin, x_max=self.grid.xmax, y_min=self.grid.ymin, y_max=self.grid.ymax
         )
         ecoscope.io.raster.RasterPy.write(
             ndarray=feature_ndarray,
@@ -160,8 +147,8 @@ def _get_ecograph(self, trajectory_gdf, individual_name, radius, cutoff, tortuos
             lines = [list(geom.iloc[i + j].coords) for j in range(tortuosity_length)]
             p1, p2, p3, p4 = lines[0][0], lines[0][1], lines[1][1], lines[1][0]
             pixel1, pixel2 = (
-                self.inverse_transform * p1,
-                self.inverse_transform * p2,
+                self.grid.inverse_transform * p1,
+                self.grid.inverse_transform * p2,
             )
 
             row1, row2 = floor(pixel1[0]), floor(pixel2[0])
@@ -296,9 +283,9 @@ def _get_feature_mosaic(self, feature, interpolation=None):
         features = []
         for individual in self.graphs.keys():
             features.append(self._get_feature_map(self, feature, individual, interpolation))
-        mosaic = np.full((self.n_cols, self.n_rows), np.nan)
-        for i in range(self.n_cols):
-            for j in range(self.n_rows):
+        mosaic = np.full((self.grid.n_cols, self.grid.n_rows), np.nan)
+        for i in range(self.grid.n_cols):
+            for j in range(self.grid.n_rows):
                 values = []
                 for feature_map in features:
                     grid_val = feature_map[i][j]
@@ -319,7 +306,7 @@ def _get_feature_map(self, feature, individual, interpolation):
 
     @staticmethod
     def _get_regular_feature_map(self, feature, individual):
-        feature_ndarray = np.full((self.n_cols, self.n_rows), np.nan)
+        feature_ndarray = np.full((self.grid.n_cols, self.grid.n_rows), np.nan)
         for node in self.graphs[individual].nodes():
             feature_ndarray[node[1]][node[0]] = (self.graphs[individual]).nodes[node][feature]
         return feature_ndarray
@@ -333,7 +320,7 @@ def _get_interpolated_feature_map(self, feature, individual, interpolation):
         for i in range(len(geom)):
             line1 = list(geom.iloc[i].coords)
             p1, p2 = line1[0], line1[1]
-            pixel1, pixel2 = self.inverse_transform * p1, self.inverse_transform * p2
+            pixel1, pixel2 = self.grid.inverse_transform * p1, self.grid.inverse_transform * p2
             row1, row2 = floor(pixel1[0]), floor(pixel2[0])
             col1, col2 = ceil(pixel1[1]), ceil(pixel2[1])
 

ecoscope/analysis/recurse.py

@@ -0,0 +1,49 @@
+from math import ceil, floor
+
+import numpy as np
+
+import ecoscope
+
+
+def get_consecutive_items_number(idxs):
+    gaps = [[s, e] for s, e in zip(idxs, idxs[1:]) if s + 1 < e]
+    edges = iter(idxs[:1] + sum(gaps, []) + idxs[-1:])
+    return len(list(zip(edges, edges)))
+
+
+def get_recursions(relocations, resolution):
+    relocations = relocations.reset_index(drop=True)
+    if not relocations["fixtime"].is_monotonic:
+        relocations.sort_values("fixtime", inplace=True)
+
+    diameter = ceil(resolution) * 2
+    utm_crs = relocations.estimate_utm_crs()
+    relocations.to_crs(utm_crs, inplace=True)
+
+    geom = relocations["geometry"]
+    eastings = np.array([geom.iloc[i].x for i in range(len(geom))]).flatten()
+    northings = np.array([geom.iloc[i].y for i in range(len(geom))]).flatten()
+
+    grid = ecoscope.base.Grid(eastings, northings, diameter)
+
+    grid_cells_dict = {}
+    for i in range(len(geom)):
+        point = geom.iloc[i]
+        grid_cell = grid.inverse_transform * (point.x, point.y)
+        row, col = floor(grid_cell[0]), ceil(grid_cell[1])
+        if (col, row) in grid_cells_dict:
+            grid_cells_dict[(col, row)].append(i)
+        else:
+            grid_cells_dict[(col, row)] = [i]
+
+    for grid_cell in grid_cells_dict:
+        grid_cells_dict[grid_cell] = get_consecutive_items_number(grid_cells_dict[grid_cell])
+
+    recursion_values = []
+    for i in range(len(geom)):
+        point = geom.iloc[i]
+        grid_cell = grid.inverse_transform * (point.x, point.y)
+        row, col = floor(grid_cell[0]), ceil(grid_cell[1])
+        recursion_values.append(grid_cells_dict[(col, row)])
+
+    return recursion_values

ecoscope/base/__init__.py

@@ -7,6 +7,7 @@
 )
 from ecoscope.base.base import EcoDataFrame, Relocations, Trajectory
 from ecoscope.base.utils import (
+    Grid,
     cachedproperty,
     create_meshgrid,
     groupby_intervals,
@@ -24,4 +25,5 @@
     "cachedproperty",
     "create_meshgrid",
     "groupby_intervals",
+    "Grid",
 ]

ecoscope/base/utils.py

@@ -1,6 +1,9 @@
+from math import ceil, floor
+
 import geopandas as gpd
 import numpy as np
 import pandas as pd
+from affine import Affine
 from shapely.geometry import box
 
 
@@ -212,3 +215,65 @@ def create_modis_interval_index(start, intervals, overlap=pd.Timedelta(0), close
     left = pd.DatetimeIndex(left)
 
     return pd.IntervalIndex.from_arrays(left=left, right=left + pd.Timedelta(days=16), closed=closed)
+
+
+class Grid:
+    """
+    A class that creates a grid covering a list of UTM points
+
+    Parameters
+    ----------
+    eastings : list
+        UTM easting coordinates
+    northings : list
+        UTM northing coordinates
+    resolution : float
+        The width/length of a grid cell in meters
+    """
+
+    def __init__(self, eastings, northings, resolution):
+        self._xmin = floor(np.min(eastings)) - resolution
+        self._ymin = floor(np.min(northings)) - resolution
+        self._xmax = ceil(np.max(eastings)) + resolution
+        self._ymax = ceil(np.max(northings)) + resolution
+
+        self._xmax += resolution - ((self._xmax - self._xmin) % resolution)
+        self._ymax += resolution - ((self._ymax - self._ymin) % resolution)
+
+        self._transform = Affine(resolution, 0.00, self._xmin, 0.00, -resolution, self._ymax)
+        self._inverse_transform = ~self._transform
+
+        self._n_rows = int((self._xmax - self._xmin) // resolution)
+        self._n_cols = int((self._ymax - self._ymin) // resolution)
+
+    @property
+    def xmin(self):
+        return self._xmin
+
+    @property
+    def ymin(self):
+        return self._ymin
+
+    @property
+    def xmax(self):
+        return self._xmax
+
+    @property
+    def ymax(self):
+        return self._ymax
+
+    @property
+    def transform(self):
+        return self._transform
+
+    @property
+    def inverse_transform(self):
+        return self._inverse_transform
+
+    @property
+    def n_rows(self):
+        return self._n_rows
+
+    @property
+    def n_cols(self):
+        return self._n_cols