Several update:

- Fix some bug around x bounds
- Numba replace more cython function
- Merge Grid dataset class
- ...

Author: AntSimi

src/py_eddy_tracker/dataset/grid.py

@@ -27,11 +27,9 @@
 """
 
 import logging
-from numpy import where, empty, array, concatenate, ma, zeros, unique, round, ones
-from scipy.ndimage import minimum_filter
+from numpy import empty, array, concatenate, ma, zeros, unique, round, ones, int_
 from matplotlib.figure import Figure
 from numba import njit, types as numba_types
-from .tools import index_from_nearest_path_with_pt_in_bbox
 
 
 class Amplitude(object):
@@ -468,7 +466,7 @@ def get_index_nearest_path_bbox_contain_pt(self, level, xpt, ypt):
         
         overhead of python is huge with numba, cython little bit best??
         """
-        index = index_from_nearest_path_with_pt_in_bbox(
+        index = index_from_nearest_path_with_pt_in_bbox_(
             level,
             self.level_index,
             self.nb_contour_per_level,
@@ -483,7 +481,7 @@ def get_index_nearest_path_bbox_contain_pt(self, level, xpt, ypt):
             xpt,
             ypt
             )
-        if index is None:
+        if index == -1:
             return None
         else:
             return self.contours.collections[level]._paths[index]
@@ -500,7 +498,7 @@ def display(self, ax, **kwargs):
         ax.autoscale_view()
 
 
-@njit(cache=True)
+@njit(cache=True, fastmath=True)
 def index_from_nearest_path_with_pt_in_bbox_(
         level_index,
         l_i,
@@ -537,7 +535,15 @@ def index_from_nearest_path_with_pt_in_bbox_(
     # We iterate over contour in the same level
     for i_elt_c in range(i_start_c, i_end_c):
         # if bbox of contour doesn't contain pt, we skip this contour
-        if y_min_per_c[i_elt_c] > ypt or y_max_per_c[i_elt_c] < ypt or x_min_per_c[i_elt_c] > xpt or x_max_per_c[i_elt_c] < xpt:
+        if y_min_per_c[i_elt_c] > ypt:
+            continue
+        if y_max_per_c[i_elt_c] < ypt:
+            continue
+        x_min = x_min_per_c[i_elt_c]
+        xpt_ = (xpt - x_min) % 360 + x_min
+        if x_min > xpt_:
+            continue
+        if x_max_per_c[i_elt_c] < xpt_:
             continue
         # Indice of first pt of contour
         i_start_pt = indices_of_first_pts[i_elt_c]
@@ -549,13 +555,15 @@ def index_from_nearest_path_with_pt_in_bbox_(
         # We do iteration on pt to check dist, if it's inferior we store
         # index of contour
         for i_elt_pt in range(i_start_pt, i_end_pt):
-            dist = (x_value[i_elt_pt] - xpt) ** 2 + (y_value[i_elt_pt] - ypt) ** 2
+            d_x = x_value[i_elt_pt] - xpt_
+            if abs(d_x) > 180:
+                d_x = (d_x + 180) % 360 - 180
+            dist = d_x ** 2 + (y_value[i_elt_pt] - ypt) ** 2
             if dist < dist_ref:
                 dist_ref = dist
                 i_ref = i_elt_c
     # No iteration on contour, we return no index of contour
     if find_contour == 0:
-        return -1
+        return int_(-1)
     # We return index of contour, for the specific level
-    return i_ref - i_start_c
-
+    return int_(i_ref - i_start_c)

src/py_eddy_tracker/eddy_feature.py

@@ -1,20 +1,20 @@
 from ..observations import EddiesObservations as Model
-from ..observations import GridDataset
+from ..dataset.grid import RegularGridDataset
 from numpy import where, meshgrid, bincount, ones, unique, bool_, arange
-import logging
+from numba import njit
 from os import path
 
 
 class CheltonTracker(Model):
-    GROUND = GridDataset(path.join(path.dirname(__file__), '/data/adelepoulle/Test/Test_eddy/20180220_high_res_mask/mask_1_60.nc'), 'lon', 'lat')
+    GROUND = RegularGridDataset(path.join(path.dirname(__file__), '/data/adelepoulle/Test/Test_eddy/20180220_high_res_mask/mask_1_60.nc'), 'lon', 'lat')
 
     @staticmethod
     def cost_function(records_in, records_out, distance):
         """We minimize on distance between two obs
         """
         return distance
 
-    def mask_function(self, other):
+    def mask_function(self, other, distance):
         """We mask link with ellips and ratio
         """
         # Compute Parameter of ellips
@@ -34,18 +34,8 @@ def mask_function(self, other):
             minor=minor, # Minor can be bigger than major??
             major=y)
 
-        # We check ratio
-        other_amplitude, self_amplitude = meshgrid(
-            other.obs['amplitude'], self.obs['amplitude'])
-        other_radius, self_radius = meshgrid(
-            other.obs['radius_e'], self.obs['radius_e'])
-        ratio_amplitude = other_amplitude / self_amplitude
-        ratio_radius = other_radius / self_radius
-        mask *= \
-            (ratio_amplitude < 2.5) * \
-            (ratio_amplitude > 1 / 2.5) * \
-            (ratio_radius < 2.5) * \
-            (ratio_radius > 1 / 2.5)
+        # We check ratio (maybe not usefull)
+        check_ratio(mask, self.obs['amplitude'], other.obs['amplitude'], self.obs['radius_e'], other.obs['radius_e'])
         indexs_closest = where(mask)
         mask[indexs_closest] = self.across_ground(self.obs[indexs_closest[0]], other.obs[indexs_closest[1]])
         return mask
@@ -87,3 +77,29 @@ def post_process_link(self, other, i_self, i_other):
             i_self = i_self[mask]
             i_other = i_other[mask]
         return i_self, i_other
+
+
+@njit(cache=True)
+def check_ratio(current_mask, self_amplitude, other_amplitude, self_radius, other_radius):
+    """
+    Only very few case are remove with selection
+    :param current_mask:
+    :param self_amplitude:
+    :param other_amplitude:
+    :param self_radius:
+    :param other_radius:
+    :return:
+    """
+    self_shape, other_shape = current_mask.shape
+    r_min = 1 / 2.5
+    r_max = 2.5
+    for i in range(self_shape):
+        for j in range(other_shape):
+            if current_mask[i, j]:
+                r_amplitude = other_amplitude[j] / self_amplitude[i]
+                if r_amplitude >= r_max or r_amplitude <= r_min:
+                    current_mask[i, j] = False
+                    continue
+                r_radius = other_radius[j] / self_radius[i]
+                if r_radius >= r_max or r_radius <= r_min:
+                    current_mask[i, j] = False

src/py_eddy_tracker/featured_tracking/old_tracker_reference.py

@@ -0,0 +1,211 @@
+# -*- coding: utf-8 -*-
+"""
+"""
+from numpy import sin, pi, cos, arctan2, empty, nan, absolute, floor, ones, linspace, interp
+from numba import njit, prange
+from numpy.linalg import lstsq
+
+
+@njit(cache=True, fastmath=True, parallel=False)
+def distance_grid(lon0, lat0, lon1, lat1):
+    """
+    Args:
+        lon0:
+        lat0:
+        lon1:
+        lat1:
+
+    Returns:
+        nan value for far away point, and km for other
+    """
+    nb_0 = lon0.shape[0]
+    nb_1 = lon1.shape[0]
+    dist = empty((nb_0, nb_1))
+    D2R = pi / 180.
+    for i in prange(nb_0):
+        for j in prange(nb_1):
+            dlat = absolute(lat1[j] - lat0[i])
+            if dlat > 15:
+                dist[i, j] = nan
+                continue
+            dlon = absolute(lon1[j] - lon0[i])
+            if dlon > 180:
+                dlon = absolute((dlon + 180) % 360 - 180)
+            if dlon > 20:
+                dist[i, j] = nan
+                continue
+            sin_dlat = sin((dlat) * 0.5 * D2R)
+            sin_dlon = sin((dlon) * 0.5 * D2R)
+            cos_lat1 = cos(lat0[i] * D2R)
+            cos_lat2 = cos(lat1[j] * D2R)
+            a_val = sin_dlon ** 2 * cos_lat1 * cos_lat2 + sin_dlat ** 2
+            dist[i, j] = 6370.997 * 2 * arctan2(a_val ** 0.5, (1 - a_val) ** 0.5)
+    return dist
+
+
+@njit(cache=True, fastmath=True)
+def distance(lon0, lat0, lon1, lat1):
+    D2R = pi / 180.
+    sin_dlat = sin((lat1 - lat0) * 0.5 * D2R)
+    sin_dlon = sin((lon1 - lon0) * 0.5 * D2R)
+    cos_lat1 = cos(lat0 * D2R)
+    cos_lat2 = cos(lat1 * D2R)
+    a_val = sin_dlon ** 2 * cos_lat1 * cos_lat2 + sin_dlat ** 2
+    return 6370997.0 * 2 * arctan2(a_val ** 0.5, (1 - a_val) ** 0.5)
+
+
+@njit(cache=True)
+def distance_vincenty(lon0, lat0, lon1, lat1):
+    """ better than haversine but buggy ??"""
+    D2R = pi / 180.
+    dlon = (lon1 - lon0) * D2R
+    cos_dlon = cos(dlon)
+    cos_lat1 = cos(lat0 * D2R)
+    cos_lat2 = cos(lat1 * D2R)
+    sin_lat1 = sin(lat0 * D2R)
+    sin_lat2 = sin(lat1 * D2R)
+    return 6370997.0 * arctan2(
+        ((cos_lat2 * sin(dlon) ** 2) + (cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_dlon) ** 2) ** .5,
+        sin_lat1 * sin_lat2 + cos_lat1 * cos_lat2 * cos_dlon)
+
+
+@njit(cache=True, fastmath=True)
+def interp2d_geo(x_g, y_g, z_g, m_g, x, y):
+    """For geographic grid, test of cicularity
+    Maybe test if we are out of bounds
+    """
+    x_ref = x_g[0]
+    y_ref = y_g[0]
+    x_step = x_g[1] - x_ref
+    y_step = y_g[1] - y_ref
+    nb_x = x_g.shape[0]
+    is_circular = (x_g[-1] + x_step) % 360 == x_g[0] % 360
+    z = empty(x.shape)
+    for i in prange(x.size):
+        x_ = (x[i] - x_ref) / x_step
+        y_ = (y[i] - y_ref) / y_step
+        i0 = int(floor(x_))
+        i1 = i0 + 1
+        if is_circular:
+            xd = (x_ - i0)
+            i0 %= nb_x
+            i1 %= nb_x
+        j0 = int(floor(y_))
+        j1 = j0 + 1
+        yd = (y_ - j0)
+        z00 = z_g[i0, j0]
+        z01 = z_g[i0, j1]
+        z10 = z_g[i1, j0]
+        z11 = z_g[i1, j1]
+        if m_g[i0, j0] or m_g[i0, j1] or m_g[i1, j0] or m_g[i1, j1]:
+            z[i] = nan
+        else:
+            z[i] = (z00 * (1 - xd) + (z10 * xd)) * (1 - yd) + (z01 * (1 - xd) + z11 * xd) * yd
+    return z
+
+@njit(cache=True, fastmath=True, parallel=True)
+def custom_convolution(data, mask, kernel):
+    """do sortin at high lattitude big part of value are masked"""
+    nb_x = kernel.shape[0]
+    demi_x = int((nb_x - 1) / 2)
+    demi_y = int((kernel.shape[1] - 1) / 2)
+    out = empty(data.shape[0] - nb_x + 1)
+    for i in prange(out.shape[0]):
+        if mask[i + demi_x, demi_y] == 1:
+            w = (mask[i:i + nb_x] * kernel).sum()
+            if w != 0:
+                out[i] = (data[i:i + nb_x] * kernel).sum() / w
+            else:
+                out[i] = nan
+        else:
+            out[i] = nan
+    return out
+
+
+@njit(cache=True)
+def fit_circle(x_vec, y_vec):
+    nb_elt = x_vec.shape[0]
+    p_inon_x = empty(nb_elt)
+    p_inon_y = empty(nb_elt)
+
+    x_mean = x_vec.mean()
+    y_mean = y_vec.mean()
+
+    norme = (x_vec - x_mean) ** 2 + (y_vec - y_mean) ** 2
+    norme_max = norme.max()
+    scale = norme_max ** .5
+
+    # Form matrix equation and solve it
+    # Maybe put f4
+    datas = ones((nb_elt, 3))
+    datas[:, 0] = 2. * (x_vec - x_mean) / scale
+    datas[:, 1] = 2. * (y_vec - y_mean) / scale
+
+    (center_x, center_y, radius), residuals, rank, s = lstsq(datas, norme / norme_max)
+
+    # Unscale data and get circle variables
+    radius += center_x ** 2 + center_y ** 2
+    radius **= .5
+    center_x *= scale
+    center_y *= scale
+    # radius of fitted circle
+    radius *= scale
+    # center X-position of fitted circle
+    center_x += x_mean
+    # center Y-position of fitted circle
+    center_y += y_mean
+
+    # area of fitted circle
+    c_area = (radius ** 2) * pi
+    # Find distance between circle center and contour points_inside_poly
+    for i_elt in range(nb_elt):
+        # Find distance between circle center and contour points_inside_poly
+        dist_poly = ((x_vec[i_elt] - center_x) ** 2 + (y_vec[i_elt] - center_y) ** 2) ** .5
+        # Indices of polygon points outside circle
+        # p_inon_? : polygon x or y points inside & on the circle
+        if dist_poly > radius:
+            p_inon_y[i_elt] = center_y + radius * (y_vec[i_elt] - center_y) / dist_poly
+            p_inon_x[i_elt] = center_x - (center_x - x_vec[i_elt]) * (center_y - p_inon_y[i_elt]) / (
+                        center_y - y_vec[i_elt])
+        else:
+            p_inon_x[i_elt] = x_vec[i_elt]
+            p_inon_y[i_elt] = y_vec[i_elt]
+
+    # Area of closed contour/polygon enclosed by the circle
+    p_area_incirc = 0
+    p_area = 0
+    for i_elt in range(nb_elt - 1):
+        # Indices of polygon points outside circle
+        # p_inon_? : polygon x or y points inside & on the circle
+        p_area_incirc += p_inon_x[i_elt] * p_inon_y[1 + i_elt] - p_inon_x[i_elt + 1] * p_inon_y[i_elt]
+        # Shape test
+        # Area and centroid of closed contour/polygon
+        p_area += x_vec[i_elt] * y_vec[1 + i_elt] - x_vec[1 + i_elt] * y_vec[i_elt]
+    p_area = abs(p_area) * .5
+    p_area_incirc = abs(p_area_incirc) * .5
+
+    a_err = (c_area - 2 * p_area_incirc + p_area) * 100. / c_area
+    return center_x, center_y, radius, a_err
+
+
+@njit(cache=True, fastmath=True)
+def uniform_resample(x_val, y_val, num_fac=2, fixed_size=None):
+    """
+    Resample contours to have (nearly) equal spacing
+       x_val, y_val    : input contour coordinates
+       num_fac : factor to increase lengths of output coordinates
+    """
+    # Get distances
+    dist = empty(x_val.shape)
+    dist[0] = 0
+    dist[1:] = distance(x_val[:-1], y_val[:-1], x_val[1:], y_val[1:])
+    # To be still monotonous (dist is store in m)
+    dist[1:][dist[1:]<1e-3] = 1e-3
+    dist = dist.cumsum()
+    # Get uniform distances
+    if fixed_size is None:
+        fixed_size = dist.size * num_fac
+    d_uniform = linspace(0, dist[-1], fixed_size)
+    x_new = interp(d_uniform, dist, x_val)
+    y_new = interp(d_uniform, dist, y_val)
+    return x_new, y_new