Do you work with geospatial datasets in a database? Are you interested in thematic mapping? Do you consider adding a spatial component to your risk analysis application? If any of this sounds interesting to you than you should join us for the TechNet Seminar in Bad Homburg, Germany.
During the day we will dive into SQL Server 2008’s spatial capabilities and the data visualization with Bing Maps.
After an introduction into the spatial data types, spatial indexing and spatial functions in SQL Server 2008 we will look into the complete spatial ETL process including nasty things like coordinate system conversions or re-projections. Once we have the data in SQL Server we can already start with interesting analysis but more often we would want to get them out and display the data on our Bing Maps. We will of course do that as well and throughout the day we will look at the do’s and don’ts.
If you are interested register here.
Daimler AG has been marketing the smart fortwo electric drive since the middle of December. When the battery is charged, the car opens a connection to the internet and starts to communicate with Windows Azure. The application was developed jointly by Daimler and Microsoft and it allows the driver to continually follow the charging process of the battery, to display the current location of the vehicle and to find the nearest charging stations on an interactive Bing Map.
The information is accessible everywhere from a desktop browser, or alternatively from a Smartphone when you’re on the move. Smart smart :-)
…is live and what an update it is. 6.5 million square kilometres of ortho-photos over the following countries:
Additionally 110,000 square kilometres of oblique aerial images
Among the many new cities is Trier with the roman Porta Nigra:
If you want to find out the exact regions of the updates have a look here. For an animated tour visited the Bing Maps World Tour and as always check out Chris Pendleton’s blog.
We are going to overlay the OpenStreetMap-Tiles as a tile-layer on top of Bing Maps. So first of all we have to place the OpenStreetMap-tiles in a virtual directory on our web server or upload them into the cloud, e.g. the Windows Azure blob storage.
Once this is done we are faced with the different naming-conventions that Bing Maps and OpenStreetMap have for their tiles. Both use the Mercator-projection and both have tiles with a size of 256x256 pixels. That’s a god start but we need to have a closer look at the naming convention. Bing Maps uses a quadkey-naming convention and assumes that all tiles are in the same directory (for more details see this article on the Bing Maps Tile System).
OpenStreetMap on the other site uses a naming convention where the URI contains subdirectories for the zoom-level which in turn contains subdirectories for the tile-x coordinate. The tile-y coordinate is then the name of the image-file, i.e. http://mytileserver/tile-x-coordinate/tile-y-coordinate.png. An example could be http://hannesve.blob.core.windows.net/osm/17/82787/72482.png where 17 is the zoom-level, 82787 is the tile-x coordinate and 72482 is the tile-y coordinate.
The good news is that all the formulas and even the complete C#-code which gets you from the tile-x/y coordinates to the quadkey and vice versa is well documented in this article on the Bing Maps Tile System.
The process of actually integrating the tile-layer is quite different for the Bing Maps Silverlight and the Bing Maps AJAX control.
The first starting points for developers who work with the Bing Maps APIs is usually the Bing Maps Developer site. In order to work with the Silverlight control you will need a Bing Maps account and you can sign-up for such an account with your Windows Live ID here. The interactive SDK provides you with the download-links for the control and samples that get you started easily.
Imports Microsoft.Maps.MapControl Public Class OpenStreetMapTileSource Inherits LocationRectTileSource Public Sub New() MyBase.new("http://az1923.vo.msecnd.net/osm/{2}/{0}/{1}.png", _ New LocationRect(New Location(-19.047, 47.379), New Location(-18.733, 47.645)), _ New Range(Of Double)(1, 17)) End Sub Public Overloads Overrides Function GetUri(ByVal x As Integer, ByVal y As Integer, _ ByVal zoomLevel As Integer) As System.Uri Return New Uri(String.Format(Me.UriFormat, x, y, zoomLevel)) End Function End Class
MainPage.ThisMap.Children.Add(MainPage.tlOSM) MainPage.tlOSM.TileSources.Add(New OpenStreetMapTileSource) MainPage.ThisMap.MapForeground.Copyright.Attributions.Clear() Dim myCopyright As New AttributionInfo myCopyright.Text = ChrW(&HA9) + " 2010 Microsoft Corporation - Map data CCBYSA 2010 OpenStreetMap.org contributors" MainPage.ThisMap.MapForeground.Copyright.Attributions.Add(myCopyright)
In the Bing Maps AJAX control we don’t have the luxury to derive from classes and override GetUri-methods. However, we can leverage the functions mentioned in this article on the Bing Maps Tile System [12] in a proxy that calculates the tile-x and-y coordinates as well as the zoom-level from the quadkey and redirects the tile-requests. In order to do that we create a generic handler which works as a proxy:
<%@ WebHandler Language="VB" Class="OsmTileServer" %> Imports System Public Class OsmTileServer : Implements IHttpHandler 'Public baseUrl As String = "http://localhost/osm/" 'Local Public baseUrl As String = "http://az1923.vo.msecnd.net/osm/" 'Azure Public quadkey As String Public lvl As Integer Public tileX As Integer Public tileY As Integer Public Sub ProcessRequest(ByVal context As HttpContext) Implements IHttpHandler.ProcessRequest 'Fetch URL-Parameters quadkey = context.Request.Params("quadkey") 'Determine Zoom-Level lvl = quadkey.Length 'Get TileXY-Coordinates QuadKeyToTileXY(quadkey, tileX, tileY, lvl) 'Build URL to OSM-tile context.Response.Redirect(baseurl + lvl.ToString + "/" + tileX.ToString + "/" + tileY.ToString + ".png") End Sub Public ReadOnly Property IsReusable() As Boolean Implements IHttpHandler.IsReusable Get Return False End Get End Property #Region "Helper Functions" Private Const EarthRadius As Double = 6378137 Private Const MinLatitude As Double = -85.05112878 Private Const MaxLatitude As Double = 85.05112878 Private Const MinLongitude As Double = -180 Private Const MaxLongitude As Double = 180 ''' <summary> ''' Clips a number to the specified minimum and maximum values. ''' </summary> ''' <param name="n">The number to clip.</param> ''' <param name="minValue">Minimum allowable value.</param> ''' <param name="maxValue">Maximum allowable value.</param> ''' <returns>The clipped value.</returns> Private Function Clip(ByVal n As Double, ByVal minValue As Double, ByVal maxValue As Double) As Double Return Math.Min(Math.Max(n, minValue), maxValue) End Function ''' <summary> ''' Determines the map width and height (in pixels) at a specified level ''' of detail. ''' </summary> ''' <param name="levelOfDetail">Level of detail, from 1 (lowest detail) ''' to 23 (highest detail).</param> ''' <returns>The map width and height in pixels.</returns> Public Function MapSize(ByVal levelOfDetail As Integer) As UInteger Return CUInt(256) << levelOfDetail End Function ''' <summary> ''' Determines the ground resolution (in meters per pixel) at a specified ''' latitude and level of detail. ''' </summary> ''' <param name="latitude">Latitude (in degrees) at which to measure the ''' ground resolution.</param> ''' <param name="levelOfDetail">Level of detail, from 1 (lowest detail) ''' to 23 (highest detail).</param> ''' <returns>The ground resolution, in meters per pixel.</returns> Public Function GroundResolution(ByVal latitude As Double, ByVal levelOfDetail As Integer) As Double latitude = Clip(latitude, MinLatitude, MaxLatitude) Return Math.Cos(latitude * Math.PI / 180) * 2 * Math.PI * EarthRadius / MapSize(levelOfDetail) End Function ''' <summary> ''' Determines the map scale at a specified latitude, level of detail, ''' and screen resolution. ''' </summary> ''' <param name="latitude">Latitude (in degrees) at which to measure the ''' map scale.</param> ''' <param name="levelOfDetail">Level of detail, from 1 (lowest detail) ''' to 23 (highest detail).</param> ''' <param name="screenDpi">Resolution of the screen, in dots per inch.</param> ''' <returns>The map scale, expressed as the denominator N of the ratio 1 : N.</returns> Public Function MapScale(ByVal latitude As Double, ByVal levelOfDetail As Integer, ByVal screenDpi As Integer) As Double Return GroundResolution(latitude, levelOfDetail) * screenDpi / 0.0254 End Function ''' <summary> ''' Converts a point from latitude/longitude WGS-84 coordinates (in degrees) ''' into pixel XY coordinates at a specified level of detail. ''' </summary> ''' <param name="latitude">Latitude of the point, in degrees.</param> ''' <param name="longitude">Longitude of the point, in degrees.</param> ''' <param name="levelOfDetail">Level of detail, from 1 (lowest detail) ''' to 23 (highest detail).</param> ''' <param name="pixelX">Output parameter receiving the X coordinate in pixels.</param> ''' <param name="pixelY">Output parameter receiving the Y coordinate in pixels.</param> Public Sub LatLongToPixelXY(ByVal latitude As Double, ByVal longitude As Double, ByVal levelOfDetail As Integer, ByRef pixelX As Integer, ByRef pixelY As Integer) latitude = Clip(latitude, MinLatitude, MaxLatitude) longitude = Clip(longitude, MinLongitude, MaxLongitude) Dim x As Double = (longitude + 180) / 360 Dim sinLatitude As Double = Math.Sin(latitude * Math.PI / 180) Dim y As Double = 0.5 - Math.Log((1 + sinLatitude) / (1 - sinLatitude)) / (4 * Math.PI) Dim mapSize__1 As UInteger = MapSize(levelOfDetail) pixelX = CInt(Clip(x * mapSize__1 + 0.5, 0, mapSize__1 - 1)) pixelY = CInt(Clip(y * mapSize__1 + 0.5, 0, mapSize__1 - 1)) End Sub ''' <summary> ''' Converts a pixel from pixel XY coordinates at a specified level of detail ''' into latitude/longitude WGS-84 coordinates (in degrees). ''' </summary> ''' <param name="pixelX">X coordinate of the point, in pixels.</param> ''' <param name="pixelY">Y coordinates of the point, in pixels.</param> ''' <param name="levelOfDetail">Level of detail, from 1 (lowest detail) ''' to 23 (highest detail).</param> ''' <param name="latitude">Output parameter receiving the latitude in degrees.</param> ''' <param name="longitude">Output parameter receiving the longitude in degrees.</param> Public Sub PixelXYToLatLong(ByVal pixelX As Integer, ByVal pixelY As Integer, ByVal levelOfDetail As Integer, ByRef latitude As Double, ByRef longitude As Double) Dim mapSize__1 As Double = MapSize(levelOfDetail) Dim x As Double = (Clip(pixelX, 0, mapSize__1 - 1) / mapSize__1) - 0.5 Dim y As Double = 0.5 - (Clip(pixelY, 0, mapSize__1 - 1) / mapSize__1) latitude = 90 - 360 * Math.Atan(Math.Exp(-y * 2 * Math.PI)) / Math.PI longitude = 360 * x End Sub ''' <summary> ''' Converts pixel XY coordinates into tile XY coordinates of the tile containing ''' the specified pixel. ''' </summary> ''' <param name="pixelX">Pixel X coordinate.</param> ''' <param name="pixelY">Pixel Y coordinate.</param> ''' <param name="tileX">Output parameter receiving the tile X coordinate.</param> ''' <param name="tileY">Output parameter receiving the tile Y coordinate.</param> Public Sub PixelXYToTileXY(ByVal pixelX As Integer, ByVal pixelY As Integer, ByRef tileX As Integer, ByRef tileY As Integer) tileX = Math.Floor(pixelX / 256) tileY = Math.Floor(pixelY / 256) End Sub ''' <summary> ''' Converts tile XY coordinates into pixel XY coordinates of the upper-left pixel ''' of the specified tile. ''' </summary> ''' <param name="tileX">Tile X coordinate.</param> ''' <param name="tileY">Tile Y coordinate.</param> ''' <param name="pixelX">Output parameter receiving the pixel X coordinate.</param> ''' <param name="pixelY">Output parameter receiving the pixel Y coordinate.</param> Public Sub TileXYToPixelXY(ByVal tileX As Integer, ByVal tileY As Integer, ByRef pixelX As Integer, ByRef pixelY As Integer) pixelX = tileX * 256 pixelY = tileY * 256 End Sub ''' <summary> ''' Converts tile XY coordinates into a QuadKey at a specified level of detail. ''' </summary> ''' <param name="tileX">Tile X coordinate.</param> ''' <param name="tileY">Tile Y coordinate.</param> ''' <param name="levelOfDetail">Level of detail, from 1 (lowest detail) ''' to 23 (highest detail).</param> ''' <returns>A string containing the QuadKey.</returns> Public Function TileXYToQuadKey(ByVal tileX As Integer, ByVal tileY As Integer, ByVal levelOfDetail As Integer) As String Dim quadKey As New StringBuilder() For i As Integer = levelOfDetail To 1 Step -1 Dim digit As Integer = 0 Dim mask As Integer = 1 << (i - 1) If (tileX And mask) <> 0 Then digit += 1 End If If (tileY And mask) <> 0 Then digit += 1 digit += 1 End If quadKey.Append(digit.ToString) Next Return quadKey.ToString() End Function ''' <summary> ''' Converts a QuadKey into tile XY coordinates. ''' </summary> ''' <param name="quadKey">QuadKey of the tile.</param> ''' <param name="tileX">Output parameter receiving the tile X coordinate.</param> ''' <param name="tileY">Output parameter receiving the tile Y coordinate.</param> ''' <param name="levelOfDetail">Output parameter receiving the level of detail.</param> Public Sub QuadKeyToTileXY(ByVal quadKey As String, ByRef tileX As Integer, ByRef tileY As Integer, ByRef levelOfDetail As Integer) tileX = InlineAssignHelper(tileY, 0) levelOfDetail = quadKey.Length For i As Integer = levelOfDetail To 1 Step -1 Dim mask As Integer = 1 << (i - 1) Select Case quadKey(levelOfDetail - i) Case "0"c Exit Select Case "1"c tileX = tileX Or mask Exit Select Case "2"c tileY = tileY Or mask Exit Select Case "3"c tileX = tileX Or mask tileY = tileY Or mask Exit Select Case Else Throw New ArgumentException("Invalid QuadKey digit sequence.") End Select Next End Sub Private Shared Function InlineAssignHelper(Of T)(ByRef target As T, ByVal value As T) As T target = value Return value End Function #End Region End Class
In our JavaScript-code we point now simply to this proxy rather than to a virtual directory:
document.getElementById("divCopyright").style.visibility = 'visible'; var tileSourceSpec = new VETileSourceSpecification('OSM', './OsmTileServer.ashx?quadkey=%4'); tileSourceSpec.MinZoomLevel = 6; tileSourceSpec.MaxZoomLevel = 17; tileSourceSpec.Opacity = 1; tileSourceSpec.ZIndex = 100; map.AddTileLayer(tileSourceSpec);
And that’s already it again. You’ll find the source code here and a live sample on Windows Azure here.
You can download the full article as PDF here.
[1] PostgreSQL / PostGIS – Download from EnterpriseDB
[2] OpenStreetMap – Explore the Maps
[3] OpenStreetMap - Wiki
[4] OpenStreetMap – Data Downloads
[5] JOSM (Java OpenStreetMap Editor)
[6] osm2pgsql – an OpenStreetMap data to PostgreSQL converter and loader
[7] Python 2.5 – Mapnik doesn’t run with newer Python versions
[8] Mapnik – for OpenStreetMap
[9] Mapnik - Download
[10] 7-Zip
[11] Slik SVN
[12] Bing Maps Tile System
[13] Bing Maps for Developers
Bing Maps has a very good coverage with roadmaps and aerial imagery. However, there are regions where our data providers have gaps. Crowed-sourcing of geospatial data might be an option to fill these gaps and one of the most active communities is around OpenStreetMap. The OpenStreetMap data is available under the Create Common Attribution-Share Alike 2.0 License which means that as long as you provide the correct attribution to the source you should be fine to integrate the data with Bing Maps.
While it is theoretically possible to link directly to the OpenStreetMap tile servers, users who have a high load are requested to run their own tile servers. This also allows you to take control over availability and scalability of the tile servers.
That’s already all we have to do here.
You can get help on osm2pgsql by typing osm2pgsql –h in the command prompt. The parameters we used are:
-m | Store data in spherical mercator -s | Store temporary data in the database. This greatly reduces the RAM usage but is much slower.-d | The name of the PostgreSQL database to connect to-U | Postgresql user name-W | Force password promptmadagascar.osm.bz2 is the name of the OpenStreetMap data file that we want to load.
Extract all 3 archives into a world_boundaries-folder into the Mapnik-directory:
#-------------------------------------------------------------------------## Change the following for different bounding boxes and zoom levels## Antanarivo, Madagascar#minZoom = 6maxZoom = 17bbox = (47.382, -19.042, 47.638, -18.735)render_tiles(bbox, mapfile, tile_dir, minZoom, maxZoom)
set HOME=C:\mapnik_0_7_0set MAPNIK_MAP_FILE=osm-local.xmlset MAPNIK_TILE_DIR=C:\tmp\TileCache\OSMgenerate_tiles.py
More on my Blog