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
Using the particle trajectory information from each ensemble member, use the RCLV method to create a labeled mask of coherent eddies.