Creating Hill Shaded Tile Overlays

The default Bing Maps road style uses a “hillshade” effect to give an impression of underlying terrain. It’s a relatively subtle, but surprisingly powerful technique to enhance the appearance of map layers, as demonstrated by comparing the following two tiles:

Without hillshading

With hillshading

In this post, I’ll describe how to create your own hillshade overlay from digital elevation model (DEM) data, using the GDAL toolset.  By creating the overlay as a set of semi-transparent tiles, rather than pre-rendered into the tiles as shown above, you can place them on top of any Bing Maps/Google Maps et al. tilelayer to represent the underlying terrain.

The process I’ve followed is based on the work of others, most notably PerryGeo, and you can find some other guides on the internet to achieve this same effect. However, I found some of the existing guides on the subject to be either out-of-date or require knowledge of Linux BASH commands etc., so I hope that some of you will find this new step-by-step guide helpful.

1.) Acquire a DEM terrain model

To start with, you’re going to need some source data about the underlying terrain of the earth from which to calculate your hillshade. There’s lots of places to acquire this data from; Perhaps the easiest to use (assuming you’ve got Google Earth installed) is to open the kmz file available from http://www.ambiotek.com/topoview. This uses Google Earth as a graphical interface for v4.1 of the  elevation dataset gathered by the Shuttle Radar Topography Mission (SRTM), from which you can click to download individual DEM tiles covering 5°x 5°, as shown below:

Alternatively, you can access these files directly from the KCL server (my former university, incidentally) at http://srtm.geog.kcl.ac.uk/portal/srtm41/srtm_data_geotiff/

The data is provided in GeoTIFF format. You can load one of these tiles up in any graphics program that can load TIFF files, but it won’t look very interesting yet. The height information is encoded in additional metadata that will be ignored by normal graphics programs, so you’ll probably just get an image like this (this is srtm_36_02.tif):

Black parts show the presence of data in the underlying file, which we’ll subsequently process using GDAL tools to create shaded images.

2.) Reproject to Spherical Mercator

Most DEM data sources, including the SRTM data I linked to above, are provided in Plate Carree projection – i.e. WGS84 coordinates of longitude are mapped directly to the x axis of the image, while latitude is mapped directly to the y axis. Before we create tiles from this data suitable for overlay on Bing Maps, Google Maps, et al. we therefore need to transform it into the Spherical Mercator projection. You can do this using gdalwarp, as follows:

gdalwarp -dstnodata 0 -tr 305.7481 305.7481 -multi -co "TILED=YES" -t_srs EPSG:3857 srtm_36_02.tif srtm_36_02_warped.tif

The full list of parameters accepted by gdalwarp are listed here,  but the options I set are as follows:

• dstnodata states what value to use to represent nodata values (the equivalent of null in a SQL database, for example). I’ve set a value of 0 (i.e. black).
• tr gives the target resolution in x and y dimensions. The SRTM data I’m using was recorded at 90m resolution, so you might think that this should be set to 90 90. However, I’m going to be using this data for display on Bing Maps at different zoom levels, which will necessarily involve resampling the image.  Therefore, you should set this value to the resolution (in metres/pixel) of the maximum zoom level on which you plan to overlay your data. (Remember that maximum zoom level will have the smallest resolution). You can obtain this value from my Bing Maps Ready Reckoner. In the case above, I’m planning overlaying my data on Zoom Level 9 and above, so I set a value of 305.7481 (in both dimensions). If I’d wanted to go to Zoom Level 10, I would have decreased this to 152.87 instead.
• multi allows parallel processing
• co “TILED=YES” is a format-specific option that states that the output TIFF file should be tiled rather than stripped (see http://www.fileformat.info/format/tiff/egff.htm for an explanation of the difference)
• t_srs gives the destination spatial reference system into which the image should be reprojected. In this case, EPSG:3857, as used by Bing Maps, Google Maps etc.

The resulting image, srtm_36_02_warped.tif, will still be a GeoTIFF file, but will now be projected as follows. The height and width of the output image will depend on the target resolution you specified in the tr parameter:

3.) Convert from DEM to Hillshade

The warped GeoTIFF file has height data encoded in it, but we want to translate that information into a visible shaded effect, and for this we can use gdaldem.

gdaldem actually provides several interesting functions related to working with DEM data, including the ability to derive contour lines, and create shaded relief maps. Maybe I’ll write about these another time, but for this example we want to use the hillshade mode. You can shade the warped image created in the previous step as follows:

gdaldem hillshade srtm_36_02_warped.tif srtm_36_02_warped_hillshade.tif -z 2 -co "TFW=YES"

This time, I’m only supplying two additional parameters:

• z is a scaling factor applied to the generated hillside image that accentuates the hills, increasing the contrast of the image. I provided a value of 2 just to enhance the effect a bit, but you might decide you don’t need this.
• co “TFW=YES” specifies that the output image should be created with an accompanying “world file”. This is a simple ASCII text file that provides additional information about the geographic extents of the created image, which we’ll need to use in a later step to line the hillshade image up with the Bing Maps tiling system. You can look up more information about world files on wikipedia.

There are additional parameters that allow you to specify the direction and the angle of the light source from which the simulated shadows will be created.

The result of executing the above code will be another TIFF file, in which the background is black, and the elevation data from the DEM has been converted into shades of grey, as follows:

At this stage, you could stop if you wanted to, and simply create a tile layer from the hillshaded image, which would look a bit like this:

Which makes the landscape of North Wales look a bit like the Moon, I think…

To make the data slightly more usable, we need to carry on with a few more tweaks.

4.) Making a Semi-Transparent Overlay

Currently, our hillshade image is opaque, with the shadows cast by terrain represented by variations in brightness of the colour used. To make this into an re-usable overlay that can be used on top of other layers, we need to make the image semitransparent, with shadows cast by terrain being represented by variations in opacity instead.

There are several ways of modifying the image data to achieve this effect. You could do it in Photoshop or another graphics program, for example, or using the graphics libraries in C# or PHP. Since I’m currently trying to learn Python, and GDAL is quite closely linked with Python, I’ll try to do it using the Python Imaging Library instead.

The following Python script makes a number of tweaks to the image above. Firstly, it converts it to a pure greyscale image (while the image above looks greyscale, it’s actually using a colour palette). It then inverts the image, turning it into a negative image. The reason for the inversion is that we then copy the (single) channel of the greyscale image into the opacity channel of a new RGBA image – areas that were very light in the source want to have very low opacity in the transparent image, and vice-versa, so the channel needs to be inverted.

Finally, we scan through the data to find instances of pixels that are pure black (RGBA value 0, 0, 0, 255) –this was the nodata value we set in step one – and replace them with pure transparent pixels (0, 0, 0, 0). The alpha channel in the tuples of any other pixels is also lightened slightly – I chose a value of 74 somewhat arbitrarily because I thought the resulting image looked good – you can choose whatever value you want, or none at all.

from PIL import Image as PImage
from PIL import ImageOps

# Load the source file
src = PImage.open("srtm_36_02_warped_hillshade.tif")

# Convert to single channel
grey = ImageOps.grayscale(src)

# Make negative image
neg = ImageOps.invert(grey)

# Split channels
bands = neg.split()

# Create a new (black) image
black = PImage.new('RGBA', src.size)

# Copy inverted source into alpha channel of black image
black.putalpha(bands[0])

# Return a pixel access object that can be used to read and modify pixels
pixdata = black.load()

# Loop through image data
for y in xrange(black.size[1]):
  for x in xrange(black.size[0]):
    # Replace black pixels with pure transparent    
    if pixdata[x, y] == (0, 0, 0, 255):
      pixdata[x, y] = (0, 0, 0, 0) 
    # Lighten pixels slightly
    else:
      a = pixdata[x, y]
      pixdata[x, y] =  a[:-1] + (a[-1]-74,)

# Save as PNG
black.save("srtm_36_02_warped_hillshade_alpha.png", "png")

(Much of the logic in this script came from here). The resulting image will be a PNG file, in which darker shadows are represented by increasingly opaque black parts, while lighter shadows are more transparent:

5.) Cut into Tiles

Now, to cut the PNG image into tiles. The only problem is that, having converted from a TIF to a PNG in the Python script above, we have lost the geo-referencing information that, up to that point, had been stored in the GeoTIFF metadata associated with the image. PNG files only contain image data, so we need to provide the associated geo information in a separate world file. Fortunately, (assuming you haven’t done any cropping, rotating, or resampling in the meantime), the PNG file above still (geographically) matches the hillshade TIFF file created in step 3, which, as you’ll recall, we created with an associated world file. So, all we need to do is take a copy of that world file and associate it with the PNG file instead:

copy srtm_36_02_warped_hillshade.tfw srtm_36_02_warped_hillshade_alpha.pngw

Naming the file with exactly the same filename as the PNG image file itself but with the .pngw extension means that GDAL will automatically pick up the associated geo-information when it tries to use the PNG file.

The final step, then, is just to use the gdal2tiles.py script to cut the PNG image into tiles:

python gdal2tiles.py srtm_36_02_warped_hillshade_alpha.png --s_srs EPSG:3857

Let it whir away until you’ve got directories of tiles stored in the TMS z/x/y structure:

And then you can create a new TileLayer that points to this set of tiles and add it to your map, on top of any other layer. As I said right at the very top of this article, the default Bing Maps road style already has hillshading effects added, so there’s not much point overlaying this tile layer on top of the road style. But here’s what it looks like placed on top of some other map styles (click to enlarge):

Open Street Map

Without hillshading overlay

With hillshading overlay

Ordnance Survey MiniScale

Without hillshading overlay

With hillshading overlay

In Defence of 3D Charts…

There was a time, not so long ago, when pretty much all business graphics looked like this:

 

The global hegemony of Microsoft Excel (particularly, Excel 97), deployed onto millions and millions of users’ desktop PCs meant that the world was awash with poorly-labelled, hideously-coloured, and sometimes just plain wrong line charts and bar graphs. I’m not blaming Excel itself – it was a fantastically successful product – but the problem was that a huge number of people were producing business reports who had no training on the software, little background in statistics itself, and certainly no explicit knowledge of data visualisation methods (which, at the time, were unheard of as an academic discipline).

But then things got even worse, as people discovered the ability to create “3d” charts. All of a sudden, the unattractive but essentially harmless business graphics above morphed into things like this:

 

The general criticisms of these “3d” charts are quite well known, so I’m not going to repeat them here. However, to give two specific comments on the examples above:

• The artificial “perspective” added to the pie chart distorts the true relationship between the values in the series. For example, what is the relative size of the large purple slice compared to the blue slice? You can’t tell because the purple slice is shown closer to the viewer, enlarging its relative area and giving it distorted prominence.
• In the line chart, the shadow effects added to the extruded ends of each column make it difficult to read across to the axis to find out the exact value of an item. Furthermore, in several places the blue columns in the foreground (inquiries) obscure the red columns behind (sales).

“3d” charts enjoyed a brief life as the popular visualisation of choice, seen to make data more exciting and, thus, making managers pay more attention to it (also perhaps in the hope that adding “depth” to a graphic added depth to the corresponding data analysis?!).

However, before too long, they were debunked as being nothing more than marketing gloss. In fact, much more than that, 3d charts were proclaimed as being misleading and dangerous, and were wholeheartedly rejected by most members of the data visualisation community. In a write-up of one of the Jen Stirrup’s sessions from SQLBits8 conference last month, Luke Merrett describes the learning point he takes away as “…you should never, EVER use 3D graphs for reporting…”.

I wonder if 3d charts, which not so long ago were being enthusiastically embraced and adopted, have now been wholeheartedly disregarded with equal lack of consideration. So, I’d like to take this opportunity to stick up for 3d graphs and explain why, in certain situations, I believe they are an appropriate way to visualise sets of data.

“3D” Graphs or 3D Graphs?

Before I do, I need to clarify something. In the preceding paragraphs, you might have noticed that I carefully referred to “3d” graphs, not 3d graphs. Why?

As far as I’m concerned, the examples of the tilted pie chart and the extruded column chart with shadows shown above are not 3d charts. They are 2d charts that have had a graphic effect applied to them to give them a sense of depth. That effect is not in any way linked to the underlying data set – it has no corresponding dimension (from a data point of view). These arbitrary graphic effects take up space in the chart without adding any further information about the data, which is pretty much the definition of “chart junk”.

I referred to these images as “3d” charts, but they should really be described as two-and-a-half dimensional charts. Clearly (unless you’ve got some sort of very expensive holographic-projection monitor), it’s impossible to create a chart that is displayed in three-dimensional space. However, it is possible to create a 2-dimensional representation of a chart that represents data in 3 dimensions, and that is what I regard as a 3d chart.

As an example, yesterday I was doing some performance testing of spatial indexes in SQL Server Denali. There are lots of factors that can influence the performance of SQL Server spatial  queries, but I decided to focus on two factors – the number of rows in the underlying base table, and the selectivity of the query window (measured as a percentage of the bounding box). I ran a test query repeatedly, changing these two independent variables and measuring the execution time of the query (my dependent variable).

At the end of my tests, I had gathered over 15,000 individual test results, and I decided to plot them using a 3d surface chart as follows:

Why do I think this 3d graph is suitable for this data?

• It’s clear to see, at a glance, the way in which, separately, the number of rows in the table and the size of the query window affects execution time, and also the way in which they affect it in combination. Since I’m using only a single set of axes for all my data, you can easily see which of the independent variables has a greater effect on the dependent variable.
• With the exception of a few spikes, the dependent variable (execution time) increases monotonically with both independent variables. That is to say, as the number of rows in the base table increases and the size of the query window increases, the execution time will always increase (not necessarily at a consistent rate, but it doesn’t fall). This means that, so long as the axes are oriented correctly, there is no risk of a higher value in the foreground of the chart obscuring a smaller value behind it.
• I don’t need to know exact values corresponding to any point on this graph. The results illustrated in this chart were obtained from a test system and are intended to be illustrative of the shape of the graph – the exact number of milliseconds taken to execute a query will depend on the exact configuration of the server and the query in question anyway, so it is not important for people to be able to tell these. (In fact, I could remove the z axis labels altogether, but I’ve left them in to demonstrate that it a linear scale has been used). At a broad level, colour is used to double-encode the z value into one of a number of distinct categories so you can see the “bands” of execution times.
• This graph plots execution times of a spatial query, and the audience who are interested in the results are familiar with examining representations of spatial data, including 3d terrain maps. I’d therefore argue that the presentation method of this data has been chosen with consideration of the likely end user (another thing often forgotten by people presenting data), and vaguely resembles a terrain map at the base of a mountain, for example.

It has to be said that the Excel implementation of a 3d surface chart has a few deficiencies – I’d like to have been able to set a semi-transparent fill on the surface, for example. Also, the two independent variables are treated as distinct categories and must be equally spaced along the x and y axes. If you’re plotting continuous independent variables, as I was in this case, you therefore need to ensure that your data is sampled at consistent intervals.

My next challenge, of course, is that I’m currently only considering two of the variables that effect performance of spatial queries. I now need to think about the effect of other factors: the size of the bounding box, the cells_per_object limit, and the resolution of the grid cells in the index, for example. To display the results taking these factors into consideration, I’m thinking about a panel of 3d surface charts… if any “3d chart”-deniers have any comments, or can propose a better solution, I’d love to hear from you!

An Unlikely Replacement for the Bing Maps 3D Control

So at the end of this year, on December 1st 2011, the Bing Maps 3D control is being retired. This is rather a shame – I kinda liked being able to visualise geography data from SQL Server on a round earth rather than having to project it onto a flat 2d map. And no more will I be able to enjoy this beautiful vista from the top of Mount Snowdon from the comfort of my own chair:

Fortunately, there is an alternative.

Given the recent announcement that Silverlight 5 beta will be released at Mix11, you might think that I’ve got some insider information on a new Silverlight 3d Bing Maps control, but no.

And you can do some pretty cool stuff with HTML5 and CSS 3D transforms (like this Katamari Damacy hack), but that’s not it either.

I’m thinking of something much more modest, requiring only a single cross-platform plugin that almost every user already has installed on their system….

Flash? Nope…. Adobe Acrobat.

Using some STRM topography data, an image from the Bing Maps imagery service, and a splash of Safe FME, you can create fully zoomable, rotatable terrain models in PDF format. Now, ok, this doesn’t quite provide the ability to navigate the whole world in 3D – it works best for only certain limited areas of geographic interest. But look at the advantages: you can distribute interactive 3d spatial data to anyone with Adobe Acrobat reader installed – no other software required, no special licences needed, no internet connection….

Here’s the Safe FME workflow I used (inspired largely by one of Don Murray’s demos at the Safe FME world tour last week):

And here’s a screenshot from the Adobe 3D PDF output, showing the Bing Maps imagery for Snowdonia national park in Wales overlaid on a STRM height map model. Download the PDF file and try it yourself.

There are several navigation tools in Adobe PDF reader allowing you to pan, rotate, and zoom around the model. There’s even a ‘Fly’ navigation mode that allows you to re-enact the 3D control style of Ricky’ Brundritt’s VE3D flight simulator.

And now I can once again enjoy views like this on my monitor:

(Note – there is a known bug with Adobe Reader 9.3 that might produce a “3D data parsing error has occurred” message – see http://forums.adobe.com/message/2651682 . This error can be resolved by upgrading to Adobe Reader X.)

Examining 3D Terrain of Bing Maps Tiles with SQL Server 2008 and WPF (Part 3)

In the previous post in this series, I selected a Bing Maps tile and converted its quadkey to a POLYGON representing the geographic extent of that tile. I then used that POLYGON as the basis for a query of the GTOPO data in SQL Server to retrieve the corresponding elevation data for the terrain on the tile. In this post I’ll combine these two into a 3d model showing the terrain of that tile.

Triangulating the Elevation Data

The GTOPO data is a set of distinct elevation recordings, recorded on a grid spaced at regular intervals, like this:

In order to construct a 3d terrain model from this elevation data, we need to construct a continuous smooth surface from those distinct points. For this, we will use triangulation.

Tringulation, as its name suggests, is the process of creating triangles from a set of data points. These triangles will tesselate together to cover the entire extent of data, with no gaps and no overlaps. There are many different triangulations of the same set of points – but I’ll use the Delauney Triangulation. I could probably write a whole new blog post on the Delauney triangulation, but for now all I’ll say is that there are many examples of code on the internet that describe how to create a set of triangle polygons from an input set of points. For example, you could try looking here, here, or here.

Triangulating the set of elevation points leads to a set of triangular polygons, like this:

At first glance it may not be obvious why having a set of triangles in which each vertex is a point in the dataset is any better than the set of points we started with. What you’ve got to remember (and what the spatial results tab in SQL Server Management Studio can’t show you) is that each of the vertices of these triangles has an associated Z value – they are all at different heights. Therefore the triangles in our dataset now form a set of connected, angled faces.

Creating the WPF 3D Mesh

To display a 3D mesh representing the surface constructed from these triangular faces, I’ll use a WPF MeshGeometry3D object. First, connect to the database and loop through the triangulated SqlGeometry polygons, adding the X, Y, and Z values to the Mesh Positions array:

MeshGeometry3D mesh = new MeshGeometry3D();

while (dataReader.Read())

{

  SqlGeometry v = (SqlGeometry)dataReader.GetValue(0);

  for (int n = 3; n >= 1; n—)

  {

    mesh.Positions.Add(new Point3D(

      (double)v.STPointN(n).STX,

      (double)v.STPointN(n).Z / 10000,

     -(double)v.STPointN(n).STY

    ));

  }

}

There’s a couple of points to note here:

• Firstly, every value returned by my DataReader is a triangular SqlGeometry polygon. Polygons must be closed, so they actually contain 4 points – the last point being the same as the first point. However, we only want the three distinct vertex locations.
• Secondly, when defining 3D objects in WPF (as in most 3D applications), the order in which you define the coordinates is important as this is used to determine the “direction” in which each face points. Unless you explicitly specify vertex normals, you should define vertices in anticlockwise order as you look at them to ensure that the associated face points towards you  Faces are single-sided, so if you look at them from behind they will appear invisible. I loop through the points of each triangle from n = 3 to n = 1 to ensure they are added to the mesh in anticlockwise order.
• I’ve applied an arbitrary scaling factor of 10,000 to the Z value just to make the mountains and valleys look aesthetically attractive compared to the horizontal dimensions. I wouldn’t have needed to do this if I had first projected my latitude and longitude coordinates into metres, so that they were measured in a unit consistent with the GTOPO elevation data.
• Notice that, when considering a flat two-dimensional object, we tend to think of the x axis as extending across towards the right of the screen, and the y coordinate coordinate extends up the screen. The z coordinate comes towards us from the screen. However, the WPF 3d coordinate system treats x coordinate as extending to the right, the z coordinate extending upwards, and the y coordinate extending forwards. So, to map the SqlGeometry coordinate values to the WPF Point3D object, I use X = X, Y = Z, and Z = –Y.

The resulting mesh, illustrating the mountains and valleys of Scotland contained on this tile, looks like this:

Adding Material to the Mesh

Now we’ve got our mesh, we can specify a material to paint it with. But, before we do, we need to specify texture coordinates for each point of the mesh. Texture coordinates are used to describe how a 2D image should be stretched over a 3D shape. Fortunately, because we’re dealing with a square tile, the texture coordinates are thankfully simple – we want to stretch the entire image over the entire mesh, with the image equally stretched at every point.

Thus, for each point added to the mesh, the associated texturecoordinate is as follows:

mesh.Positions.Add(new Point3D(x, z, -y));

mesh.TextureCoordinates.Add(new Point(x, -y));

Once we’ve specified the appropriate texturecoordinates for each point in the mesh, we can create a material based on an ImageBrush to superimpose the original Bing Maps tile image over the top.

ImageBrush imgBrush = new ImageBrush();

imgBrush.ImageSource = new BitmapImage(new Uri(@"http://ecn.t2.tiles.virtualearth.net/tiles/r031133.png?g=603&mkt=en-us", UriKind.Absolute));

Material mat = new DiffuseMaterial(imgBrush);

this.meshGeometry.Material = mat;

Here’s what the terrain mesh looks like with an aerial tile image brush:

Or, with a road tile:

Smoothing the Mesh

Currently, my mesh lists each point of each triangle separately – it contains 2,099 faces and 6,297 vertices (3 distinct vertices for each triangle). However, many of the vertices in the mesh are actually duplicates – any vertex at which two or more triangles meet is currently listed multiple times. The effect of listing points in this manner is to cause each triangle to be rendered individually, causing the “hard edges” between faces shown in the preceding image.

As an alternative, we could list only unique points in the Positions list of the mesh, and then use the TriangleIndices property to list the index positions of each vertex that forms a triangle. For example, when looping through each point, x:

int x = mesh.Positions.IndexOf(p);

// If this point is not already in the mesh

if (x == -1)

{

  mesh.Positions.Add(p1);

  mesh.TextureCoordinates.Add(p.X, -p.Y));

  mesh.TriangleIndices.Add(mesh.Positions.Count - 1);

}

// If it already exists

else

{

  mesh.TriangleIndices.Add(x);

}

The effect this has is to smooth the mesh, as follows:

And, when overlaid with a road tile: