How to use the mercantile.xy_bounds function in mercantile

arr = numpy.transpose(arr, [2, 0, 1])

                    tbounds = mercantile.xy_bounds(leftTile)
                    data, mask = tile_read(
                        "cogeo.tif", tbounds, 256, resampling_method="nearest"
                    )
                    numpy.testing.assert_array_equal(data, arr)

                    # Bottom right tile
                    mime_type, tile = cog.get_tile(4, 3, 0)
                    tile_length = 256 * 256 * 3
                    t = struct.unpack_from("{}b".format(tile_length), tile)
                    arr = numpy.array(t).reshape(256, 256, 3).astype(numpy.uint8)
                    arr = numpy.transpose(arr, [2, 0, 1])

                    tbounds = mercantile.xy_bounds(rightTile)
                    data, mask = tile_read(
                        "cogeo.tif", tbounds, 256, resampling_method="nearest"
                    )
                    numpy.testing.assert_array_equal(data, arr)

                    # Low resolution (overview 1)
                    # Top Left tile
                    # NOTE: overview internal tile size is 128px
                    # We need to stack two internal tiles to compare with
                    # the 256px mercator tile fetched by rio-tiler
                    # ref: https://github.com/cogeotiff/rio-cogeo/issues/60
                    mime_type, tile = cog.get_tile(1, 0, 1)
                    tile_length = 128 * 128 * 3
                    t = struct.unpack_from("{}b".format(tile_length), tile)
                    arr1 = numpy.array(t).reshape(128, 128, 3).astype(numpy.uint8)
                    arr1 = numpy.transpose(arr1, [2, 0, 1])

def read_tile(self, z, x, y):
        """Read raster tile data and mask."""
        mercator_tile = mercantile.Tile(x=x, y=y, z=z)
        tile_bounds = mercantile.xy_bounds(mercator_tile)

        data, mask = tile_read(
            self.path,
            tile_bounds,
            self.tiles_size,
            indexes=self.indexes,
            nodata=self.nodata,
        )
        data = (data[0] + data[1]) / 2
        return data.astype(numpy.uint8), mask

im = Image.fromarray(img_arr)
        x = 0
        y = 0
        im_path = join(self.tmp_dir, '{}/{}/{}.png'.format(zoom, x, y))
        make_dir(im_path, use_dirname=True)
        im.save(im_path)

        tile_schema = join(self.tmp_dir, '{z}/{x}/{y}.png')

        # Get NW and SE corner of central half of tile.
        nw_bounds = mercantile.xy_bounds(0, 0, zoom)
        nw_merc_y = nw_bounds.top - (nw_bounds.top - nw_bounds.bottom) / 4
        nw_merc_x = nw_bounds.left + (nw_bounds.right - nw_bounds.left) / 4
        nw_lng, nw_lat = merc2lnglat(nw_merc_x, nw_merc_y)

        se_bounds = mercantile.xy_bounds(0, 0, zoom)
        se_merc_y = se_bounds.bottom + (se_bounds.top - se_bounds.bottom) / 4
        se_merc_x = se_bounds.right - (se_bounds.right - se_bounds.left) / 4
        se_lng, se_lat = merc2lnglat(se_merc_x, se_merc_y)

        # min_lat, min_lng, max_lat, max_lng = bounds
        bounds = [se_lat, nw_lng, nw_lat, se_lng]
        output_uri = join(self.tmp_dir, 'output.tif')
        _zxy2geotiff(tile_schema, zoom, bounds, output_uri)

        with rasterio.open(output_uri) as dataset:
            tiff_arr = dataset.read()
            self.assertEqual(tiff_arr.shape, (3, 128, 128))
            exp_arr = np.transpose(img_arr, (2, 0, 1))[:, 64:-64, 64:-64]
            np.testing.assert_array_equal(tiff_arr, exp_arr)

# ref: https://github.com/cogeotiff/rio-cogeo/issues/60
                    mime_type, tile = cog.get_tile(1, 0, 1)
                    tile_length = 128 * 128 * 3
                    t = struct.unpack_from("{}b".format(tile_length), tile)
                    arr1 = numpy.array(t).reshape(128, 128, 3).astype(numpy.uint8)
                    arr1 = numpy.transpose(arr1, [2, 0, 1])

                    mime_type, tile = cog.get_tile(2, 0, 1)
                    tile_length = 128 * 128 * 3
                    t = struct.unpack_from("{}b".format(tile_length), tile)
                    arr2 = numpy.array(t).reshape(128, 128, 3).astype(numpy.uint8)
                    arr2 = numpy.transpose(arr2, [2, 0, 1])
                    arr = numpy.dstack((arr1, arr2))

                    lowTile = mercantile.Tile(118594, 60034, 17)
                    tbounds = mercantile.xy_bounds(lowTile)
                    data, mask = tile_read(
                        "cogeo.tif", tbounds, 256, resampling_method="nearest"
                    )
                    data = data[:, 128:, :]
                    numpy.testing.assert_array_equal(data, arr)

blocks = list(set(out_dst.block_shapes))
                assert len(blocks) == 1
                ts = blocks[0][0]
                assert not out_dst.width % ts
                assert not out_dst.height % ts
                max_zoom = get_max_zoom(out_dst)

                bounds = list(
                    transform_bounds(
                        *[src_dst.crs, "epsg:4326"] + list(src_dst.bounds),
                        densify_pts=21
                    )
                )

                leftTile = mercantile.tile(bounds[0], bounds[3], max_zoom)
                tbounds = mercantile.xy_bounds(leftTile)
                west, north = tbounds.left, tbounds.top
                assert out_dst.transform.xoff == west
                assert out_dst.transform.yoff == north

                rightTile = mercantile.tile(bounds[2], bounds[1], max_zoom)
                tbounds = mercantile.xy_bounds(rightTile)
                east, south = tbounds.right, tbounds.bottom

                lrx = round(
                    out_dst.transform.xoff + out_dst.transform.a * out_dst.width, 6
                )
                lry = round(
                    out_dst.transform.yoff + out_dst.transform.e * out_dst.height, 6
                )
                assert lrx == round(east, 6)
                assert lry == round(south, 6)

def _expand_bounds(self, bounds):
        if bounds is None:
            return bounds
        min_tile_x, min_tile_y, max_tile_x, max_tile_y = self._tile_coords(bounds)

        ul = box(*mercantile.xy_bounds(mercantile.Tile(z=self.zoom_level, x=min_tile_x, y=max_tile_y)))
        lr = box(*mercantile.xy_bounds(mercantile.Tile(z=self.zoom_level, x=max_tile_x, y=min_tile_y)))

        return ul.union(lr).bounds

def merc2pixel(tile_x, tile_y, zoom, merc_x, merc_y, tile_sz=256):
    """Convert Web Mercator point to pixel coordinates.

    This is within the coordinate frame of a single ZXY tile.

    Args:
        tile_x: (int) x coordinate of ZXY tile
        tile_y: (int) y coordinate of ZXY tile
        zoom: (int) zoom level of ZXY tile
        merc_x: (float) Web Mercator x axis of point
        merc_y: (float) Web Mercator y axis of point
        tile_sz: (int) size of ZXY tile
    """
    tile_merc_bounds = mercantile.xy_bounds(tile_x, tile_y, zoom)
    pix_y = int(
        round(tile_sz * ((tile_merc_bounds.top - merc_y) /
                         (tile_merc_bounds.top - tile_merc_bounds.bottom))))
    pix_x = int(
        round(tile_sz * ((merc_x - tile_merc_bounds.left) /
                         (tile_merc_bounds.right - tile_merc_bounds.left))))
    return (pix_x, pix_y)

"""Render a tile into Web Mercator.

    Arguments:
        tile {mercantile.Tile} -- Tile to render.
        sources {list} -- Sources to render from.

    Keyword Arguments:
        transformation {Transformation} -- Transformation to apply. (default: {None})
        format {function} -- Output format. (default: {None})
        scale {int} -- Output scale factor. (default: {1})
        expand {bool} -- Whether to expand single-band, paletted sources to RGBA. (default: {True})

    Returns:
        (dict, bytes) -- Tuple of HTTP headers (dict) and bytes.
    """
    bounds = Bounds(mercantile.xy_bounds(tile), WEB_MERCATOR_CRS)
    shape = tuple(map(int, Affine.scale(scale) * TILE_SHAPE))

    return render(
        bounds,
        shape,
        WEB_MERCATOR_CRS,
        sources=sources,
        format=format,
        transformation=transformation,
        expand=expand,
    )

if band not in LANDSAT_BANDS:
            raise InvalidBandName("{} is not a valid Landsat band name".format(band))

    scene_params = _landsat_parse_scene_id(sceneid)
    meta_data = _landsat_get_mtl(sceneid).get("L1_METADATA_FILE")
    landsat_address = "{}/{}".format(LANDSAT_BUCKET, scene_params["key"])

    wgs_bounds = toa_utils._get_bounds_from_metadata(meta_data["PRODUCT_METADATA"])

    if not utils.tile_exists(wgs_bounds, tile_z, tile_x, tile_y):
        raise TileOutsideBounds(
            "Tile {}/{}/{} is outside image bounds".format(tile_z, tile_x, tile_y)
        )

    mercator_tile = mercantile.Tile(x=tile_x, y=tile_y, z=tile_z)
    tile_bounds = mercantile.xy_bounds(mercator_tile)

    def _tiler(band):
        address = "{}_B{}.TIF".format(landsat_address, band)
        if band == "QA":
            nodata = 1
        else:
            nodata = 0

        return utils.tile_read(
            address, bounds=tile_bounds, tilesize=tilesize, nodata=nodata, **kwargs
        )

    with futures.ThreadPoolExecutor(max_workers=MAX_THREADS) as executor:
        data, masks = zip(*list(executor.map(_tiler, bands)))
        data = np.concatenate(data)
        mask = np.all(masks, axis=0).astype(np.uint8) * 255

"""Burn tile with features.

    Args:
      tile: the mercantile tile to burn.
      features: the geojson features to burn.
      size: the size of burned image.

    Returns:
      image: rasterized file of size with features burned.
    """

    # the value you want in the output raster where a shape exists
    burnval = 1
    shapes = ((geometry, burnval) for feature in features for geometry in feature_to_mercator(feature))

    bounds = mercantile.xy_bounds(tile)
    transform = from_bounds(*bounds, size, size)

    return rasterize(shapes, out_shape=(size, size), transform=transform)