Largest closed contour + mask of the interior

gistfile1.txt

def get_largest_clcontour(cset):    
    """
    This function iterates over all contours in a plt.contour() instance to find the largest closed one. Found in:
    
    https://stackoverflow.com/questions/38249949/finding-largest-closed-3d-contour-plot
    
    Parameters:
    -----------
        cset: plt.contour() instance with the defined levels over which this function will iterate.
   
    Returns:
    --------
        coords: [:, 3]. In [:, 0] and [:, 1] are the x-axis and y-axis coordinates respectively of the largest closed contour. [:, 2] is the level.
    """
    
    maxsize = 0
    
    # Iterate over all the contour segments and find the largest
    for i, segs in enumerate(cset.allsegs):
        for j, seg in enumerate(segs):
            # First make sure it's closed
            if (seg[0]-seg[-1]).any():
                continue
            # Now get it's size
            size = distance_matrix(seg, seg).max()
            if size > maxsize:
                maxsize = size
                maxseg = (i, j)
                segment = seg

    # Now highlight the "biggest" closed contour
    cset.collections[maxseg[0]].set_color('r')
    cset.collections[maxseg[0]].set_lw(5)
    pts2d = cset.allsegs[maxseg[0]][maxseg[1]]
    z = cset.levels[maxseg[0]]
    coords = np.c_[ pts2d, z*np.ones(pts2d.shape[0]) ]
    
    return coords

def interior_mask(variable, coords):
    """
    Masks a variable using the coordinates of a closed contour obtained with get_largest_clcontour(). Masks the exterior, but can be easily modified to mask the interior.
    Parameters:
    -----------
        variable: 2-d variable that was used to find the largest closed contour.
        coords: the coordinates of said contour, output of get_largest_contour()
   
    Returns:
    --------
        variable_masked: variable with all gridpoints outside the contour defined by coords masked
    """
    
    mask = np.full(np.shape(variable), True)
    polygon = Polygon(list(coords[:, 0:2]))
    for lat in range(0, len(mask[:, 0]), 1):
        for lon in range(0, len(mask[0, :]), 1):
            pnt = Point(variable['x-coordinate'][lon], variable['y-coordinate'][lat])
            if polygon.contains(pnt):
                mask[lat, lon] = False
    variable_mask = np.ma.array(variable, mask = mask)
    
    return variable_mask