April 13, 2012

## Gridding Geometries (or, “Creating patchwork animals in SQL Server”)

This is one of those things that I can’t imagine anybody would ever really want to do but, seeing as I haven’t posted anything for a while, I thought I’d write about it just in case it’s of use to someone…

First, suppose you had an arbitrary geometry. Here’s an example:

```DECLARE @snake geometry = geometry::Parse('CIRCULARSTRING(0 0, 5 -5, 10 0, 15 5, 20 0, 25 -5, 29 -1)').STBuffer(2).STDifference(geometry::Parse('MULTIPOINT(28.5 -0.5, 29.75 -0.8)').STBuffer(0.6)).STUnion(geometry::Parse('MULTIPOINT(28.5 -0.5, 29.75 -0.8)').STBuffer(0.2)).STUnion(geometry::Parse('POLYGON((29.2 0.5, 29.5 2.7, 29.79 2.2, 30.22 2.75, 30 0.5, 29.2 0.5))'));
SELECT @snake;
```

In in the immortal words of Rolf Harris, “can you guess what it is yet?”.

Yes, that’s right – it’s a snake and it looks like this:

Now suppose that you wanted to divide that geometry up into a series of individual polygonal tiles on a regular-spaced grid. You could do so using the following T-SQL:

```-- Determine the extent of the grid required to cover the geometry
DECLARE @grid geometry = @snake.STEnvelope();
DECLARE @gridwidth float = @grid.STPointN(3).STX - @grid.STPointN(1).STX;
DECLARE @gridheight float = @grid.STPointN(3).STY - @grid.STPointN(2).STY;

-- Work out the lower-left corner of the grid
DECLARE @xoffset float = @grid.STPointN(1).STX,
@yoffset float = @grid.STPointN(1).STY;

-- Specify the number of columns and rows in the grid
DECLARE @numcolumns float = 80, @numrows float = 35;

-- Work out the height and width of each tile
DECLARE @cellwidth float = @gridwidth / @numcolumns;
DECLARE @cellheight float = @gridheight / @numrows;

-- Create a table variable to hold the tiles
DECLARE @gridcells table (id int, geom geometry)

-- Now, loop through and cut up the geometry into tiles
DECLARE @x int = 0, @y int = 0;
WHILE @y < @numrows
BEGIN
WHILE @x < @numcolumns
BEGIN
INSERT INTO @gridcells VALUES(
@y * @numcolumns + @x,
'POLYGON((' + cast(@xoffset + (@x * @cellwidth) AS varchar(32)) + ' ' + cast(@yoffset + (@y * @cellheight) AS varchar(32)) + ','
+ cast(@xoffset + ((@x + 1)  * @cellwidth) AS varchar(32)) + ' ' + cast(@yoffset + (@y * @cellheight) AS varchar(32)) + ','
+ cast(@xoffset + ((@x + 1)  * @cellwidth) AS varchar(32)) + ' ' + cast(@yoffset + ((@y + 1) * @cellheight) AS varchar(32)) + ','
+ cast(@xoffset + (@x * @cellwidth) AS varchar(32)) + ' ' + cast(@yoffset + ((@y + 1) * @cellheight) AS varchar(32)) + ','
+ cast(@xoffset + (@x * @cellwidth) AS varchar(32)) + ' ' + cast(@yoffset + (@y * @cellheight) AS varchar(32)) + '))'
)
SET @x = @x + 1;
END
SET @x = 0;
SET @y = @y + 1;
END

-- Finally, select the tiles from the table variable
SELECT geom.STIntersection(@snake) FROM @gridcells ORDER BY id;

```

My snake has now been chopped up into 800 little tiles, and here’s what the result looks like:

I wonder if this is how Elmer the Elephant was created?

March 28, 2012

## Rotating Bing Maps

The question of how to “rotate” a Bing Map comes up a surprising number of times on the MSDN Bing Maps forums.

The Bing Maps Silverlight Control has a Heading property, which, according to the documentation: “when overridden in a derived class, gets or sets the directional heading of the map”. However, I certainly can’t get it to have any effect, and it seems I’m not the only one.

 map.Heading = 0; map.Heading = 180; (Notice the difference? Nope, me neither.)

The Bing Maps AJAX v7 control also has a heading property, which is described here as “The directional heading of the map. The heading is represented in geometric degrees with 0 or 360 = North, 90 = East, 180 = South, and 270 = West”. What the documentation doesn’t mention is that the heading is only taken into account when using Birdseye view. So, in the examples below you can see four different aspects of Norwich castle in Birdseye view by changing the “heading” to represent north, east, south, and west:

However, you can choose from these four headings only in Birdseye view. If you’re using the aerial or road map style then the “heading” property, as in the Silverlight control, has no effect.

So, how can you create an aerial or road Bing Map with an arbitrary rotation, oriented so that “up” was southwest, or 40°, say?

## Rotating the Bing Maps Silverlight control

If you’re using the Silverlight control you can apply a RotateTransform to the Map element. The basic XAML syntax is as follows:

```<UserControl x:Class="BingMapsSilverlightRotation.MainPage"
xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation"
xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml"
xmlns:m="clr-namespace:Microsoft.Maps.MapControl;assembly=Microsoft.Maps.MapControl">

<Grid x:Name="LayoutRoot">
<m:Map CredentialsProvider="{StaticResource MyCredentials}">
<m:Map.RenderTransform>

<RotateTransform Angle="15" />

</m:Map.RenderTransform>
</m:Map>
</Grid>
</UserControl>
```

The map control (and all the elements contained within it) will then be rotated clockwise by the Angle specified in the RotateTransform. The result in this case is:

Note that, as for any of the techniques described here, both the background map image and the labels will be rotated. This is unavoidable because the labels are hard-baked into the tile images themselves (though you could of course use one of the unlabelled map styles and then add your own map labels from OpenStreetMap, say).

## Rotating the Bing Maps AJAX v7 control

To rotate the Bing Maps AJAX control, you can use the CSS3 transform property. Unfortunately (as is so often the case with web development), there is little standardisation between how different browsers implement the transform feature, and you’ll have to apply a range of variations on the theme to create a cross-browser solution. For example, the CSS styles required to rotate any element by 15 degrees clockwise across different browsers are as follows:

```transform:rotate(15deg); /* CSS3 "standard", but not actually used by any mainstream browser! */
-ms-transform:rotate(15deg); /* Internet Explorer */
-moz-transform:rotate(15deg); /* Firefox */
-webkit-transform:rotate(15deg); /* Safari and Chrome */
-o-transform:rotate(15deg); /* Opera */

```

You can apply these transforms directly onto the div element containing the map. The problem is that, as shown below, this will not only rotate the map image, but also the navbar, scalebar, and copyright notices:

To rotate only the map tiles, you need a more specific CSS selector to target those elements that should be rotated. Unfortunately, the various elements of the Bing Map control generally do not have unique IDs or classes, so you have to rely on child selectors. Fortunately, the navbar itself does have a class (MicrosoftNav), so you can ensure this is not rotated by using the not() selector. As at the time of writing, the following CSS selector targets all the divs that contain image tiles in a Bing Maps control created in the “mapDiv” div:

```#mapDiv > div.MicrosoftMap > div:not(.MicrosoftNav)

```

Usual caveats apply – this may break at any time in the future if the Bing Maps control gets updated yada, yada… but right now it works, as demonstrated in the code listing below tested in IE9, Firefox 11 and Chrome 17:

```<!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">
<title>Bing Maps Rotation</title>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<style type="text/css">
#mapDiv > div.MicrosoftMap > div:not(.MicrosoftNav)
{
transform:rotate(15deg);
-ms-transform:rotate(15deg); /* IE 9 */
-moz-transform:rotate(15deg); /* Firefox */
-webkit-transform:rotate(15deg); /* Safari and Chrome */
-o-transform:rotate(15deg); /* Opera */
}
</style>
<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: "AtGUpcoIpD5G3ty3r0wnfr!ggingk3yjg779XlZPsKWfhP",
center: new Microsoft.Maps.Location(54, -4),
mapTypeId: Microsoft.Maps.MapTypeId.aerial,
zoom: 4
});
}
</script>
<div id="mapDiv" style="position:relative; width:400px; height:300px;" />
</body>
</html>

```

Some further examples:

March 14, 2012

## Creating “The Legend of Zelda” map tiles from Bing Maps aerial imagery

Time for something light-hearted, I feel. Following on from last year’s Bing Maps tribute to Grand Theft Auto, I thought I’d do another videogame-themed map hack inspired by perhaps my most-loved game of all time – the Legend of Zelda.

Almost the entire year of my 10-year old life was spent playing the Legend of Zelda on the Nintendo Entertainment System. In those days, games rarely had their own in-game maps, and it was long before GameFAQs or Prima strategy guides existed, so the only way to keep track of where you’d been in an adventure game was to create your own maps. My mum would bring me home sheets of graph paper from work, and I would spend hours laboriously drawing the layout of each room in every dungeon, and how each of the screens in the overworld connected to each other. (If you’ve ever tried to do this yourself, you’ll appreciate the ingenuity of the “lost woods”, which are very hard to map on paper!).

As a result, I’ve always personally found the overworld maps of Zelda games to be a thing of beauty. Here’s the world map for the original Legend Of Zelda:

and here’s The Legend of Zelda: A Link to the Past:

The world of Hyrule is undeniably beautiful – its deserts, rivers, and forests becoming ever more intricate in each successive game in the Zelda series. But wouldn’t it be cool to create an 8bit-style tileset of the real world that you navigate around?

Now, I’m no graphic artist and, even if I was, the thought of drawing the 275 billion tiles required to cover the world at zoom level 19 would probably put me off the idea anyway. So, I need some way to automatically determine whether a given area should be represented as a water tile, or a grass tile, say. And I wondered if it was possible to do this based on a simple image analysis of the corresponding aerial tile from Bing Maps.

## Converting Bing Maps Aerial Imagery to 8bit tiles

So, here’s the plan: I build a Bing Maps application that uses a custom tile layer. The handler for the tile layer requests aerial imagery from Bing Maps. But, instead of returning the tile image direct to the client, it converts the features on that tile into appropriate 8bit sprites and returns a new tile composed of those sprites to the client instead.

To choose which tile images to use, I’ll divide each 256px x 256 tile retrieved from Bing Maps into smaller 32px x 32px subtiles. I’ll then do a simple RGB analysis to determine the “average” colour of each subtile.

 If the average colour of that subtile is mostly green (that is, when comparing the RGB component values, G > B > R) then I’ll insert a grass tile: If the average colour is mostly blue (that is, B > R and B > G) then insert a water tile: If the average colour is mostly dark brown (R >= G and R >= B and B <= 90) then insert a rock tile: If the average colour is mostly light brown (R >= G and R >= B and B > 90) then insert a dirt/sand tile: To add a bit of variety, I’ll also make it so that every grass tile has a 1/50 chance of being replaced with a flower tile:

You can obviously tweak these, or add more rules to fit the tileset as you see fit.

Here’s the code to implement these rules in a C# handler:

```<%@ WebHandler Language="C#" Class="Handler" %>

using System;
using System.Web;
using System.Drawing;
using System.Data.SQLite;
using System.IO;
using System.Drawing.Imaging;

public class Handler : IHttpHandler
{

public void ProcessRequest(HttpContext context)
{

// Create a random number generator to vary grass tiles
Random random = new Random();

// Retrieve the original aerial tile image
System.Net.WebResponse response = request.GetResponse();
System.IO.Stream responseStream = response.GetResponseStream();
Bitmap mapImage = new Bitmap(responseStream);

// Set up a graphics class based on the tile bitmap
using (Graphics graphics = System.Drawing.Graphics.FromImage(mapImage))
{
int tileWidth = 32, tileHeight=32;
int numRows = 8, numCols = 8;

// Loop through the aerial tile in 16x16 squares
for (int row = 0; row < numRows; row++)
{
for (int col = 0; col < numCols; col++)
{
// Get the bitmap data for this square
Rectangle tileRect = new Rectangle(col * tileWidth, row * tileHeight, tileWidth, tileHeight);
System.Drawing.Imaging.BitmapData bmData = mapImage.LockBits(
tileRect,
mapImage.PixelFormat);

// Get the average R,G,B values of pixels in this square
long[] totals = new long[] { 0, 0, 0 };
unsafe
{
// 24bit image so 3 bytes per pixel (PNG + transparency would be 4)
int PixelSize = 3;
for (int y = 0; y < tileHeight; y++)
{
byte* p = (byte*)bmData.Scan0 + (y * bmData.Stride);
for (int x = 0; x < tileWidth; x++)
{
totals[0] += p[x * PixelSize]; // Blue
totals[1] += p[x * PixelSize + 1]; // Green
totals[2] += p[x * PixelSize + 2]; // Red
}
}
}
mapImage.UnlockBits(bmData);

// Work out the average RGB colour value
int avgB = (int)(totals[0] / (bmData.Width * bmData.Height));
int avgG = (int)(totals[1] / (bmData.Width * bmData.Height));
int avgR = (int)(totals[2] / (bmData.Width * bmData.Height));

// Determine average colour of this tile
Color fillColour = Color.FromArgb(avgR, avgG, avgB);

// Snow
if (fillColour.G >= 225 && fillColour.B >= 225 && fillColour.R >= 225)
{
graphics.FillRectangle(Brushes.White, tileRect);
}

// Grass
else if (fillColour.G >= fillColour.B && fillColour.G >= fillColour.R)
{
// Random chance of being a flower tile
if (random.Next(0, 50) < 1)
{
Bitmap grassImage = (Bitmap)Image.FromFile(@"C:\Users\Alastair\Documents\Visual Studio 2008\WebSites\BingMaps\v7\ZeldaMap\flowers.png");
graphics.DrawImage(grassImage, tileRect);
}
else
{
Bitmap grassImage = (Bitmap)Image.FromFile(@"C:\Users\Alastair\Documents\Visual Studio 2008\WebSites\BingMaps\v7\ZeldaMap\grass.png");
graphics.DrawImage(grassImage, tileRect);
}
}

// Water
else if (fillColour.B >= fillColour.G && fillColour.B >= fillColour.R)
{
graphics.FillRectangle(new SolidBrush(Color.FromArgb(255,48,89,166)), tileRect);
}

// Dirt/Sand
else if (fillColour.R >= fillColour.G && fillColour.R >= fillColour.B && fillColour.B <= 90)
{
Bitmap dirtImage = (Bitmap)Image.FromFile(@"C:\Users\Alastair\Documents\Visual Studio 2008\WebSites\BingMaps\v7\ZeldaMap\rock.png");
graphics.DrawImage(dirtImage, tileRect);
}

// Rock
else if (fillColour.R >= fillColour.G && fillColour.R >= fillColour.B && fillColour.B > 90)
{
Bitmap rockImage = (Bitmap)Image.FromFile(@"C:\Users\Alastair\Documents\Visual Studio 2008\WebSites\BingMaps\v7\ZeldaMap\dirt.png");
graphics.DrawImage(rockImage, tileRect);
}

// Default - just fill tile in single average colour
else {
graphics.FillRectangle(new SolidBrush(fillColour), tileRect);
}
}
}
}

// Send the resulting 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);
}

public bool IsReusable
{
get
{
return false;
}
}

}

```

(Note that, because I’m accessing aerial tiles directly from the Bing Maps tile servers rather than through the API, this method is unsupported and probably not allowed in a production environment. But then, why on earth would you ever want to do this in a production environment? It’s just for fun!)

Based on the rules above, let’s see how certain tiles are rendered:

 0331 0331 1 01 01

That’s pretty much the effect I was after, and here’s what it looks like then you create a custom tile layer based on the above handler displayed on a Bing Maps control at zoom level 1. Pretty neat, huh?:

To finish this example off, I want to add a player to the world as well. He’ll stay in the middle of the map, and the tiles will scroll under his feet as you press the cursor keys for him to navigate around. I’ll create a few simple frames of animation to show him walking in different directions, as shown in the sprite strip to the left.

The player sprite will appear in a div placed in the centre of the map, and I’ll use javascript to change the background of the div to the appropriate frame depending on whether the player is moving up, right, down, or left. Since there are three frames of animation for each direction, I’ll create a simple frame matrix that will loop through the relevant frames.

var frameMatrix = [1, 2, 0, 4, 5, 3, 7, 8, 6, 10, 11, 9];

In addition to animating the background of the player sprite div, pressing a cursor key will also cause the map to scroll in the appropriate direction. For this, I’ve added a custom Pan() method that extends the Microsoft.Maps.Map prototype.

Here’s the code listing for the HTML page:

```<!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">
<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">

// Define some globals for handling key input
var keys = { UP: 38, DOWN: 40, LEFT: 37, RIGHT: 39 };
var keysPressed = [];
// These used to track player animation frames
var frameMatrix = [1, 2, 0, 4, 5, 3, 7, 8, 6, 10, 11, 9];
var currentFrame = 0;
// The map. Obviously.
var map = null;
// Add a convenient method to pan the map
Microsoft.Maps.Map.prototype.Pan = function(x, y) {
var pos = this.tryPixelToLocation(new Microsoft.Maps.Point(this.getWidth() / 2 + x, this.getHeight() / 2 - y), Microsoft.Maps.PixelReference.control);
this.setView({ center: pos });
};
// Set the frame rate of how often the map moves/player animates
var int = window.setInterval("animate()", 150);

// This function gets called once per animation loop
function animate() {
// Only look at the first key pressed
switch (keysPressed[0]) {
case keys.UP:
// If we've reached the end of the "up" animation frames, start again
if (currentFrame >= 2) { currentFrame = 0; }
// otherwise, go on to the next "up" animation frame
else { currentFrame++ };
// Pan the map up
map.Pan(0, 10);
break;
case keys.DOWN:
// If we've reached the end of the "down" animation frames, start again
if (currentFrame < 3 || currentFrame >= 5) { currentFrame = 3; }
// otherwise, go on to the next "down" animation frame
else { currentFrame++ };
// Pan the map down
map.Pan(0, -10);
break;
case keys.LEFT:
// If we've reached the end of the "left" animation frames, start again
if (currentFrame < 6 || currentFrame >= 8) { currentFrame = 6; }
// otherwise, go on to the next "left" animation frame
else { currentFrame++ };
// Pan the map left
map.Pan(-10, 0);
break;
case keys.RIGHT:
// If we've reached the end of the "right" animation frames, start again
if (currentFrame < 9 || currentFrame >= 11) { currentFrame = 9; }
// otherwise, go on to the next "right" animation frame
else { currentFrame++ };
// Pan the map right
map.Pan(10, 0);
break;
}

// Update the offset of the background of the player div
document.getElementById("player").style.backgroundPosition = "0 " + -currentFrame * 32 + "px";
}

// This function keeps track of the key(s) currently being pressed.
// Uses an array so that, for example, can use diagonal movement when
// two keys are pressed.
function handleKeyDown(e) {
switch (e.keyCode) {
case keys.UP:
if (!inArray(keys.UP, keysPressed) && !inArray(keys.DOWN, keysPressed)) {
keysPressed.push(keys.UP);
}
break;
case keys.DOWN:
if (!inArray(keys.DOWN, keysPressed) && !inArray(keys.UP, keysPressed)) {
keysPressed.push(keys.DOWN);
}
break;
case keys.LEFT:
if (!inArray(keys.LEFT, keysPressed) && !inArray(keys.RIGHT, keysPressed)) {
keysPressed.push(keys.LEFT);
}
break;
case keys.RIGHT: // Right:
if (!inArray(keys.RIGHT, keysPressed) && !inArray(keys.LEFT, keysPressed)) {
keysPressed.push(keys.RIGHT);
}
break;
}
e.handled = true;
}

// Remove key from the array when they are released
function handleKeyUp(e) {
keysPressed = removeFromArray(e.keyCode, keysPressed);
}

// Utility methods
function inArray(element, arr) {
for (var i = 0; i < arr.length; i++) {
if (arr[i] === element) {
return true;
}
}
return false;
}
function removeFromArray(element, arr) {
for (var i = 0; i < arr.length; i++) {
if (arr[i] == element)
arr.splice(i, 1);
}
return arr;
}

// Create the map
function GetMap() {
map = new Microsoft.Maps.Map(document.getElementById("mapDiv"),
{ credentials: "GETYEROWNBINGMAPSKEY!",
center: new Microsoft.Maps.Location(50, 0),
mapTypeId: Microsoft.Maps.MapTypeId.mercator,
zoom: 5,
inertiaIntensity: 0.00001,
tileBuffer:3,
disableKeyboardInput: true
});

// Create a custom tile source that points to the .NET handler
var tileSource = new Microsoft.Maps.TileSource({ uriConstructor: location.href.substring(0, location.href.lastIndexOf('/')) + "/" + "handler.ashx?q={quadkey}" });

// 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);

// Position the player sprite in the middle of the map
document.getElementById("player").style.backgroundImage = "url('player.png')";
document.getElementById("player").style.left = map.getWidth() / 2 + "px";
document.getElementById("player").style.top = map.getHeight() / 2 + "px";

// Attach event handlers
}
</script>
<div id='mapDiv' style="position:relative; width:640px; height:480px;"></div>
<div id='player' style="position: absolute; z-index:10; width:20px; height:32px; background-repeat:no-repeat;"></div>
</body>
</html>

```

Finally, here’s a video of the finished product. Watch me walk from somewhere near Paris to Valencia and then onto Naples, all in glorious dynamically-generated 8-bit graphics!

March 11, 2012

## “Point-of-Maximal Distance” queries (AKA Variations on a Theme of Nearest-Neighbours)

Nearest neighbour queries (or, for my international readers, nearest-neighbor queries), are a pretty common query pattern used in location-based applications and spatial analysis:- given a set of known locations of restaurants, customers, stations etc.,  “What/where is the nearest X to a given location?” I’ve already written about performing nearest-neighbour queries with SQL Server spatial datatypes here and here.

Occasionally, you’ll see this pattern reversed – “What is the X furthest away from my current location?” – which could be referred to as a farthest-neighbour query.

However, there’s another variant that I’ve not often seen discussed, which is: “Given a set of known locations, where is the point on Earth that maximises the distance to the closest X.” In other words, where would you have to go to get as far away as possible from any item in the given set of data? The key difference to appreciate about this query structure is that, unlike the first two examples, the answer is not selected from one of a set of known locations – it is location that must be determined programmatically to maximise the distance from X. This makes the query a lot more complicated, but also more interesting. (I don’t know if there is a generic name for this type of query or not…  a “point-of-maximal distance” query?).

I can think of several potential practical applications for this type of query – a fire response service may want to identify the location at increased risk because it is furthest away from any of their response stations, for example. Or, a mobile telecommunications operator might want to identify the “darkest spot” in their network coverage – the point that lies furthest from any of their transmitters.

As ever, I think a practical example might help demonstrate. One of my recent posts described how to import a set of point data representing Tescos stores in the UK. People that know me will know that I’m a Green, tree-hugging, anti-capitalist hippie, so let’s suppose that I wanted to move house to get as far away from Tescos as possible (for my international readers, Tesco is a multi-billion pound supermarket chain in the UK often presented as one of the faces of modern global capitalism – I think you have similar companies like Wal-Mart or Subway).

I’m going to set an additional constraint that, while I want to live far from Tescos, I still want to live in mainland UK. So, my problem is to find the point that lies somewhere inside this Polygon that is furthest away from any of the dots (which are the locations of Tescos stores):

So…. how do that?

## Method 1: Iterative T-SQL

To start, I’m going to create the shape of the UK as a geography MultiPolygon called @UK. I’ve already got a table containing various parts of the UK that I can aggregate, but I want to exclude the Scottish Islands and other small polygon geometries that will confuse matters (Scottish islands = good for single malt whiskey, bad for spatial data… they’re just too craggy!). So, I’ll create a simple numbers table in a CTE and join to that table in order to loop through each element in the UK table and only aggregate those polygons that are greater than 10,000km2 in area:

```DECLARE @UK geography;

-- Create a MultiPolygon of the UK, omitting any islands
-- with area less than 10,000km2
WITH Numbers(n) AS (
SELECT ROW_NUMBER() OVER (ORDER BY number) FROM master..spt_values
)
SELECT @UK = geography::UnionAggregate(Shape.STGeometryN(n))
FROM UK JOIN Numbers ON Numbers.n  10000000000;

```

I’ll also create a geography MultiPoint of all Tescos stores called @Tescos. The only slight complication here is that the Garmin POI dataset I’m using seems to have some erroneous data, so I’ll put a WHERE condition to filter only those records that lie within the approximate geographic bounds of the UK:

```DECLARE @Tescos geography;

-- Create a MultiPoint of Tescos stores
SELECT @Tescos = geography::UnionAggregate(Location)
FROM TescoStores
WHERE Location.Long BETWEEN -10 AND 3
AND Location.Lat BETWEEN 50 AND 60;
```

Now that we’ve got the basics set up, how do we go about finding the part of @UK that lies furthest from @Tescos? To explain one approach, consider an analogy in which the shape of the UK is cut out of blotting paper. You then arrange a set of pipettes filled with ink, so that each is positioned above the location of a Tescos store. And then, all at the same time, you start dropping ink from the pipettes onto the location of the stores.

The ink blots centred on the location of each store will expand outwards at the same rate, gradually filling up the blotting paper shape of the UK. The last part of the paper to be covered by ink is that which lies furthest from any of the ink sources, and must be the solution.

We can recreate this approach using SQL Server’s spatial functions:

• To mimic the expanding ink blots we can call STBuffer() on the @Tescos MultiPoint with increasing buffer sizes to create expanding circular zones centred around each store.
• To determine the area of blotting paper that is still not covered by ink, we can use STDifference() to compare the inked area to the original shape.
• If the part of the UK yet to be covered by ink is still of a sufficiently large area, we simply add more ink by looping through again with an increased buffer size. If, however, the whole UK has been covered by ink then we’ve put too much ink on and we need to reduce the buffer size slightly.
• When the area yet to be covered reaches a sufficiently small size (this approach is unlikely to ever converge on a single Point), declare that to be the solution.

Here’s the query:

```DECLARE @x float = 100000; -- Declare an initial buffer size. I'm guessing that the furthest point in the UK lies more than 100km from Tescos
WHILE @UK.STDifference(@Tescos.STBuffer(@x)).STArea() > 10000 -- Stop when we find an area less than 10km2 in area
BEGIN

-- Print a little debug message to see what's happening on this iteration
PRINT CAST(@UK.STDifference(@Tescos.STBuffer(@x)).STArea() AS varchar(32)) +
'm2 of the UK lies more than ' + CAST(@x AS varchar(32)) + 'm from the closest Tescos';

-- Amount by which to increase the buffer size
DECLARE @y float;

-- Fairly arbitrary, but enlarge the buffer by the square root of the remaining area
-- Intention here is to close in on the appropriate buffer value relatively quickly
SET @y = SQRT(@UK.STDifference(@Tescos.STBuffer(@x)).STArea());

-- Check whether increasing the buffer by this amount would be too great
-- and lead to the whole UK being covered in ink!
WHILE @UK.STDifference(@Tescos.STBuffer(@x + @y)).STArea() = 0 BEGIN

-- Reduce the amount by which the buffer is increased
SET @y = @y/2;
END

-- Apply the increase in buffer
SET @x = @x + @y;

END

```

Once this loop has finished, @x will contain a value, measured in metres, that represents the radius to which each of the ink blots must expand so that the area left uncovered has an area of less than 10km2 . Here’s what the map looks like when a circular buffer of radius @x is created centred on the location of each Tesco store:

```SELECT @Tescos.STBuffer(@x)
UNION ALL SELECT @Tescos
UNION ALL SELECT @UK;
```

As you can tell, the whole of the UK is pretty much covered – the iteration only stops when the area remaining uncovered is less than 10km2 in in size (WHILE @UK.STDifference(@Tescos.STBuffer(@x)).STArea() > 10000), which means that, at this zoom level, it’s going to be pretty hard to perceive.

To find where the remaining un-inked area is, we can SELECT @UK.STDifference(@Tescos.STBuffer(@x)), and, to make it more obvious, I’ll also add a buffer to enlarge the solution area. And the answer? (as you might have been able to tell from the previous image) is the small blue polygon shown below:

Perhaps unsurprisingly, it’s an area at the tip of the highlands of Scotland, and it’s 119,591 metres from the closest Tesco store. Ok, so it’s not an exact point, but the way that this iteration works it would be very hard to reach such a precise location.

## Method 2: SQLCLR Triangulation?

My gut instinct told me that there should be a way of coming up with the exact solution to the problem, in a much simpler manner, using triangulation. For example, suppose you were to create the Delaunay triangulation of the set of points such that the vertices of each triangle were formed from the location of the Tesco stores:

Since the triangles are formed by connecting vertices of Tesco stores, in those areas where stores are further apart, the circumcircle (the point of equal distance from all three vertices) of the triangle created between them will be larger. In fact, you’d think the point that lies at the centre of the triangle of largest circumradius would be the point that lies furthest from any Tesco store. And this would be true, were it not for (quite literally) “edge cases”…

… when, as in this case, the solution lies on the edge of the area in question, it is not covered by the extent of the triangulation (which only creates triangles that are interior to the set of data). Only by artificially extending the triangles out from the set of Tesco stores towards infinity could we ensure that the whole UK had been covered, and it turns out the SQL Server is not too happy about creating infinite-sized triangles (who’d have thought?).

So that’s a bit disappointing… the iterative T-SQL method still works and, with a bit of tweaking to choose how the buffer should be refined on each iteration, could probably be made to work more efficiently. However, I was hoping to be able to get something exact using SQLCLR.

I should just point out that the fact the triangulation method doesn’t extend beyond the bounds of the supplied dataset can still make it useful for lots of other purposes. For example, creating an aggregate union of the triangles creates a concave hull around the location of all Tescos stores in the UK (not a convex hull as could be created by STConvexHull()). The result looks scarily like the shape of the UK itself:

And, deriving a Voronoi tessellation from the set of triangles creates the “catchment area” around each store – the area formed from all those points closest to that store than any other. You can think of this as the materialised form of the results of a nearest-neighbour query:

(Note that, as expected, the point that lies furthest away from a Tescos store is contained within the largest Voronoi polygon – the pinkish one at the top of the map). At first, the Voronoi tessellation and the Delaunay triangulation may appear similar – indeed, they are “dual structures” – while Tesco stores form the vertices of each triangle in the Delaunay triangulation, they are the sites around which each polygon is created in the Voronoi tessellation.

If anybody can suggest any alternative approaches to identify the “point of maximal distance” from a set of data, I’d love to hear it…

Note: The SQLCLR code for creating the Delaunay triangulation, concave hull, and Voronoi tessellation demonstrated in this post is taken from “Pro Spatial with SQL Server 2012”, which also contains a chapter explaining these topics in more detail.