RCLV identification

[hscannell]

Using the particle trajectory information from each ensemble member, use the RCLV method to create a labeled mask of coherent eddies.

# Eddy detection

lav_abs = np.abs(particle_history[:,2].copy())  # absolute value of vorticity on particles
lx,ly = particle_history[:,0].copy(),particle_history[:,1].copy() # particle trjectories

# lq = np.zeros([Nhold,Npart])
# for i in range(Nhold):
#     lq[i] = lpa.interpolate_gridded_scalar(lx[i], ly[i], q[i,0]) # get the PV on particles (if needed)
    
lx.shape = (Nhold, pNy, pNx)
ly.shape = (Nhold, pNy, pNx)
lav_abs.shape = (Nhold, pNy, pNx)
# lq.shape = (Nhold, pNy, pNx)

lx1 = np.unwrap(lx[::-1]*2*pi/m.W, axis=0)*m.W/(2*pi)  # unwrap trajectories to correct the jump of displacement at the periodic boundaries
ly1 = np.unwrap(ly[::-1]*2*pi/m.L, axis=0)*m.L/(2*pi)
lx1 = lx1[::-1]  # switch back to backward timing
ly1 = ly1[::-1]

ndays1 = 60  # lifetime of RCLV
lavd1 = lav_abs[-ndays1:].mean(axis=0)  #LAVD field 
del(lav_abs)

cov1 = 10.0  # convexity deficincy threshold, here uses 10.0 to disable it
CI = -0.75   # threshold for coherency index
min_area = 200  # minimum area (number of pixels) of RCLV

kwargs = dict(min_distance=20,  # minimum distance between LAVD maxima (pixel units)
              max_width=100,    # maximum width of the search window (pixel units)
              CI_th=CI, CI_tol=0.01, convex_def=cov1, convex_def_tol=0.001, 
              init_contour_step_frac=0.64,  # the value with which to increment the initial contour level
              min_limit_diff=1e-10, min_area=min_area, periodic=(True, True))

contours1 = list(rclv.find_convex_contours(lavd1, 
                                           lx1[-1:-ndays1-1:-1], ly1[-1:-ndays1-1:-1], 
                                           dx=dx, dy=dy, # particle spacing, must have the same unit as lx1 and ly1
                                           **kwargs))

mask = rclv.label_points_in_contours(lavd1.shape, [c[1] for c in contours1])  #label particles in RCLVs