Point-in-Polygon for the z=0 projection. This function enhances
        the performance of ``Polygon.contains()`` by verifying only the 
        points which are inside the bounding box of the polygon. To do 
        it fast, it needs the points array to be already sorted by one
        column.
        
        :param points: list of *(x, y, z) or (x, y)* coordinates of the
            points to check. (The z value will not be taken into 
            account).
        :type points: ndarray (shape=(N, 2 or 3))
        :param sorted_col: Index of the sorted column (0 or 1).
        :type sorted_col: int
        :param radius: Enlarge Polygons domain by a specified quantity.
        :type radius: float
        :returns: Which points are inside the polygon.
        :rtype: ndarray (dtpye=bool)
        
        .. warning:: By default pip considers that the set of points is
            currently sorted by the first column.
        .. warning:: This method only works if the polygon has been 
            locked (:func:`lock`).

def pip(self, points, sorted_col=0, radius=0):
        """
        Point-in-Polygon for the z=0 projection. This function enhances
        the performance of ``Polygon.contains()`` by verifying only the 
        points which are inside the bounding box of the polygon. To do 
        it fast, it needs the points array to be already sorted by one
        column.
        
        :param points: list of *(x, y, z) or (x, y)* coordinates of the
            points to check. (The z value will not be taken into 
            account).
        :type points: ndarray (shape=(N, 2 or 3))
        :param sorted_col: Index of the sorted column (0 or 1).
        :type sorted_col: int
        :param radius: Enlarge Polygons domain by a specified quantity.
        :type radius: float
        :returns: Which points are inside the polygon.
        :rtype: ndarray (dtpye=bool)
        
        .. warning:: By default pip considers that the set of points is
            currently sorted by the first column.
        .. warning:: This method only works if the polygon has been 
            locked (:func:`lock`).
        """
        xy = points[:, :2]
        n_points = xy.shape[0]
        index = np.arange(n_points, dtype = int)
        b = self.domain
        b[0] = b[0] - radius
        b[1] = b[1] + radius
        
        # Slicing the sorted column
        k = np.searchsorted(xy[:, sorted_col],
                            (b[0, sorted_col], b[1, sorted_col]+1e-10))
        xy = xy[k[0]:k[1]]
        index = index[k[0]:k[1]]
        
        # solution
        k = index[self.path.contains_points(xy, radius=radius)]
        sol = np.zeros(n_points, dtype=bool)
        sol[k] = True
                
        return sol