Getting map tiles for polygon?

Here's solution in C#:

    private const double MinLatitude = -85.05112878;
    private const double MaxLatitude = 85.05112878;
    private const double MinLongitude = -180;
    private const double MaxLongitude = 180;

    public static Tuple<int, int> LatLongToTileXY(double latitude, double longitude, int z)
    {
        int tileX;
        int tileY;
        latitude = Clip(latitude, MinLatitude, MaxLatitude);
        longitude = Clip(longitude, MinLongitude, MaxLongitude);

        double x = (longitude + 180) / 360;
        double sinLatitude = Math.Sin(latitude * Math.PI / 180);
        double y = 0.5 - Math.Log((1 + sinLatitude) / (1 - sinLatitude)) / (4 * Math.PI);

        uint mapSize = MapSize(z);
        tileX = (int)Clip(x * mapSize + 0.5, 0, mapSize - 1) / 256;
        tileY = (int)Clip(y * mapSize + 0.5, 0, mapSize - 1) / 256;

        return Tuple.Create(tileX, tileY);
    }


/// <summary>
/// Determines the map width and height (in pixels) at a specified level
/// of detail.
/// </summary>
/// <param name="levelOfDetail">Level of detail, from 1 (lowest detail)
/// to 23 (highest detail).</param>
/// <returns>The map width and height in pixels.</returns>
public static uint MapSize(int levelOfDetail)
{
    return (uint)256 << levelOfDetail;
}

/// <summary>
/// Clips a number to the specified minimum and maximum values.
/// </summary>
/// <param name="n">The number to clip.</param>
/// <param name="minValue">Minimum allowable value.</param>
/// <param name="maxValue">Maximum allowable value.</param>
/// <returns>The clipped value.</returns>
private static double Clip(double n, double minValue, double maxValue)
{
    return Math.Min(Math.Max(n, minValue), maxValue);
}

private static Tuple<double, double> XY2Deg(int xtile, int ytile, int zoom)
{
    var n = Math.Pow(2.0, zoom);
    double lon_deg = xtile / n * 360.0 - 180.0;
    double lat_deg = (180 / Math.PI) * Math.Atan(Math.Sinh(Math.PI * (1 - 2 * ytile / n)));
    return Tuple.Create(lat_deg, lon_deg);
}


/// <summary>
/// Converts a pixel from pixel XY coordinates at a specified level of detail
/// into latitude/longitude WGS-84 coordinates (in degrees).
/// </summary>
/// <param name="pixelX">X coordinate of the point, in pixels.</param>
/// <param name="pixelY">Y coordinates of the point, in pixels.</param>
/// <param name="levelOfDetail">Level of detail, from 1 (lowest detail)
/// to 23 (highest detail).</param>
/// <param name="latitude">Output parameter receiving the latitude in degrees.</param>
/// <param name="longitude">Output parameter receiving the longitude in degrees.</param>
public static Tuple<double, double> TileXYToLatLong(int tileX, int tileY, int z)
{
    double latitude;
    double longitude;
    int pixelX = tileX * 256;
    int pixelY = tileY * 256;
    double mapSize = MapSize(z);
    double x = (Clip(pixelX, 0, mapSize - 1) / mapSize) - 0.5;
    double y = 0.5 - (Clip(pixelY, 0, mapSize - 1) / mapSize);

    latitude = 90 - 360 * Math.Atan(Math.Exp(-y * 2 * Math.PI)) / Math.PI;
    longitude = 360 * x;
    return Tuple.Create(latitude, longitude);
}




//get the range of tiles that intersect with the bounding box of the polygon    
private static Tuple<Tuple<int, int>, Tuple<int, int>> GetTileRange(DbGeography area, int zoom)
{
    //minimum bounding region (xm, ym, xmx, ymx)
    string dbGeography = area.AsText();
    var dbGeometry = DbGeometry.FromText(dbGeography);
    var bnds = dbGeometry.Envelope;

    double xm = bnds.PointAt(1).XCoordinate.Value;
    double xmx = bnds.PointAt(3).XCoordinate.Value;
    double ym = bnds.PointAt(1).YCoordinate.Value;
    double ymx = bnds.PointAt(3).YCoordinate.Value;

    var starting = LatLongToTileXY(ym, xm, zoom);
    var ending = LatLongToTileXY(ymx, xmx, zoom);

    var x_range = Tuple.Create(starting.Item1, ending.Item1);
    var y_range = Tuple.Create(ending.Item2, starting.Item2);

    return Tuple.Create(x_range, y_range);
}

private static Boolean DoesTileIntersects(int x, int y, int z, DbGeography area)
{
    ////Zoom tolerance; Below these zoom levels, only check if tile intersects with bounding box of polygon
    //if (z < 10)
    //    return true;
    DbGeography tile = GetTileASpolygon(x, y, z);
    bool intersects = area.Intersects(tile);
    return intersects;
}

//to get the tile as a polygon object
private static DbGeography GetTileASpolygon(int x, int y, int z)
{
    var nw = TileXYToLatLong(x, y, z);
    var se = TileXYToLatLong(x + 1, y + 1, z);

    return DbGeography.FromText( string.Format("POLYGON(({0} {1}, {0} {2}, {3} {2}, {3} {1}, {0} {1}))",
                 nw.Item2,
                 nw.Item1,
                 se.Item1,
                 se.Item2), 4326);
}


/// <summary>
/// Converts tile XY coordinates into a QuadKey at a specified level of detail.
/// </summary>
/// <param name="tileX">Tile X coordinate.</param>
/// <param name="tileY">Tile Y coordinate.</param>
/// <param name="levelOfDetail">Level of detail, from 1 (lowest detail)
/// to 23 (highest detail).</param>
/// <returns>A string containing the QuadKey.</returns>
public static string TileXYToQuadKey(int x, int y, int z)
{
    StringBuilder quadKey = new StringBuilder();
    for (int i = z; i > 0; i--)
    {
        char digit = '0';
        int mask = 1 << (i - 1);
        if ((x & mask) != 0)
        {
            digit++;
        }
        if ((y & mask) != 0)
        {
            digit++;
            digit++;
        }
        quadKey.Append(digit);
    }
    return quadKey.ToString();
}

/// <summary>
/// Converts a QuadKey into tile XY coordinates.
/// </summary>
/// <param name="quadKey">QuadKey of the tile.</param>
/// <param name="tileX">Output parameter receiving the tile X coordinate.</param>
/// <param name="tileY">Output parameter receiving the tile Y coordinate.</param>
/// <param name="levelOfDetail">Output parameter receiving the level of detail.</param>
public static Tuple<int, int, int> QuadKeyToTileXY(string quadKey)
{
    int tileX;
    int tileY;
    int zoom;

    tileX = tileY = 0;
    zoom = quadKey.Length;
    for (int i = zoom; i > 0; i--)
    {
        int mask = 1 << (i - 1);
        switch (quadKey[zoom - i])
        {
            case '0':
                break;

            case '1':
                tileX |= mask;
                break;

            case '2':
                tileY |= mask;
                break;

            case '3':
                tileX |= mask;
                tileY |= mask;
                break;

            default:
                throw new ArgumentException("Invalid QuadKey digit sequence.");
        }
    }
    return Tuple.Create(tileX, tileY, zoom);
}


// entry point
public static List<string> GetTiles(DbGeography area)
{
    var tiles = new List<string>();
    for (int z = 1; z <= 16; z++)
    {
        var ranges = FeatureHelper.GetTileRange(area, z);
        var x_range = ranges.Item1;
        var y_range = ranges.Item2;


        for (int y = y_range.Item1; y < y_range.Item2 + 1; y++)
        {
            for (int x = x_range.Item1; x < x_range.Item2 + 1; x++)
            {
                if (FeatureHelper.DoesTileIntersects(x, y, z, area))
                    tiles.Add(TileXYToQuadKey(x, y, z));
            }
        }
    }
    return tiles; }