January 6, 2011

## 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:

## 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:

Tags: , , ,
December 26, 2010

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

In the last post, I downloaded some elevation data from the GTOPO digital elevation model and loaded it into a SQL Server table. In this post, I’m now going to select a subset of that data in order to create a terrain map of a particular Bing Maps tile.

## The Bing Maps Tile System

As you probably already know, Bing Maps are constructed from a set of raster “tiles” – prerendered 256px x 256px gif, jpg, or png images that form the background of each map. There are different sets of tiles that correspond to the different styles of map.

Here’s a few examples:

Aerial Map Tile

Ordnance Survey Map Tile

Every tile is numbered with a quadkey – a unique identifier that exactly describes the position and zoom level of where that tile should be placed. I’m not going to explain the quadkey system here, because there’s already an excellent article on MSDN.

The important fact is that, by knowing only the quadkey of a tile, it is possible to determine the exact coordinate boundaries of every corner of that tile. In SQL Server, that means we can construct a POLYGON that represents the extent of that tile, and retrieve any data that lies on the tile (or, in this example, the elevation data that corresponds to the terrain of that tile). The following method demonstrates how to return a SqlGeography polygon from a Bing Maps quadkey:

`public static SqlGeography GeographyPolygonFromQuadKey(string quadKey) {`

`      int tilex = 0, tiley = 0;`

`      int tilewidth = 256, tileheight = 256;`

`      int zoomLevel = quadKey.Length;`

`      // Work out the x and y (grid) position of this tile`

`      for (int i = zoomLevel; i > 0; i--)`

`      {`

`        int mask = 1 << (i - 1);`

`        switch (quadKey[zoomLevel - 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.");`

`            }`

`        }`

` `

`        // From the grid position and zoom, work out the min and max Latitude / Longitude values of this tile`

`        float minLongitude = (float)(tilex * tilewidth) * 360 / (float)(tilewidth * Math.Pow(2, zoomLevel)) - 180;`

`        float maxLatitude = (float)Math.Asin((Math.Exp((0.5 - (tiley * tileheight) / (tileheight) / Math.Pow(2, zoomLevel)) * 4 * Math.PI) - 1) / (Math.Exp((0.5 - (tiley * tileheight) / 256 / Math.Pow(2, zoomLevel)) * 4 * Math.PI) + 1)) * 180 / (float)Math.PI;`

`        float maxLongitude = (float)((tilex + 1) * tilewidth) * 360 / (float)(tilewidth * Math.Pow(2, zoomLevel)) - 180;`

`        float minLatitude = (float)Math.Asin((Math.Exp((0.5 - ((tiley + 1) * tileheight) / (tileheight) / Math.Pow(2, zoomLevel)) * 4 * Math.PI) - 1) / (Math.Exp((0.5 - ((tiley + 1) * tileheight) / 256 / Math.Pow(2, zoomLevel)) * 4 * Math.PI) + 1)) * 180 / (float)Math.PI;`

` `

`        // Create a new SqlGeography instance representing the extent of this tile`

`        SqlGeographyBuilder gb = new SqlGeographyBuilder();`

`        gb.SetSrid(4326);`

`        gb.BeginGeography(OpenGisGeographyType.Polygon);`

`        gb.BeginFigure(minLongitude, maxLatitude);`

`        gb.AddLine(minLongitude, minLatitude);`

`        gb.AddLine(maxLongitude, minLatitude);`

`        gb.AddLine(maxLongitude, maxLatitude);`

`        gb.AddLine(minLongitude, maxLatitude);`

`        gb.EndFigure();`

`        gb.EndGeography();`

` `

`      return gb.ConstructedGeography;`

`    }`

Now, I just need to choose a particular tile whose terrain to examine…. I’ve chosen the tile at zoom level 6 that portrays most of Scotland – my ancestral home and also home to some interesting geographic features (i.e. mountains and lochs!). Here’s the tile viewed using Bing Maps’ road map style:

The quadkey for this tile is 031133 and, using the method above, I can calculate the corresponding SqlGeography as:

POLYGON((-5.625 55.7765730186677, 0 55.7765730186677, 0 58.8137417157078, –5.625 58.8137417157078, -5.625 55.7765730186677))

Therefore, we can retrieve the relevant GTOPO30 elevation data for this tile from SQL Server using the following query:

SELECT * FROM W020N90 WHERE geog4326.STIntersects(geography::STPolyFromText(‘POLYGON((-5.625 55.7765730186677, 0 55.7765730186677, 0 58.8137417157078, –5.625 58.8137417157078, -5.625 55.7765730186677))’, 4326) = 1;

Here’s the results:

In the next post, I’ll use this elevation data to create a 3d surface onto which to overlay the Bing Maps tile image.

Tags: , , ,
December 17, 2010

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

Following on from my last last post, I thought I’d continue to explore the topic of altitude (or elevation) and how it relates to Bing Maps.

A map is, by definition, two-dimensional – it is an image that has been projected onto a flat surface. However, it is still possible, and on many occasions very important, to portray three-dimensional properties of features on that map. One way of doing this is by presenting an oblique, “birdseye” view rather than a pure top-down view. This is certainly an effective way to present the height of individual distinct features such as buildings, as in the default Bing Maps control. But what about displaying variations in height of the underlying terrain? In this first of a set of posts I’ll look at a way of combining altitude data, SQL Server, and WPF to create a “3D” view of the elevation of terrain on a Bing Maps tile.

## Obtaining Altitude Data

Bing Maps does not expose altitude data, so the first step is to acquire this data from somewhere else. There have been various scientific surveys of global elevation data, produced by government agencies and scientific organisations such as NASA. Three common sets of data are as follows:

• GTOPO30 – global elevation grid at resolution of 30 arcseconds (approx 1km).
• SRTM – elevation grid at resolution 3 arcseconds (approx 90m) between latitudes –56 and +60.
• ASTER – elevation grid with resolution 1 arcsecond (approx 30m) between latitudes –83 and 83.

### Elevation Webservices

The data from various elevation surveys is exposed via webservices provided by the USGS, Geonames and Google (amongst others). For example, to determine the altitude of the point located at (53.06839056444848, -4.076281785964964), you could call one of the following URLs

These services are very useful if you want to retrieve the elevation of only one, or a small number of points. However, for this example I’m going to use several thousand elevation points spread out across the surface of a Bing Maps tile, so these services aren’t really suitable.

### Elevation Datasets

Each dataset is also available to download in full. I’ll be using the GTOP30 dataset, which you can download from http://eros.usgs.gov/#/Find_Data/Products_and_Data_Available/gtopo30_info. It’s split into several smaller tiles, so if you only want to download a particular area that’s fine. Tiles are named according to the coordinates at the top left corner. I’m going to use the W020N90 tile, which starts at a longitude of –20 and latitude of 90, and covers most of Western Europe, as highlighted below:

The W020N90 data is provided as a gzip archive, which is about 8.6Mb in size.

The GTOPO data is supplied as a DEM image, so the next step is to convert it into a more useable format. I’m going to use the excellent free MicroDEM program to crop a section of the W020N90 DEM file and save it as ASCII text file.

### Converting the DEM Binary file to ASCII XYZ

2. Once the data has loaded, use the Subset & zoom tool from the MicroDEM menu bar to select the part of the map to crop (if desired).
3. Select File –> Save DEM –> Current subset, MD DEM.
4. Close the current file and re-open the newly saved cropped file.
5. Save the cropped file in text format by selecting File –> Save DEM –> ASCII XYZ.

### Tidying Up the ASCII XYZ file

Once saved out as a ASCII XYZ file, the elevation data is stored as a very simple dataset with three columns as follows:

LAT        LONG        ELEVATION

52.837501  -8.479166   55
53.670834  -8.479166   77

The only slightly annoying thing about the XYZ export from MicroDEM is that there isn’t a consistent separator used to delimit the columns – it seems that any length of space is used as a separator. To fix this, I then opened up Visual Studio and performed a quick replace (Ctrl + H) using regular expressions:

• First, to remove leading spaces on some of the lines, search for ^[ \t]+ and replace with an empty string.
• Then, to replace the variable-length spaces between columns, search for [ ][ ]* and replace with ,

We’ve now got a nice comma-separated CSV file that is easy to import into SQL Server:

LAT,LONG,ELEVATION

52.837501,-8.479166,55
53.670834,-8.479166,77

### Importing into SQL Server

Once you’ve imported the CSV file into a new table in SQL Server (Use the import/export wizard), we can add a geography column to that table, and populate it with a Point in each row. Since we want our Points to contain Z values, we have to use one of the WKT static methods (the geography Point() method, the xxxFromWKB() and GeomFromGml() methods can only be used to create two-dimensional coordinates):

`ALTER TABLE W020N90`

`ADD GEOG4326 geography;`

` `

`UPDATE W020N90`

`SET GEOG4326 = geography::STPointFromText('POINT(' + LONG + ' ' + LAT + ' ' + ELEVATION + ')', 4326);`

(note that WKT requires longitude/latitude/altitude coordinate ordering, not latitude/longitude/altitude as you might expect).

At this stage, you can check the GTOPO elevation point data (well, the first 5,000 rows of it, anyway) by executing a SELECT query from the SQL Server table and checking the Spatial Results tab in SQL Server Management Studio, as shown below:

In the next post, I’ll examine how you can convert this table of points into a set of polygonal faces that can be used as the basis for a 3D terrain model for any particular geographic area.

Tags: , , ,