def METHOD_NAME(
    pts: Coordinates,
    N: int = 100,
    margin: float = 0.4,
    margin_x: float | None = None,
    margin_y: float | None = None,
    nb_pixels_x: int | None = None,
    nb_pixels_y: int | None = None,
    min_pixel_size: float = 0.4,
    snap_res: float = 0.05,
) -> Coordinates:
    """Pixelates a shape (as 2d array) onto an NxN grid.
    Arguments:
        pts: The 2D array to be pixelated.
        N: The number of pixels on an edge of the grid.
    Returns:
        A list of pixel bounding boxes
    """
    shape = points_to_shapely(pts)  # convert to shapely
    if not shape:
        return []
    if margin_x is None:
        margin_x = margin
    if margin_y is None:
        margin_y = margin
    if nb_pixels_x is None:
        nb_pixels_x = N
    if nb_pixels_y is None:
        nb_pixels_y = N
    west, south, east, north = shape.bounds
    width = east - west + 2 * margin_x
    height = north - south + 2 * margin_y
    if min_pixel_size is not None:
        max_nb_pixels_x = int(np.ceil(width / min_pixel_size))
        nb_pixels_x = min(nb_pixels_x, max_nb_pixels_x)
    if min_pixel_size is not None:
        max_nb_pixels_y = int(np.ceil(height / min_pixel_size))
        nb_pixels_y = min(nb_pixels_y, max_nb_pixels_y)
    w = width / nb_pixels_x + snap_res
    h = height / nb_pixels_y + snap_res
    ax = margin_x
    ay = margin_y
    xs = np.linspace(west + w / 2 - ax, east - w / 2 + ax, nb_pixels_x)
    ys = np.linspace(south + h / 2 - ay, north - h / 2 + ay, nb_pixels_y)
    xs = _snap_to_resolution(xs, snap_res)
    ys = _snap_to_resolution(ys, snap_res)
    grid = it.product(xs, ys)
    pixels = []
    for x, y in grid:
        _w, _s, _e, _n = x - w / 2, y - h / 2, x + w / 2, y + h / 2
        newpix = geometry.box(_w, _s, _e, _n)
        if shape.intersects(newpix):
            r = (_w, _s, _e, _n)
            pixels.append(r)
    return pixels