Click here to Skip to main content
15,867,308 members
Articles / Programming Languages / Visual Basic

Georeferencing Open Street Map Tiles to Use with MapWinGIS

Rate me:
Please Sign up or sign in to vote.
5.00/5 (5 votes)
19 Sep 2012CPOL3 min read 55.5K   3.5K   15   9
How to add an Open Street Map overlay to a Map using MapWinGIS.

Introduction 

One of the problems viewing GIS data, e.g. Shape Files, is that it is hard to navigate around them as there are often no common landmarks. An obvious solution to this is to add a transparent map as an overlay, referenced to the same coordinate scheme. The Open Street Map data is available, for free, as rendered tiles in a number of resolutions, covering virtually the whole planet. MapWinGIS supports overlaying image files onto map as long as they are geo-referenced, so it is relatively straightforward to add a dynamic map overlay which scales with the MapWinGIS window, to maintain a OpenStreetMap overlay at all resolutions:

Image 1

Background 

Details on Open Street Map Tiles (Slippy Maps) can be found here: http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Tile_servers.

Details on World Files, used to Geo-reference image files, can be found here: http://en.wikipedia.org/wiki/World_file.

Locating the Slippy Tiles

The two functions below, from the Slippy Tile Wiki, convert between Slippy Map Tile references and WGS-84 Latitude and Longitude. These are used to work out the tiles required to cover the map area and then to Geo-reference each tile before adding it to the Map.

Details of the method and how Slippy Tiles are numbered can be found here: http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Derivation_of_tile_names  

VB
Private Function CalcTileXY(ByVal LatDeg As Single, _
                            ByVal LonDeg As Single, _
                            ByVal ZoomLevel As Long) As PointF
 
 CalcTileXY.X = CLng(Math.Floor((LonDeg + 180) / 360 * 2 ^ ZoomLevel))
 CalcTileXY.Y = CLng(Math.Floor((1 - Math.Log(Math.Tan(LatDeg * Math.PI / 180)+ 1
                                             / Math.Cos(LatDeg * Math.PI / 180)) 
                                             / Math.PI) / 2 * 2 ^ ZoomLevel))
End Function
 
   
Private Function CalcTileTLCoord(ByVal TileCoord As PointF, _
                                 ByVal ZoomLevel As Long) As PointF
    Dim n As Double = Math.PI - ((2.0 * Math.PI * TileCoord.Y) / Math.Pow(2.0, ZoomLevel))
    CalcTileTLCoord.X = (TileCoord.X / Math.Pow(2.0, ZoomLevel) * 360) - 180
    CalcTileTLCoord.Y = 180 / Math.PI * Math.Atan(Math.Sinh(n))
 
End Function

Creating a Map Overlay

In order to create a map overlay, we need to work out the area to cover. This can be found by looking at the Extents of the AxMapWinGIS Map: 

VB
Dim MapExtents As MapWinGIS.Extents = MapMain.Extents
Dim LonMin As Double = MapExtents.xMin
Dim LonMax As Double = MapExtents.xMax
Dim LatMin As Double = MapExtents.yMin
Dim LatMax As Double = MapExtents.yMax

NB This method will only work if the AxMapWinGIS Map is using the same coordinate system as OSM Slippy Tiles i.e., WGS-84. If not, you will have to convert WGS-84 into the Map's coordinate system and vica-versa.   

Once we have the area define by Min and Max Lat and Lon, we can use one of the Slippy Map functions to convert from the coordinates to Tile references for the Top Left corner and Bottom right corner of the Map:  

VB
TLTile = CalcTileXY(LatMax, LonMin, ZoomLevel)      'Get the Tile reference for the Top Left Tile
BRTile = CalcTileXY(LatMin, LonMax, ZoomLevel)      'Get the Tile reference for the Bottom Right Tile

Then it is a simple process to loop through all the tiles in between the two corners, to create a mosaic of map tiles which cover the whole map area: 

VB
' Now add each tile in turn to the map
Dim H As Integer, V As Integer          ' Loop parameters
Dim TileRef As Point
 
For H = TLTile.X To BRTile.X            ' H is looping through Longitude
    For V = TLTile.Y To BRTile.Y        ' V is looping through Latitude
         TileRef.X = H
         TileRef.Y = V
         Call LoadOSMTile(MapMain, TileRef, ZoomLevel)
     Next V
Next H

Downloading and Geo-referencing a Tile

Each tile is stored, on the Slippy Tile Server, in a directory structure: 

VB
Const ServerName As String = "http://a.tile.openstreetmap.org/"
Dim URL As String = ServerName & ZoomLevel & "/" & TileRef.X & "/" & TileRef.Y & ".png"

To download a tile, we use Net.HttpWebRequest:

VB
'Get the image from the web server
Dim request As Net.HttpWebRequest = DirectCast(Net.HttpWebRequest.Create(URL), Net.HttpWebRequest)
Dim response As Net.HttpWebResponse = DirectCast(request.GetResponse, Net.HttpWebResponse)
Dim Image As System.Drawing.Image = Image.FromStream(response.GetResponseStream())
response.Close()

Once we have downloaded the tile, we save it in the local cache:

VB.NET
'Save the image in the local cache
Image.Save(CacheTileDirectory & "\" & CacheFileName)

And Geo-reference it by creating a World File for the tile (using the .pgw extension for a World file associated with a .png file):

VB
' We georeference the image using the Top Left Coordinate 
Dim TileTLCoord As PointF = CalcTileTLCoord(TileRef, ZoomLevel)

' We also need the width and height of the image in the map coordinates (degrees)
' We get this by subtracting the difference between this tile and the next tile
Dim Point2 As PointF
Point2.X = TileRef.X + 1
Point2.Y = TileRef.Y + 1
Dim TileWidthDeg As Double = CalcTileTLCoord(Point2, ZoomLevel).X - CalcTileTLCoord(TileRef, ZoomLevel).X
Dim TileHeightDeg As Double = CalcTileTLCoord(TileRef, ZoomLevel).Y - CalcTileTLCoord(Point2, ZoomLevel).Y
Debug.Print("Tile Width = " & TileWidthDeg & ", Height = " & TileHeightDeg)
 
'Create World file to geo-reference the PNG. See: http://en.wikipedia.org/wiki/World_file
Dim oWrite As System.IO.StreamWriter = IO.File.CreateText(CacheTileDirectory & _
    "\" & Left(CacheFileName, CacheFileName.Length - 4) & ".pgw")
oWrite.WriteLine(FormatNumber(TileWidthDeg / 256, 10))           'width of each pixel in map units (degrees)
oWrite.WriteLine("0.0")                                          'rotation parameter (not used in this case)
oWrite.WriteLine("0.0")                                          'rotation parameter (not used in this case)
oWrite.WriteLine(FormatNumber(-1 * TileHeightDeg / 256, 10))     'height of each pixel in map units (degrees)
oWrite.WriteLine(FormatNumber(TileTLCoord.X, 10))                'x coordinate of the Upper LH cell (degrees)
oWrite.WriteLine(FormatNumber(TileTLCoord.Y, 10))                'y coordinate of the Upper LH cell (degrees)
oWrite.Close()

Note that the Height is negative, to indicate pixels going down vertically from the top left corner (the geo-reference point).

Once we have saved the World File to the same directory as the PNG image, we can add the PNG file to the Map. MapWinGIS will automatically detect the world file, as long as it is in the same directory:

VB
'Load the geo-referenced PNG Slippy Tile (MapWinGIS auto detects 
the World file and loads that at the smae time)
Dim ImageType As MapWinGIS.ImageType = MapWinGIS.ImageType.PNG_FILE
Dim MWImage As New MapWinGIS.Image
If Not MWImage.Open(CacheTileDirectory & "\" & CacheFileName, ImageType) Then
     Debug.Print("Image load failed: " & MWImage.ErrorMsg(MWImage.LastErrorCode))
     Exit Sub
End If
 
'Set Transparency to 50%
MWImage.TransparencyPercent = 0.5
 
'Add to the Map as a new layer
Dim PNGLayer As Integer = MapMain.AddLayer(MWImage, True)

Using the code    

The complete Visual Studio 2010 Project is included which comprises the code above and a Windows Form wrapper which lets the user navigate a Shape file, zooming in and zooming out. The mouse position is displayed in a text window at the bottom of the screen, so you can check the overlay is correct by locating known features etc:

Image 2 

The only external reference required is a Shape File to load, and the code currently uses one of the MapWindow example project files, which installs with MapWindow.

"C:\Program Files\MapWindow\Sample Projects\World\Shapefiles\WORLD30.shp"

So, either make sure this is present, or choose another file to play with.

History

First version 19 Sep 2012.

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)


Written By
Architect CBNL
United Kingdom United Kingdom
I dabble in a bit of programming as part of my job.

Comments and Discussions

 
QuestionOn zoom under 8 the formulas result in a deviation in latitude Pin
Member 1143892121-May-16 0:41
Member 1143892121-May-16 0:41 
QuestionWrong Projection... Pin
AndyW9997-Oct-13 2:02
AndyW9997-Oct-13 2:02 
QuestionEPSG 3857 Pin
Member 978460228-Jan-13 21:28
Member 978460228-Jan-13 21:28 
QuestionOSGB36 Projection Pin
tvtiger28-Jan-13 21:02
tvtiger28-Jan-13 21:02 
QuestionHelp me please? Pin
i0019-Sep-12 17:06
i0019-Sep-12 17:06 
AnswerRe: Help me please? Pin
Ben Freeman CBNL19-Sep-12 20:07
Ben Freeman CBNL19-Sep-12 20:07 
GeneralRe: Help me please? Pin
cooljd26-Sep-12 17:46
cooljd26-Sep-12 17:46 
GeneralRe: Help me please? Pin
richardw4827-Sep-12 7:52
richardw4827-Sep-12 7:52 
AnswerRe: Help me please? Pin
John Fairbanks6-Dec-12 16:55
John Fairbanks6-Dec-12 16:55 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.