Posts tagged ‘Tile Layer’

November 28, 2012

Masking out particular areas of interest on Bing Maps

A recent question on the MSDN Bing Maps forums asked if it was possible to display a Bing Map in which only a particular country was visible, painting out everything else. There’s certainly nothing built into the API to do this – the tiles from which maps are constructed are pre-rendered images and know nothing about country outlines displayed on them.

However, if you have coordinate data representing the area of interest you want to display – be that a country or any other area of interest – I don’t think it should be too hard to achieve the effect desired using a dynamic tilelayer. The principle is exactly the same as if you were creating a raster tilehandler that dynamically drew polygons representing an area of interest. The only difference is that, instead of shading the polygon defined by your data, you invert the graphics to shade everything else, leaving the polygon of interest transparent.

If you’re using .NET to draw your tiles, what you’d do is create a graphics object representing the tile as usual, but instead of calling FillRegion() on the polygon drawn on the tile, you’d create a region covering the whole tile and then use Region.Exclude() to “remove” the area covered by the polygon, before filling the rest of the region.

As an example, I grabbed the coordinates for the island of Corsica from GADM, and put them into an array, as:

PointF[] selectedArea = {
  new PointF(9.17986f, 41.36541f),
  new PointF(9.22097f, 41.36625f),
  new PointF(9.2198...
  ...
}

If drawn using a conventional tile rasteriser, this would look like:

image

However, the addition of the following two lines inverts the region to be shaded:

// Define a new region that spans the whole tile
Region region = new Region(new Rectangle(0, 0, 256, 256));

// Exclude the polygon from the region
region.Exclude(polygonPath);

Re-rendering the map (and changing the fill colour to blue) now gives:

image

And, if you tweak the shading options a bit and draw a line around the region you can create something a bit more aesthetically appealing. I’ve only added the masking layer with opacity 0.7 so that you can see the other countries slightly in the background, but you could make the whole rest of the map opaque if you wanted.

So here’s a Bing Map that masks everything other than Corsica. Although the screenshot below is static, this solution uses a normal tile handler so the map retains all the normal panning/zooming functionality, and only Corsica will remain in clear view:

image

Here’s the code for the HTML page shown above:


<!DOCTYPE HTML PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html  xmlns="http://www.w3.org/1999/xhtml">
<head>
<title></title>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<script type="text/javascript" src="http://ecn.dev.virtualearth.net/mapcontrol/mapcontrol.ashx?v=7.0"></script>
<script type="text/javascript">
function GetMap() {

var map = new Microsoft.Maps.Map(document.getElementById("mapDiv"),
{ credentials: "ENTERYOURBINGMAPSKEYHERE",
center: new Microsoft.Maps.Location(42, 9),
mapTypeId: Microsoft.Maps.MapTypeId.road,
zoom: 8
});

// Create the tile layer source
var tileSource = new Microsoft.Maps.TileSource(
{  uriConstructor: 'TileHandler.ashx?q={quadkey}',
width: 256,
height: 256
});

// Construct the layer using the tile source
var tilelayer = new Microsoft.Maps.TileLayer({ mercator: tileSource, opacity: 0.8 });

// Push the tile layer to the map
map.entities.push(tilelayer);
}
</script>
</head>
<body onload="GetMap();">
<div id='mapDiv' style="position:relative; width:800px; height:600px;"></div>
</body>
</html>

And here’s the code for the TileHandler.ashx:


<%@ WebHandler Language="C#" Debug="true" %>

using System;
using System.Web;
using System.Drawing;
using System.Drawing.Drawing2D;
using Microsoft.SqlServer.Types;
using System.Data.SqlClient;

public class Handler : IHttpHandler
{

public void ProcessRequest(HttpContext context)
{

// Retrieve the requested quadkey
string quadKey = context.Request.QueryString["q"];
int zoomLevel;

// Convert quadkey to tile coordinates
int tileX, tileY;
QuadKeyToTileXY(quadKey, out tileX, out tileY, out zoomLevel);

// Calculate the pixel coordinates of the top left corner of this tile
int TLpixelX, TLpixelY;
TileXYToPixelXY(tileX, tileY, out TLpixelX, out TLpixelY);

// Use the map image as a graphics object
Bitmap mapImage = new System.Drawing.Bitmap(256, 256);
System.Drawing.Graphics graphics = System.Drawing.Graphics.FromImage(mapImage);

// Define Long/Lats of an area of interest (Corsica, in this example)
PointF[] selectedArea = {new PointF(9.17986f,41.36541f),new PointF(9.22097f,41.36625f),new PointF(9.21986f,41.37486f),new PointF(9.25958f,41.41347f),new PointF(9.26486f,41.42847f),new PointF(9.21986f,41.40514f),new PointF(9.22902f,41.42624f),new PointF(9.22152f,41.43986f),new PointF(9.21319f,41.44041f),new PointF(9.24347f,41.44652f),new PointF(9.26652f,41.46625f),new PointF(9.27791f,41.46374f),new PointF(9.26625f,41.47041f),new PointF(9.28430f,41.47735f),new PointF(9.28819f,41.49263f),new PointF(9.28347f,41.50069f),new PointF(9.26597f,41.49875f),new PointF(9.28486f,41.51652f),new PointF(9.27236f,41.52319f),new PointF(9.27347f,41.53125f),new PointF(9.28708f,41.52875f),new PointF(9.30291f,41.54569f),new PointF(9.30902f,41.54264f),new PointF(9.30874f,41.54958f),new PointF(9.34402f,41.55902f),new PointF(9.35069f,41.56652f),new PointF(9.34708f,41.57513f),new PointF(9.36763f,41.58930f),new PointF(9.36847f,41.59791f),new PointF(9.34069f,41.59208f),new PointF(9.33736f,41.60152f),new PointF(9.31958f,41.60514f),new PointF(9.29263f,41.58152f),new PointF(9.28152f,41.59791f),new PointF(9.28736f,41.60874f),new PointF(9.30736f,41.61569f),new PointF(9.30763f,41.62847f),new PointF(9.32597f,41.61875f),new PointF(9.32069f,41.61513f),new PointF(9.33402f,41.62152f),new PointF(9.35569f,41.61708f),new PointF(9.34847f,41.63791f),new PointF(9.35646f,41.64070f),new PointF(9.35569f,41.64097f),new PointF(9.36680f,41.64430f),new PointF(9.35646f,41.64070f),new PointF(9.37069f,41.63569f),new PointF(9.37735f,41.64847f),new PointF(9.38624f,41.64541f),new PointF(9.38180f,41.65152f),new PointF(9.38930f,41.65652f),new PointF(9.38125f,41.66402f),new PointF(9.39680f,41.66902f),new PointF(9.37624f,41.67013f),new PointF(9.37236f,41.67930f),new PointF(9.40124f,41.69235f),new PointF(9.39930f,41.70291f),new PointF(9.40819f,41.71374f),new PointF(9.40180f,41.71708f),new PointF(9.40847f,41.76653f),new PointF(9.39569f,41.77485f),new PointF(9.40041f,41.78902f),new PointF(9.39458f,41.79902f),new PointF(9.40652f,41.81708f),new PointF(9.40541f,41.85791f),new PointF(9.39735f,41.86125f),new PointF(9.39680f,41.87541f),new PointF(9.41597f,41.93097f),new PointF(9.40041f,41.94625f),new PointF(9.40708f,41.95152f),new PointF(9.39930f,41.95180f),new PointF(9.41125f,41.95375f),new PointF(9.41236f,41.93930f),new PointF(9.42069f,41.96736f),new PointF(9.54902f,42.10319f),new PointF(9.55486f,42.12430f),new PointF(9.55958f,42.19652f),new PointF(9.55208f,42.23208f),new PointF(9.56041f,42.28263f),new PointF(9.53236f,42.37041f),new PointF(9.54319f,42.42680f),new PointF(9.54180f,42.45902f),new PointF(9.52930f,42.49124f),new PointF(9.53319f,42.54597f),new PointF(9.51097f,42.58930f),new PointF(9.44736f,42.65652f),new PointF(9.44680f,42.68652f),new PointF(9.45819f,42.70264f),new PointF(9.46652f,42.76097f),new PointF(9.49180f,42.79374f),new PointF(9.49319f,42.80736f),new PointF(9.48041f,42.83708f),new PointF(9.48625f,42.85125f),new PointF(9.48208f,42.86680f),new PointF(9.47291f,42.87486f),new PointF(9.47819f,42.88208f),new PointF(9.46930f,42.93624f),new PointF(9.45208f,42.96152f),new PointF(9.46402f,42.98625f),new PointF(9.42013f,43.01208f),new PointF(9.41291f,43.00541f),new PointF(9.35736f,43.00708f),new PointF(9.34013f,42.99430f),new PointF(9.35041f,42.96791f),new PointF(9.34597f,42.95625f),new PointF(9.35875f,42.94541f),new PointF(9.34930f,42.92958f),new PointF(9.35958f,42.92375f),new PointF(9.32180f,42.89902f),new PointF(9.32958f,42.87208f),new PointF(9.34013f,42.86541f),new PointF(9.30958f,42.83319f),new PointF(9.34263f,42.79708f),new PointF(9.33708f,42.76458f),new PointF(9.34513f,42.73513f),new PointF(9.32208f,42.71652f),new PointF(9.32236f,42.69652f),new PointF(9.29791f,42.68152f),new PointF(9.30041f,42.67597f),new PointF(9.28541f,42.67541f),new PointF(9.27069f,42.70041f),new PointF(9.25430f,42.70513f),new PointF(9.25680f,42.71875f),new PointF(9.22902f,42.71819f),new PointF(9.23236f,42.72347f),new PointF(9.22208f,42.73569f),new PointF(9.19986f,42.72513f),new PointF(9.16763f,42.73708f),new PointF(9.13569f,42.72874f),new PointF(9.12347f,42.73263f),new PointF(9.09847f,42.71513f),new PointF(9.08513f,42.71569f),new PointF(9.08458f,42.70486f),new PointF(9.07013f,42.69291f),new PointF(9.05875f,42.69652f),new PointF(9.05986f,42.68291f),new PointF(9.05319f,42.68152f),new PointF(9.06013f,42.66125f),new PointF(9.02902f,42.65375f),new PointF(9.01486f,42.64014f),new PointF(9.00236f,42.64430f),new PointF(8.94652f,42.63319f),new PointF(8.93319f,42.64680f),new PointF(8.93625f,42.64152f),new PointF(8.91263f,42.62930f),new PointF(8.87958f,42.62930f),new PointF(8.86902f,42.60791f),new PointF(8.85041f,42.61236f),new PointF(8.83125f,42.60152f),new PointF(8.82708f,42.60680f),new PointF(8.79847f,42.60263f),new PointF(8.81041f,42.59486f),new PointF(8.79847f,42.58291f),new PointF(8.80513f,42.57208f),new PointF(8.78375f,42.55680f),new PointF(8.75986f,42.55847f),new PointF(8.76319f,42.56902f),new PointF(8.75569f,42.57180f),new PointF(8.72763f,42.56180f),new PointF(8.72430f,42.58458f),new PointF(8.70625f,42.57013f),new PointF(8.71791f,42.57069f),new PointF(8.71569f,42.56013f),new PointF(8.72264f,42.55569f),new PointF(8.71597f,42.55097f),new PointF(8.71874f,42.54013f),new PointF(8.70986f,42.53625f),new PointF(8.71958f,42.52541f),new PointF(8.69847f,42.52736f),new PointF(8.68402f,42.51514f),new PointF(8.66847f,42.51763f),new PointF(8.66430f,42.49319f),new PointF(8.64874f,42.47986f),new PointF(8.65041f,42.47319f),new PointF(8.67652f,42.47625f),new PointF(8.67986f,42.46958f),new PointF(8.66402f,42.45763f),new PointF(8.66764f,42.44569f),new PointF(8.64763f,42.44319f),new PointF(8.66319f,42.43402f),new PointF(8.65597f,42.41708f),new PointF(8.64569f,42.41347f),new PointF(8.62319f,42.42263f),new PointF(8.60791f,42.41764f),new PointF(8.60069f,42.40874f),new PointF(8.61041f,42.40486f),new PointF(8.60152f,42.39902f),new PointF(8.60847f,42.38625f),new PointF(8.57625f,42.38402f),new PointF(8.57819f,42.37486f),new PointF(8.56874f,42.37041f),new PointF(8.55402f,42.37180f),new PointF(8.54736f,42.38124f),new PointF(8.54541f,42.36986f),new PointF(8.53430f,42.37236f),new PointF(8.53847f,42.36458f),new PointF(8.55541f,42.36458f),new PointF(8.55097f,42.34513f),new PointF(8.55930f,42.34319f),new PointF(8.55291f,42.33291f),new PointF(8.57319f,42.33708f),new PointF(8.59152f,42.35291f),new PointF(8.61513f,42.34958f),new PointF(8.61763f,42.34124f),new PointF(8.62847f,42.34013f),new PointF(8.62736f,42.33180f),new PointF(8.60124f,42.32208f),new PointF(8.60124f,42.30875f),new PointF(8.62625f,42.31263f),new PointF(8.63958f,42.29902f),new PointF(8.66097f,42.30208f),new PointF(8.67458f,42.29375f),new PointF(8.67374f,42.28319f),new PointF(8.68874f,42.28013f),new PointF(8.69180f,42.26930f),new PointF(8.63513f,42.25291f),new PointF(8.60986f,42.25402f),new PointF(8.55958f,42.23569f),new PointF(8.53847f,42.23736f),new PointF(8.55125f,42.22680f),new PointF(8.57319f,42.22819f),new PointF(8.56735f,42.21874f),new PointF(8.57708f,42.21402f),new PointF(8.56819f,42.21347f),new PointF(8.56569f,42.20486f),new PointF(8.58180f,42.20680f),new PointF(8.57736f,42.18791f),new PointF(8.58624f,42.18097f),new PointF(8.56069f,42.17180f),new PointF(8.59375f,42.16903f),new PointF(8.55736f,42.14541f),new PointF(8.59097f,42.14902f),new PointF(8.59375f,42.14347f),new PointF(8.57930f,42.12791f),new PointF(8.61569f,42.13347f),new PointF(8.62291f,42.12208f),new PointF(8.65125f,42.11875f),new PointF(8.65819f,42.10069f),new PointF(8.68097f,42.10486f),new PointF(8.68986f,42.11541f),new PointF(8.70041f,42.11291f),new PointF(8.70986f,42.09791f),new PointF(8.70180f,42.08513f),new PointF(8.71069f,42.08708f),new PointF(8.72208f,42.07208f),new PointF(8.71569f,42.06347f),new PointF(8.73597f,42.06597f),new PointF(8.74847f,42.04958f),new PointF(8.74263f,42.04124f),new PointF(8.72486f,42.04319f),new PointF(8.72458f,42.03402f),new PointF(8.68569f,42.02763f),new PointF(8.65541f,42.01013f),new PointF(8.65874f,41.99152f),new PointF(8.67652f,41.98930f),new PointF(8.64958f,41.96930f),new PointF(8.61430f,41.97180f),new PointF(8.59208f,41.96402f),new PointF(8.59847f,41.95208f),new PointF(8.60736f,41.95208f),new PointF(8.60986f,41.93903f),new PointF(8.62347f,41.93569f),new PointF(8.60708f,41.89374f),new PointF(8.64125f,41.90986f),new PointF(8.71847f,41.90736f),new PointF(8.74124f,41.91541f),new PointF(8.74791f,41.93402f),new PointF(8.76458f,41.92180f),new PointF(8.78319f,41.92541f),new PointF(8.80430f,41.89569f),new PointF(8.79097f,41.88430f),new PointF(8.77847f,41.88430f),new PointF(8.78513f,41.87958f),new PointF(8.78069f,41.87319f),new PointF(8.79125f,41.86819f),new PointF(8.78902f,41.85097f),new PointF(8.74875f,41.84541f),new PointF(8.78402f,41.83263f),new PointF(8.77236f,41.82291f),new PointF(8.77430f,41.81236f),new PointF(8.74764f,41.81124f),new PointF(8.73736f,41.79680f),new PointF(8.71180f,41.80263f),new PointF(8.71013f,41.79458f),new PointF(8.73235f,41.77958f),new PointF(8.69569f,41.75152f),new PointF(8.66124f,41.75041f),new PointF(8.65930f,41.73847f),new PointF(8.70680f,41.73874f),new PointF(8.70236f,41.72569f),new PointF(8.70930f,41.72180f),new PointF(8.77597f,41.74125f),new PointF(8.78430f,41.73180f),new PointF(8.77097f,41.71347f),new PointF(8.78736f,41.70875f),new PointF(8.78125f,41.69930f),new PointF(8.81235f,41.71458f),new PointF(8.84319f,41.69652f),new PointF(8.91458f,41.69097f),new PointF(8.91375f,41.67986f),new PointF(8.88180f,41.67096f),new PointF(8.87208f,41.64625f),new PointF(8.85124f,41.64819f),new PointF(8.82041f,41.62958f),new PointF(8.80569f,41.64208f),new PointF(8.79402f,41.63263f),new PointF(8.79513f,41.62347f),new PointF(8.78513f,41.61791f),new PointF(8.79125f,41.61069f),new PointF(8.78652f,41.59708f),new PointF(8.77374f,41.59152f),new PointF(8.80402f,41.57541f),new PointF(8.78291f,41.56708f),new PointF(8.78736f,41.55680f),new PointF(8.81069f,41.55652f),new PointF(8.82097f,41.54375f),new PointF(8.84458f,41.54680f),new PointF(8.84541f,41.54013f),new PointF(8.85347f,41.54680f),new PointF(8.85597f,41.53513f),new PointF(8.84208f,41.52597f),new PointF(8.84236f,41.51819f),new PointF(8.87652f,41.52485f),new PointF(8.89069f,41.51597f),new PointF(8.88319f,41.50485f),new PointF(8.91930f,41.50541f),new PointF(8.92041f,41.48847f),new PointF(8.93264f,41.49624f),new PointF(8.93847f,41.48875f),new PointF(8.96041f,41.49180f),new PointF(8.96763f,41.48875f),new PointF(8.96513f,41.48208f),new PointF(8.97791f,41.48291f),new PointF(8.98097f,41.47347f),new PointF(8.99402f,41.48764f),new PointF(9.00263f,41.47486f),new PointF(9.01930f,41.47986f),new PointF(9.01402f,41.47208f),new PointF(9.02402f,41.46125f),new PointF(9.03791f,41.47152f),new PointF(9.04041f,41.45624f),new PointF(9.07902f,41.47652f),new PointF(9.06541f,41.45430f),new PointF(9.07347f,41.44374f),new PointF(9.09430f,41.44930f),new PointF(9.10847f,41.44013f),new PointF(9.12458f,41.44374f),new PointF(9.10541f,41.42847f),new PointF(9.10958f,41.42069f),new PointF(9.09458f,41.41236f),new PointF(9.09347f,41.39375f),new PointF(9.13347f,41.39930f),new PointF(9.12652f,41.39430f),new PointF(9.16958f,41.38375f),new PointF(9.17986f,41.36541f)};

// Convert shape from Lat/Lngs to pixel coordinates relative to this tile
Point[] points = new Point[selectedArea.Length];
for (int k = 0; k < selectedArea.Length; k++)
{
int x = LongitudeToXAtZoom(selectedArea[k].X, zoomLevel);
int y = LatitudeToYAtZoom(selectedArea[k].Y, zoomLevel);
points[k] = new Point(x-TLpixelX, y-TLpixelY);
}

// Create a polygon from these points
GraphicsPath polygonPath = new GraphicsPath();
polygonPath.AddPolygon(points);

// Define a new region that spans this tile
Region region = new Region(new Rectangle(0, 0, 256, 256));

// Exclude the polygon from the region
region.Exclude(polygonPath);

// Now fill everything but the polygon
graphics.FillRegion(Brushes.Black, region);

graphics.SmoothingMode = SmoothingMode.AntiAlias;
graphics.DrawLines(
new System.Drawing.Pen(System.Drawing.Color.FromArgb(255, 255, 255, 255), 4),
points);

// Send the tile image back to the client
context.Response.ContentType = "image/png";
System.IO.MemoryStream memStream = new System.IO.MemoryStream();
mapImage.Save(memStream, System.Drawing.Imaging.ImageFormat.Png);
memStream.WriteTo(context.Response.OutputStream);

}

///
/// COMMON UTILITY FUNCTIONS
/// See http://msdn.microsoft.com/en-us/library/bb259689.aspx
///

private const double earthRadius = 6378137; //The radius of the earth - should never change!
private const double earthCircum = earthRadius * 2.0 * Math.PI; //calulated circumference of the earth
private const double earthHalfCirc = earthCircum / 2; //calulated half circumference of the earth
private const int pixelsPerTile = 256;

public static void QuadKeyToTileXY(string quadKey, out int tileX, out int tileY, out int levelOfDetail)
{
tileX = tileY = 0;
levelOfDetail = quadKey.Length;
for (int i = levelOfDetail; i > 0; i--)
{
int mask = 1 << (i - 1);
switch (quadKey[levelOfDetail - 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.");
}
}
}

public static void TileXYToPixelXY(int tileX, int tileY, out int pixelX, out int pixelY)
{
pixelX = tileX * 256;
pixelY = tileY * 256;
}

public static int LatitudeToYAtZoom(double lat, int zoom)
{
int y;
double arc = earthCircum / ((1 << zoom) * pixelsPerTile);
double sinLat = Math.Sin(DegToRad(lat));
double metersY = earthRadius / 2 * Math.Log((1 + sinLat) / (1 - sinLat));
y = (int)Math.Round((earthHalfCirc - metersY) / arc);
return y;
}

public static int LongitudeToXAtZoom(double lon, int zoom)
{
int x;
double arc = earthCircum / ((1 << zoom) * pixelsPerTile);
double metersX = earthRadius * DegToRad(lon);
x = (int)Math.Round((earthHalfCirc + metersX) / arc);
return x;
}

private static double DegToRad(double d)
{
return d * Math.PI / 180.0;
}

public bool IsReusable
{
get
{
return false;
}
}

}
July 8, 2011

Bing Maps Ready Reckoner

Compared to more fully-featured and complex GIS systems, the maths behind Bing Maps is relatively straightforward:

  • Bing Maps are always projected using the same Spherical Mercator projection (EPSG: 3857).
  • Bing Maps are always displayed at one of a fixed number of zoom levels. At zoom level 1, the entire world fits onto a 2 x 2 grid of 256px x 256px tiles. Then, at every subsequent zoom level, the resolution in both the vertical and horizontal dimensions doubles, up to a maximum of zoom level 19.
  • Vector coordinate data is always measured using latitude and longitude coordinates from the WGS84 reference system (EPSG:4326).
  • Raster data is always supplied as a set of 256px x 256px tiles, referenced according to the quadkey tile numbering system.

The advantage of operating within these fixed parameters is that it’s possible to create some rules of thumb that can be applied fairly consistently when developing any Bing map application. The following “Ready Reckoner” table(s) gives some figures I find useful when creating datasets for use on Bing maps. Since they use the same projection, these same figures can also be used as guides for Google Maps or Openlayers, say.

Firstly, some constants:

  • The full extents of every Bing Map range from –180° to +180° longitude, and from –85.05112878° to +85.05112878° latitude. This makes the projected image appear perfectly square at any zoom level
  • The straight line distance right across the middle of the map (i.e. the circumference of the Earth at the equator) is 40,075,018 metres.
  • The width/height of the map at any zoom level, z,  can be calculated as 256 * (2z)

Part 1: Map Scale and Resolution

The following table gives the total map dimensions, scale, and resolution at different zoom levels:

Zoom Level Width/Height of map (pixels) Scale (m/pixel) Scale
(°
longitude/pixel)
Scale
(°
latitude/pixel)
Coordinate Resolution
1 512 78271.52 0.703125 0.332230972 0.1
2 1024 39135.76 0.3515625 0.166115486 0.1
3 2048 19567.88 0.17578125 0.083057743 0.01
4 4096 9783.94 0.087890625 0.041528871 0.01
5 8192 4891.97 0.043945313 0.020764436 0.01
6 16384 2445.98 0.021972656 0.010382218 0.01
7 32768 1222.99 0.010986328 0.005191109 0.001
8 65536 611.5 0.005493164 0.002595554 0.001
9 131072 305.75 0.002746582 0.001297777 0.001
10 262144 152.87 0.001373291 0.000648889 0.0001
11 524288 76.44 0.000686646 0.000324444 0.0001
12 1048576 38.22 0.000343323 0.000162222 0.0001
13 2097152 19.11 0.000171661 8.11111E-05 0.00001
14 4194304 9.55 8.58307E-05 4.05555E-05 0.00001
15 8388608 4.78 4.29153E-05 2.02778E-05 0.00001
16 16777216 2.39 2.14577E-05 1.01389E-05 0.00001
17 33554432 1.19 1.07288E-05 5.06944E-06 0.000001
18 67108864 0.6 5.36442E-06 2.53472E-06 0.000001
19 134217728 0.3 2.68221E-06 1.26736E-06 0.000001

Notes

Scale: The various columns of scale are all measured at the equator – the scale on any projected Mercator map varies with latitude. Notice that, given the projected map image is square, yet the range of possible longitude values is roughly double that of the range of latitude values, the map is approximately 2x as precise in the latitude/y axis than it is in the longitude/x axis. Or, to put it another way, you have to move your mouse twice as far up the screen to cover one degree of latitude than you have to move it across the screen to cover one degree of longitude.

Coordinate resolution: The “coordinate resolution” column gives an indication of how many decimal places any coordinate value can be rounded to and still maintain the same visual acuity at a given zoom level. In other words, at zoom level 3, you may as well round any coordinate values to 2 decimal places – any coordinate stated with accuracy beyond that will likely be resolved to the same pixel location anyway.

Part 2: Raster Tiling

When creating raster tilesets (using Safe FME, GDAL2Tiles, or similar), it’s handy to have an idea of roughly how many tiles are required to cover a given area, and how long it’s likely to take to render those tiles! In the tables that follow, I’ve given information on the tiles required to cover three different-sized areas at different zoom levels.

Entire Map

The following table gives figures relating to a complete tileset covering the entire map:

Zoom Level Number of tiles required Time To Render
(at 1 tile/sec)
Uncompressed Filesize
(at 20kb/tile)
Compressed Filesize
(at 6kb/tile)
1 4 4 seconds 80 Kb 24 Kb
2 16 16 seconds 320 Kb 96 Kb
3 64 1 minute 1.25 Mb 384 Kb
4 256 4 minutes 5 Mb 1.5 Mb
5 1,024 20 minutes 20 Mb 6 Mb
6 4,096 1 hour 80 Mb 24 Mb
7 16,384 4 hours 320 Mb 96 Mb
8 65,536 20 hours 1.25 Gb 384 Mb
9 262,144 3 days 5 Gb 1.5 Gb
10 1,048,576 12 days 20 Gb 6 Gb
11 4,194,304 7 weeks 80 Gb 24 Gb
12 16,777,216 27 weeks 320 Gb 96 Gb
13 67,108,864 2 years 1.25 Tb 384 Gb
14 268,435,456 8 years 5 Tb 1.5 Tb
15 1,073,741,824 Give up. 20 Tb 6 Tb
16 4,294,967,296 No, really. 80 Tb 24 Tb
17 17,179,869,184 320 Tb 96 Tb
18 68,719,476,736 1.25 Pb 384 Tb
19 274,877,906,944 5 Pb 1.5 Pb

10km x 10km Square

Now, just to tile just an arbitrary 10km x 10km area:

Zoom Level Number of Tiles Required Time To Render Uncompressed Filesize
1 1 1 second 20 Kb
2 1 1 second 20 Kb
3 1 1 second 20 Kb
4 1 1 second 20 Kb
5 1 1 second 20 Kb
6 1 1 second 20 Kb
7 1 1 second 20 Kb
8 1 1 second 20 Kb
9 1 1 second 20 Kb
10 1 1 second 20 Kb
11 2 1 second 41.6 kB
12 4 2 seconds 124 kb
13 16 8 seconds 320 kb
14 64 30 seconds 900 Kb
15 225 2 minutes 3 Mb
16 840 10 minutes * 9.23 Mb *
17 3000 40 minutes 35 Mb
18 12000 3 hours 100 Mb
19 48000 12 hours 500 Mb

100 km x 100km Square

And, to cover a 100km x 100km square (not the same as just multiplying the last table by 100! – see notes):

Zoom Level Number of Tiles Time To Render Uncompressed Filesize
1 1 1 second 20 Kb
2 1 1 second 20 Kb
3 1 1 second 20 Kb
4 1 1 second 20 Kb
5 1 1 second 20 Kb
6 1 1 second 20 Kb
7 2 1 second 40 Kb
8 4 3 seconds 160 Kb
9 9 15 seconds 600 Kb
10 25 1 minute 2.2 Mb
11 90 2 minutes 7.8 Mb
12 324 5 minutes 25 Mb
13 1,190 20 minutes 78 Mb
14 4,624 1 hour 240 Mb
15 18,360 4.5 hours 700 Mb
16 72,900 15 hours * 1.9 Gb *
17 291,600 3 days 5.6 Gb
18 1,166,400 10 days 22 Gb
19 4,665,600 1.5 months 100 Gb

Notes

Number of Tiles Required: A 10km x 10km area is sufficient to cover the extent of a small to medium-sized town/city. If you’re familiar with the National Grid of Great Britain, then 10km x 10km corresponds to an OS Grid Reference square of two letters and two numbers, e.g. TG20, shown below:

image

A 100km x 100km square is sufficient to cover an expansive metropolis and its surrounding area, or a British county or two, say (though is still much smaller than the average US or Australian state). A 100km x 100km square in the Great Britain OS Grid reference system is denoted by just a two letter code, e.g. TQ shown below:

image

 

The number of tiles required to tile the entire world in the table above is calculated from mathematical properties of the map. The number of tiles for 10km x 10km and 100km x 100km are based on real values that I’ve collected from tiling Ordnance Survey grid data, which, as shown in the images above, is not axis-aligned to the Bing Maps projection. This explains why the number of tiles required does not always scale consistently with each zoom level. The values with asterisks (zoom level 16) in the tables are the actual figures I got just this morning from rendering OS VectorMap tiles for the exact two areas shown above – TG20 and TQ – on a single core of my laptop.

Time To Render: This figure varies wildly depending on your particular system configuration, whether your source data is an existing raster file or whether it’s being generated dynamically (e.g. from a database query), and many other factors. The column here gives timings at a rate of 1 tile per second simply as a convenient figure from which to perform further approximate calculations. If your system renders 6 tiles per second, you can divide these times by 6 ;)

Note also that the time to render a set of tiles at zoom levels 14 – 16, say, is not necessarily the sum of the times taken to render at each of those zoom levels independently. Many tile rendering engines will render the tiles at the deepest zoom in full, and then aggregate those tiles together in groups of 4 to render the subsequent higher layers more quickly than rendering them directly from source.

Filesize: I find that, by default, the PNG images output by applications such as Safe FME and GDAL are not compressed, and a 256 x 256px tile may, on average be around 20kB in size, which is the baseline used in the “uncompressed” column. Running these tiles through one or more PNG compressions programs can reduce these filesizes (without losing quality) to around 6kB each, which is the baseline used in the “compressed” column. Also, the total filesize of a tileset varies a lot with the density distribution of your data – at close zoom levels, sparse datasets will have a lot of empty tiles, meaning that filesize does not increase linearly with number of tiles generated.

July 6, 2011

Converting TMS Tile Coordinates to Google/Bing/OSM Tile Coordinates

In addition to referencing tiles by quadkey, Bing Maps also refers to tiles by x, y, and z coordinates. This is the way you reference tiles when overriding the GetUri() method of the Silverlight Microsoft.Maps.MapControl.TileSource class, for example, and the x, y, z properties of a tile are passed to the uriConstructor function specified in the options for an AJAX v7 Microsoft.Maps.TileSource.

In the x, y, z tile identification system:

  • z is the zoom level of the grid
  • x and y are the column and row number at which this tile should be placed, measured from an origin at the top left of the map.

Thus, at zoom level 1, you can describe the four tiles required to make a complete Bing Map using x, y, z coordinates as follows:

image

Bing Maps is not entirely unusual in this respect – this exact same referencing system is used by Google Maps, Open Street Maps, ESRI, and many others. You’d therefore be forgiven for thinking that it was some sort of standard.

In fact, this tile numbering sounds a lot like (and is often incorrectly described as) a  TMS index. The Tile Map Service (TMS) Specification defined by the Open Source Geospatial Foundation is a standard for serving map tiles that, like the system described above, places tiles on a grid and refers to their position using x, y, and z coordinates, where z is the zoom level, and x and y refer to column and row positions. However, according to the TMS specifications, “The x-coordinate of the tile numbers increases with the x-coordinate of the spatial reference system, and the y-coordinate of the tile numbers also increases with the y-coordinate of the spatial reference system.”. In other words, tile (0, 0), at any zoom level should always be placed at the bottom left of the map, not the top left.

Using the TMS system, the tile indexes at zoom level 1 become as follows:

image

The problem is that, since these systems are identical in almost every other respect, people tend to assume that the system used by Google/Bing/OSM et al. is the standard and systems that output tiles using it sometimes mistakenly refer to it as the TMS format. Then,  every now and again, you come across a piece of software or service that genuinely outputs tiles numbered according to the TMS standard and, if you don’t correct the tile origin, your tiles all end upside-down on the wrong side of the world :(

Fortunately, it’s a very simple correction to make. At any given map zoom level, zoom, you can invert the y index of a tile from TMS to Google/Bing as follows:

var ymax = 1 << zoom;
var y = ymax - y - 1;

Assuming that you have a set of TMS tiles stored in a subdirectory structure that follows the pattern /z/y/x.png (as described in the Tile Resources section of the OSGeo TMS standard), then here’s a full example showing how to add a tilelayer of TMS tiles to Bing Maps v7:

<!DOCTYPE HTML PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html  xmlns="http://www.w3.org/1999/xhtml">
<head>
<title></title>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<script type="text/javascript" src="http://ecn.dev.virtualearth.net/mapcontrol/mapcontrol.ashx?v=7.0"></script>
<script type="text/javascript">
  var map = null;
  function GetMap() {
    // Create a basic map
    map = new Microsoft.Maps.Map(document.getElementById("mapDiv"),
      { credentials: "ENTERYOURBINGMAPSKEY",
        center: new Microsoft.Maps.Location(52.6, 1.26),
        zoom: 12
      });

    // Create the tile source
    var tileSource = new Microsoft.Maps.TileSource({ uriConstructor: getTMSTilePath });

    // Construct the layer using the tile source
    var tilelayer = new Microsoft.Maps.TileLayer({ mercator: tileSource, opacity: 1 });

    // Push the tile layer to the map
    map.entities.push(tilelayer);
  }

  function getTMSTilePath(tile) {

    var x = tile.x;
    var z = tile.levelOfDetail;
    // Invert tile y origin from top to bottom of map
    var ymax = 1 << z;
    var y = ymax - tile.y - 1;

    return location.href.substring(0, location.href.lastIndexOf('/')) + "/" + z + "/" + x + "/" + y + ".png";
  }
</script>
</head>
<body onload="GetMap();">
  <div id='mapDiv' style="position:relative; width:1024px; height:768px;"></div>
</body>
</html>

As a final note, just to add to the confusion, the Web Map Tile Service (WMTS) recently published by the OGC provides an alternative standard to TMS, but, like Google/Bing, uses an origin at the top-left hand corner.  For a discussion of the differences between TMS and WMTS, see http://2010.foss4g.org/presentations/3653.pdf

June 16, 2011

Tile Rendering with Mapnik (Part 1 of n, where n is large)

In my last blog post, I explained a little about the process by which you can overlay raster data on Bing Maps, by geo-referencing the source data (if necessary), projecting into the appropriate spherical Mercator projection, then cutting the resulting image into 256px x 256px tiles named according to the quadkey tile numbering system.

The application I used in the last example to do this was Microsoft Mapcruncher. Mapcruncher has got a lot of benefits:

  • It’s free
  • It’s a windows application with no special dependencies and easy to install
  • It’s got a nice GUI that’s very easy to use
  • It will perform geo-referencing/warping/tile-cutting for you as part of a single process

Mapcruncher is great if you have a relatively small raster image that you want to integrate into Bing Maps / Google Maps as part of a one-off process. However, there are many occasions when you might need a little more than that. As part of a current project, for example, I need to create a tileset that covers an area of approximately 60km x 50km (not that small), based on data held in SQL Server 2008 (not raster), that will be updated approximately every month (not one-off).

So I set out to evaluate alternatives, and the first I considered was mapnik.

Installing Mapnik

Mapnik is an open-source mapping toolkit. It’s what openstreetmap uses to render the tiles used in their base map imagery. They’ve got over 220Gb of XML data in the planet.osm file, covering the whole world, with thousands of updates every day, so if mapnik can render that I’m pretty certain it should be able to cope with my modest requirements.

Mapnik tiles in Open Street Map

Installing and configuring mapnik, unfortunately, turned out to be quite a challenge, and has taken me a significant amount of time to get a working installation. In fact, had I realised quite how many steps would be involved, I’d probably have taken more care to document them carefully. As it is, I’m going to try to note down in this post what I did while they’re still relatively fresh in my mind.

What about using pre-compiled Windows binaries?

Mapnik, like many open source applications, is primarily targeted at a UNIX stack. That means that a large part of the documentation will refer to commands like sudo apt-get, which will look fairly alien to many Windows users. There’s also a lot of dependencies on other packages – python, libboost, libpng, etc. which you may not be familiar with.

Now, fortunately, some kind people have taken the time to prepare pre-compiled Windows binaries of these various tools, and even packaged them together in convenient download format.

For example, the OSGeo4W package contains windows binary executables for Mapnik, along with the GDAL library (used for importing, converting of various spatial data formats), QGIS (open source desktop GIS application), and many other useful open source spatial goodies. The problem is that, with trying to keep track of all those separate dependencies, the package itself can become out-of-date quite quickly. The latest build of OSGeo4W, for example, is still based on python 2.5.2-1. The current build of the 2.x branch of python is 2.7.1, and even that represents the end-of-life release for the 2.x branch. The latest version of python is actually 3.2. I didn’t really want to start a project on software that had already been deprecated.

Likewise, the excellent FWTools project, which also bundles together many of the same packages as OSGeo4W still comes bundled with Python 2.3.4 and version 1.7 of the GDAL library. Crucially, for me, GDAL introduced support for SQL Server as a spatial data source only in version 1.8, so using FWTools wasn’t an option either.

What’s more, Mapnik itself seems to have forked into two versions – the current stable release being 0.7.1, but with many comments being made about breaking changes in the new 2.x development version. The only precompiled windows binaries I could find were of the 0.7 version and, again, I didn’t want to invest a lot of time setting up a project based on software that was about to go out of date.

So, precompiled binaries was, at this stage, a no-go.

Build-It-Yourself

I decided I was going to have to build my own installation, and here was my wish list:

  • Python 3.2
  • GDAL 1.8
  • Mapnik 2.x

Python was (thankfully) easy to install – there’s an installer package available from http://python.org/ftp/python/3.2/python-3.2.msi

GDAL was also not too bad – there are x86 and x64 packages (bundled with mapserver) that you can download from http://www.gisinternals.com/sdk/

Now onto Mapnik. And it is here, with retrospect, that I wish I’d stated taking notes. You can download the source for mapnik from here. However, before compiling mapnik itself, you need to download and/or compile its required dependencies. These are: proj4, boost, zlib, freetype, icu, libxml2, libpng, libjpeg, and libtiff.

Now, you could download the source for each of these and build them separately but fortunately, once again, some kind soul has done much of this work for us. If you go to the gnuwin32 project on sourceforge, you’ll find links to download most of these libraries. For those packages not included in gbuwin32, ICU is available from here, and you can get the latest zlib from here.

Once you’ve got all the dependencies sorted out, it’s onto the configuration changes. If you follow the article at http://trac.mapnik.org/wiki/Python3k you’ll see a number of steps required to rebuild the python bindings with Mapnik to target Python 3.x. After some fiddling about with these and a bit of guesswork, I managed eventually to get everything built.

Testing it Out

Mapnik comes with a test script so, gingerly, I tried running it. Lots of errors – couldn’t find xxx etc. I realised this was because I hadn’t set the environment variables and paths correctly so, after sorting this out, I had another go. This time looked a lot more promising:

image

And, lo and behold, here was the (beautiful) example image generated:

demo_high

Over-confident of my new found ability, I then tried altering the example rundemo.py script to point at a SQL Server datasource. Mapnik supports OGR datasources, so I first created a virtual layer that connected to my SQL Server. For the purposes of testing, I decided to select a set of data from the OS VectorMap District settlement area data (note that I’m not sure if OGR can deal with SQL Server’s native binary geography/geometry format, so I use STAsText() to get the WKT representation and then specify WKT encoding in the GeometryField):

<OGRVRTDataSource>
  <OGRVRTLayer name="AASQLlayer">
    <SrcDataSource>MSSQL:server=.\SQLEXPRESS;database=OSVectorMap;trusted_connection=yes</SrcDataSource> 
    <SrcSQL>SELECT geom27700.STAsText() AS geomWKT FROM TG11_Settlement_Area</SrcSQL>
    <GeometryField encoding="WKT" field="geomWKT"/>
  </OGRVRTLayer>
</OGRVRTDataSource>

Testing the virtual layer with ogrinfo seemed to suggest that everything was working ok:

image

So then I modifed the python script to add the new layer. Note that OS Vectormap data is defined using the OSGB British National Grid coordinate system (EPSG:27700), and mapnik expects the parameters for the srs property to use PROJ4 syntax, which you can get from http://spatialreference.org/ref/epsg/27700/proj4/:

vectormap_lyr = mapnik.Layer('OS Vectormap')
vectormap_lyr.srs = "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs"

vectormap_lyr.datasource = mapnik.Ogr(file='MSSQL.ovf',layer="AASQLlayer")

Unfortunately, here I hit another problem:

image

And it seems here my luck has run out, because I simply can’t work out how to get round this one. The error message is simply “Failed to open datasource”, yet ogrinfo confirms that the datasource is fine, and other GDAL/OGR components can read from it, so I don’t know if it’s because I built mapnik wrong or forgot to change/include a particular setting.

I did just notice that there’s a Google Summer of Code Project to improve Mapnik installation on Windows, and I’m really hoping it’s successful because the results generated by Mapnik are beautiful.

As a workaround, I actually used OGR2OGR to export my data from SQL Server into a shapefile called MSSQL_export.shp, and then used that as a datasource in Mapnik by changing the python datasource to:

vectormap_lyr.datasource = mapnik.Shapefile(file='MSSQL_export')

Finally, after a bit of XML styling, I was able to get the following image (click for full size):

image

I’m actually really pleased with the image quality, but until I can get Mapnik to retrieve data directly from SQL Server there’s not much point proceeding with the tilecutting process – I can’t really justify an additional step of exporting from SQL Server to shapefile just to get Mapnik to load it.

If anyone has had any success of getting Mapnik and SQL Server to play nicely together, please let me know!

Follow

Get every new post delivered to your Inbox.

Join 53 other followers