Merge pull request #27 from jk-rieck/26-detection-of-horseshoe-shapes-needs-refinement

26 detection of horseshoe shapes needs refinement

Author: jk-rieck

eddytools/detection.py

@@ -377,18 +377,20 @@ def detect_OW_core(data, det_param, OW, vort, t, OW_thr, e1f, e2f,
             del eddi[e]
             continue
         if (det_param['no_long'] or det_param['no_two']):
-            min_width = int(np.around(np.sqrt(region_Npix / np.pi)))
+            min_width = int(np.floor(np.sqrt(region_Npix / np.pi)))
             X_cen = int(np.around(np.mean(index[1])))
             Y_cen = int(np.around(np.mean(index[0])))
         if det_param['no_two']:
             if len(np.shape(data[det_param['OW_thr_name']])) > 1:
                 peak_thr = OW_thr.values[interior].mean()
             else:
                 peak_thr = OW_thr
-            X_peaks = len(find_peaks(-OW.isel(time=t).values[Y_cen,
-                          iimin:iimax], height=-peak_thr)[0])
-            Y_peaks = len(find_peaks(-OW.isel(time=t).values[ijmin:ijmax,
-                          X_cen], height=-peak_thr)[0])
+            X_peak_info = find_peaks(-OW.isel(time=t).values[Y_cen,
+                              iimin:iimax], height=-peak_thr)
+            Y_peak_info = find_peaks(-OW.isel(time=t).values[ijmin:ijmax,
+                              X_cen], height=-peak_thr)
+            X_peaks = len(X_peak_info[0])
+            Y_peaks = len(Y_peak_info[0])
             if ((X_peaks > 1) | (Y_peaks > 1)):
                 Ypix_cen1 = get_width(OW.isel(time=t).values[ijmin:ijmax,
                                 X_cen], peak_thr)
@@ -398,14 +400,43 @@ def detect_OW_core(data, det_param, OW, vort, t, OW_thr, e1f, e2f,
                                 iimin:iimax], peak_thr)
                 Xpix_cen2 = get_width(OW.isel(time=t).values[Y_cen,
                                 iimax-1:iimin2:-1], peak_thr)
-                eddy_no_horseshoe = (((Xpix_cen1 > min_width)
-                                    & (Ypix_cen1 > min_width)) |
-                                     ((Xpix_cen2 > min_width)
-                                    & (Ypix_cen2 > min_width)) |
-                                     ((Xpix_cen2 > min_width)
-                                    & (Ypix_cen1 > min_width)) |
-                                     ((Xpix_cen1 > min_width)
-                                    & (Ypix_cen2 > min_width)))
+                tmp_no_horseshoe = (((Xpix_cen1 > min_width)
+                                   & (Ypix_cen1 > min_width)) |
+                                    ((Xpix_cen2 > min_width)
+                                   & (Ypix_cen2 > min_width)) |
+                                    ((Xpix_cen2 > min_width)
+                                   & (Ypix_cen1 > min_width)) |
+                                    ((Xpix_cen1 > min_width)
+                                   & (Ypix_cen2 > min_width)))
+                if (X_peaks > 1):
+                    X_peak1 = iimin + X_peak_info[0][0]
+                    X_peak2 = iimin + X_peak_info[0][1]
+                    Xmin = -np.max(X_peak_info[1]['peak_heights'])
+                    Xmax = np.max(OW.isel(time=t).values[Y_cen,
+                                                         X_peak1:X_peak2])
+                    Xdiff = Xmin - Xmax
+                    Xpix_ratio = Xpix_cen1 / Xpix_cen2
+                    if ((Xpix_ratio > 2.5) | (Xpix_ratio < 0.4)):
+                        Xdiff = 0
+                else:
+                    Xdiff = 0
+                    Xmin = 1
+                if (Y_peaks > 1):
+                    Y_peak1 = ijmin + Y_peak_info[0][0]
+                    Y_peak2 = ijmin + Y_peak_info[0][1]
+                    Ymin = -np.max(Y_peak_info[1]['peak_heights'])
+                    Ymax = np.max(OW.isel(time=t).values[Y_peak1:Y_peak2,
+                                                         X_cen])
+                    Ydiff = Ymin - Ymax
+                    Ypix_ratio = Ypix_cen1 / Ypix_cen2
+                    if ((Ypix_ratio > 2.5) | (Ypix_ratio < 0.4)):
+                        Ydiff = 0
+                else:
+                    Ydiff = 0
+                    Ymin = 1
+                Xdiff_small = np.abs(Xdiff) < 0.75 * np.abs(Xmin)
+                Ydiff_small = np.abs(Ydiff) < 0.75 * np.abs(Ymin)
+                eddy_no_horseshoe = tmp_no_horseshoe & Xdiff_small & Ydiff_small
             else:
                 eddy_no_horseshoe = True
         else: