Archive for June, 2011

June 29, 2011

Creating a Custom Bing Maps Scale Bar (A.K.A. “Yes, but how many Space Needles is that?”)

Earlier in the year, during the Microsoft MVP summit, myself and a group of fellow delegates spent a very enjoyable Sunday afternoon sightseeing around Seattle (although February in Seattle was a bit cold for Soul Solutions’ John O’Brien, who is used to the slightly more clement weather down under!)

Whilst ascending the glass elevator during the guided tour of the Seattle Space Needle, the tour guide took pride in pointing out various objects that could be seen in the distance, and comparing their relative size in terms of “space needles”.

This then became a bit of a running joke for the duration of the summit. The Boeing factory was big, but how big was it in space needles? Likewise, the Safeco field, home of the Seattle Mariners. Bigger or smaller than a space needle?

Fortunately, you need wonder no more. The following simple (and slightly silly) code listing demonstrates how you can add your own custom scale bar to Bing Maps AJAX v7. It makes use of the getMetersPerPixel() method of the Map class to obtain the resolution in metres/pixel at the centre of the current map view. If you also know the length (or height, in this case) or any arbitrary object, you can then use this information to create your own units of map measurement.

<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html>
<head>
  <title></title>
  <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
  <style type="text/css">
    #SpaceNeedleScaleBar
    {
      white-space: nowrap;
      color: #FFFFFF;
      font-size: 80%;
      position: absolute;
      z-index: 1;
      width: 100px;
      right: 167px;
      bottom: 30px;
      height: 34px;
    }
    #SpaceNeedleScaleBarbg
    {
      position: absolute;
      top: 1px;
      right: 0;
    }
    #SpaceNeedleScaleBarfg
    {
      position: absolute;
      top: 0;
      right: 1px;
      color: #000000;
    }
    #SpaceNeedleimg
    {
      position: absolute;
      width: 80px;
      bottom: 0;
      right: 3px;
    }
  </style>
  <script type="text/javascript" src="http://ecn.dev.virtualearth.net/mapcontrol/mapcontrol.ashx?v=7.0&amp;mkt=en-gb"></script>
  <script type="text/javascript">
    var map = null;
    function GetMap() {
      // Create a basic map
      map = new Microsoft.Maps.Map(
         document.getElementById("mapDiv"),
         { credentials: "INSERTYOURBINGMAPSKEY",
           center: new Microsoft.Maps.Location(54, -3),
           zoom: 5
         });

      // Update the scale bar every time the map view changes
      Microsoft.Maps.Events.addHandler(map, 'viewchangeend', updateScaleBar);
    }

    function updateScaleBar() {
      var spaceNeedleHeight = 184; // Actual height of space needle (metres)
      var spaceNeedleImage = 60; // Size of space needle image (pixels)
      var mapResolution = map.getMetersPerPixel(); // Resolution at centre of map (metres/pixel)

      // Calculate how many space needles will fit into the length of the image at this zoom level 
      var howmanyspaceneedles = (mapResolution / spaceNeedleHeight * spaceNeedleImage).toFixed(2) + ' space needles';

      // Update the display
      document.getElementById('SpaceNeedleScaleBarbg').childNodes[0].nodeValue = howmanyspaceneedles;
      document.getElementById('SpaceNeedleScaleBarfg').childNodes[0].nodeValue = howmanyspaceneedles;
    }
  </script>

</head>
<body onload="GetMap();">
  <div id='mapDiv' style="position: relative; width: 480px; height: 640px;">
    <div id="SpaceNeedleScaleBar">
      <div id="SpaceNeedleScaleBarbg">space needles</div>
      <div id="SpaceNeedleScaleBarfg">space needles</div>
      <div id="SpaceNeedleimg"><img alt="Seattle Space Needle" src="spaceneedle.png" /></div>
    </div>
  </div>
</body>
</html>

The result, showing the new “space needles” scale alongside the more traditional imperial and metric units at the bottom of the map, looks like that below. Or you can try out an example here.

image

Tags:
June 27, 2011

Geodesics on Bing Maps v7

This month’s MSDN magazine has an article describing how to create curved lines on the Bing Maps AJAX control. While I don’t want to criticise the author at all, there are two comments I would make on the article:

  • Firstly, it’s written using v6.3 of the AJAX control – v7.0 has been available for well over 6 months now and (despite some teething problems) this latest version is recommended for all new development.
  • Secondly, the article describes how to draw arbitrary Bezier curves on the projected plane of the map. Whilst this is an interesting exercise (and the author goes on to describe important concepts such as how to test the routine), it’s not actually that useful. More often, when we see curved lines on a map, we expect them to represent geodesics – the shortest path between two points on the surface of the earth. Although this was never the intention of the article, Bing Maps evangelist Chris Pendleton mistakenly drew this conclusion and tweeted a link to the article stating that it demonstrated geodesics, when in fact it does not.

Therefore, I thought that responding to this article would provide a good prompt for me to dust off and update my own v6.x geodesic curve script from several years ago (originally published here).

What’s a Geodesic, anyway?

A geodesic is a “locally-length minimising curve” (Mathworld) – it’s the shortest path between any two points on a given surface. On a flat plane, like a piece of paper, a geodesic is a straight line. On a sphere, a geodesic is a great circle. When dealing with geospatial data, a geodesic is the shortest distance between two points on the surface of the earth.

Generally speaking, Bing Maps has no regard for geodesic shapes relative to the earth’s surface – instead it draws shapes directly onto the projected map image. Drawing a straight line between two points on a map represents the shortest path between those points in the projected plane of the map, but it generally does not represent the shortest path between those same two points on the surface of the earth.

For example, consider a polyline drawn between Munich and Seattle, both of which lie at a latitude of approximately 48 degrees. You can define a polyline connecting these two points as follows:

Microsoft.Maps.Polyline([
  new Microsoft.Maps.Location(48, 11.5),
  new Microsoft.Maps.Location(48, -122)]);

When displayed on the map, this polyline will follow a constant latitude between the two points, like this:image

However, this is certainly not the shortest route between Munich and Seattle. If you are unsure why, consider how this same line would appear when viewed on a 3-dimensional model of the earth. In the screenshot below, the line that follows a constant latitude of 48 degrees, as shown above, is plotted in red, while the geodesic line that represents the true shortest line connecting the destinations is shown in blue:

image

Notice how, rather than being parallel to the equator, the geodesic route goes much further north over the top of the UK, then over Iceland, Greenland, and much of Canada before turning south again into Seattle. You can try this yourself on a globe – the geodesic route is the path that a piece of string follows when held tight between two locations. (For those readers familiar with SQL Server, the red line above is equivalent to a route calculated using the geometry datatype, while the blue line is equivalent to using the geography datatype)

As shown above, the shortest “straight line” route on a map is not the shortest direct path between two points on a globe. Likewise, the shortest geodesic route between two locations on the globe does not generally correspond to a straight line on a map. This is why, when airline companies show maps illustrating the flightpaths to various destinations, the lines appear curved – because they’re representing the geodesic path on the surface of the earth, which, when projected onto a map, will generally not be straight lines:

image

Drawing Geodesic curves in Bing Maps

Geodesics are clearly very useful if you want to visualise the path of shortest distance between two points. So how do you go about drawing geodesic curves in Bing Maps? Well, Bing Maps does not support curved geometries directly, so instead we must approximate the shape of a geodesic curve by creating a polyline containing several small segments. Using a larger number of segments will make the polyline appear more smooth and more closely resemble the shape of the smooth curve, but will also increase its complexity. I find that using 32 segments is more than sufficient accuracy for most maps. We’ll call this value n.

var n = 32;

Then, we need to determine the overall extent of the route, which we’ll call d. The shortest distance between any two points on a sphere is the great circle distance. Assuming that the coordinates of the start and end points are (lat1, lon1) and (lat2, lon2) respectively, measured in Radians, then we can work out the great circle distance between them using the Haversine formula, as follows:

var d = 2 * asin(sqrt(pow((sin((lat1 - lat2) / 2)), 2) + cos(lat1) * cos(lat2) * pow((sin((lon1 - lon2) / 2)), 2)));

We then determine the coordinates of the endpoints of each segment along the geodesic path. If f is a value from 0 to 1, which represents the percentage of the route travelled from the start point (lat1,lon1) to the end point (lat2,lon2), then the latitude and longitude coordinates of the point that lies at f proportion of the route can be calculated as follows:

var A = sin((1 - f) * d) / sin(d);
var B = sin(f * d) / sin(d);

// Calculate 3D Cartesian coordinates of the point
var x = A * cos(lat1) * cos(lon1) + B * cos(lat2) * cos(lon2);
var y = A * cos(lat1) * sin(lon1) + B * cos(lat2) * sin(lon2);
var z = A * sin(lat1) + B * sin(lat2);

// Convert these to latitude/longitude
var lat = atan2(z, sqrt(pow(x, 2) + pow(y, 2)));
var lon = atan2(y, x);

By repeating the above with different values of f, (the number of repetitions set according to the number of segments in the line), we can construct an array of latitude and longitude coordinates at set intervals along the geodesic curve from which a polyline can be constructed.

The following code listing wraps this all together in a reusable function, ToGeodesic, that returns an array of points that approximate the geodesic path between the supplied polyline or polygon locations.

To demonstrate the use of the function, I’ve added two entity collections to the map. The simple layer acts a regular shape layer, containing polylines and polygons displayed in red. Whenever an entity is added to this collection, the entityAdded event handler fires, which converts the added entity into the equivalent geodesic shape and adds it to the geodesicShapeLayer, displayed in blue. By maintaining two layers you can continue to deal with shapes as per normal in the layer collection – the additional points needed to simulate the geodesic display of each shape only get added in the copy of the shape added to the geodesicShapeLayer. You may then, for example, choose not to display the non-geodesic layer and only use it as a method to manage shapes, while the geodesic layer is used to actually display them on the map.

<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html>
<head>
  <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&amp;mkt=en-gb"></script>
  <script type="text/javascript">
    var map = null;
    var geodesicLayer;
    var layer;
    var baseLayerOptions = {
      strokeColor: new Microsoft.Maps.Color(255, 200, 0, 0),
      fillColor: new Microsoft.Maps.Color(100, 100, 0, 0)
    };
    var geodesicLayerOptions = {
      strokeColor: new Microsoft.Maps.Color(255, 0, 0, 100),
      fillColor: new Microsoft.Maps.Color(100, 0, 0, 100)
    };

    function GetMap() {
      map = new Microsoft.Maps.Map(
         document.getElementById("mapDiv"),
         { credentials: "INSERTYOURBINGMAPSKEY",
           center: new Microsoft.Maps.Location(40, 0),
           zoom: 2
         });

      // Add a simple layer and attach event handler
      layer = new Microsoft.Maps.EntityCollection();
      map.entities.push(layer);
      Microsoft.Maps.Events.addHandler(layer, 'entityadded', entityAdded);

      // Add the geodesic layer to the map
      geodesicLayer = new Microsoft.Maps.EntityCollection();
      map.entities.push(geodesicLayer);

      // Add a polyline to the basic layer
      var Line = new Microsoft.Maps.Polyline([new Microsoft.Maps.Location(48, 11.5), new Microsoft.Maps.Location(48, -122)], baseLayerOptions);
      layer.push(Line);

      // Add a polygon to the basic layer
      var Polygon = new Microsoft.Maps.Polygon([
        new Microsoft.Maps.Location(0, 60),
        new Microsoft.Maps.Location(30, 120),
        new Microsoft.Maps.Location(60, 90),
        new Microsoft.Maps.Location(60, 30),
        new Microsoft.Maps.Location(0, 60)
      ], baseLayerOptions);
      layer.push(Polygon);
    }

    // Event handler gets called for every new shape added to the map layer
    function entityAdded(e) {
      // Create an array of locations representing the geodesic shape
      // of the shape added
      var geodesicLocs = ToGeodesic(e.entity.getLocations(), 32);

    // Create a new shape based on these locations
      switch (e.entity.toString()) {
        case '[Polyline]':
          var geodesicShape = new Microsoft.Maps.Polyline(geodesicLocs);
          break;
        case '[Polygon]':
          var geodesicShape = new Microsoft.Maps.Polygon(geodesicLocs);
          break;
      }
      // Set a few styling options
      geodesicShape.setOptions(geodesicLayerOptions);

      // Put the shape on the geodesic shape layer
      geodesicLayer.push(geodesicShape);
    }

    // Creates geodesic approximation of the lines drawn between an array
    // of points, by dividing each line into a number of segments.
    function ToGeodesic(points, n) {
      if (!n) { n = 32 }; // The number of line segments to use
      var locs = new Array();
      for (var i = 0; i < points.length - 1; i++) {
        with (Math) {
          // Convert coordinates from degrees to Radians
          var lat1 = points[i].latitude * (PI / 180);
          var lon1 = points[i].longitude * (PI / 180);
          var lat2 = points[i + 1].latitude * (PI / 180);
          var lon2 = points[i + 1].longitude * (PI / 180);
          // Calculate the total extent of the route
          var d = 2 * asin(sqrt(pow((sin((lat1 - lat2) / 2)), 2) + cos(lat1) * cos(lat2) * pow((sin((lon1 - lon2) / 2)), 2)));
          // Calculate  positions at fixed intervals along the route
          for (var k = 0; k <= n; k++) {
            var f = (k / n);
            var A = sin((1 - f) * d) / sin(d);
            var B = sin(f * d) / sin(d);
            // Obtain 3D Cartesian coordinates of each point
            var x = A * cos(lat1) * cos(lon1) + B * cos(lat2) * cos(lon2);
            var y = A * cos(lat1) * sin(lon1) + B * cos(lat2) * sin(lon2);
            var z = A * sin(lat1) + B * sin(lat2);
            // Convert these to latitude/longitude
            var lat = atan2(z, sqrt(pow(x, 2) + pow(y, 2)));
            var lon = atan2(y, x);
            // Create a Location (remember to convert back to degrees)
            var p = new Microsoft.Maps.Location(lat / (PI / 180), lon / (PI / 180));
            // Add this to the array
            locs.push(p);
          }
        }
      }
      return locs;

    }
  </script>
</head>
<body onload="GetMap();">
  <div id='mapDiv' style="position: absolute; width: 800px; height: 600px;">
  </div>
</body>
</html>

Here’s the results, showing both the flat (red) and geodesic (blue) layers of a polyline and a polygon:

image

Tags:
June 20, 2011

Adding a Touch of Style to Vector Shapes on Bing Maps

In addition to their obvious functional purpose, maps can be objects of great aesthetic beauty. On my bookshelf I have a book (The Map Book, by Peter Barber) full of beautiful illustrations of maps throughout history, such as these:

image image
image image

Perhaps it is unreasonable to compare a simple web-mapping tool such as Bing Maps to the art of a master cartographer, but one thing is certain; namely, this is not beautiful:

image

I try to convince myself that Microsoft chose these horrible default colours so that people would look at them, think “My God – those colours are awful – I must change them!”, and then look up how to set the FillColor, StrokeColor, and StrokeThickness properties of the Microsoft.Maps.Polygon object, but in practice this doesn’t always happen. Every time I see a default-coloured polygon like that above, I’m sure Marinus of Tyre must turn in his grave.

But it doesn’t have to be like this. Applying some simple style rules to change the fill and stroke properties can have a dramatic improvement on the appearance of your map. As an example, I took around 3,000 features from the OS Vector Map dataset in an area around Norwich and plotted them on Bing Maps as polylines and polygons. Every element contained an additional feature code property, which described the type of feature using the same numbering system as with the Ordnance Survey’s own dataset. For example, feature code 25710 is a motorway; code 25999 is an area of woodland, and 25301 is a single-track railway.

Then I created a set of “stylesheets” – simple arrays of different Microsoft.Maps.Color values that should be applied to features based on their type. For example:

var style1 = {
  // Building
  25014: {
    strokeColour: new Microsoft.Maps.Color(alpha, 126, 119, 98),
    strokeWidth: 1,
    fillColour: new Microsoft.Maps.Color(alpha, 232, 214, 176)
  },
  // Glasshouse
  25013: {
    strokeColour: new Microsoft.Maps.Color(alpha, 232, 214, 176),
    strokeWidth: 1,
    fillColour: new Microsoft.Maps.Color(alpha, 255, 255, 255)
  },
  // Electricity Transmission Line
  25102: {
    strokeColour: new Microsoft.Maps.Color(alpha, 158, 170, 158),
    strokeWidth: 1
  },
  ...
}

And finally, a simple function that looped through the elements on the map and applied the correct “style” to each element based on its feature code, using the setOptions() method. The results are shown below, or you can try a live demo here.

image
Default
image
Style 1
image
Style 2
image
Black and White
image
Random (very Mondrian-esque, I think!)
 

Now aren’t they much prettier?

Tags: ,
June 18, 2011

Importing Spatial Data to SQL Server with OGR2OGR–now even easier!

A few months back, I posted an article explaining how to import spatial data into SQL Server 2008 from any format supported by the OGR library (including ESRI shapefiles, GML, and TIGER data), using OGR2OGR. That article was written using OGR2OGR from v1.7 of the GDAL 1.7 library, which doesn’t support SQL Server 2008 directly, so I instead used OGR2OGR to create a CSV file containing spatial data in Well-Known Text format and then parsed that data in SQL Server using the STGeomFromText() method.

The good news is that things have become a bit easier since then, and version 1.8 of the GDAL library now has a MSSQLSpatial driver that can interface directly with geometry and geography data in SQL Server 2008.

The bad news is that most of the places that offer pre-compiled GDAL binaries for Windows have yet to update to the new version. FWTools, for example, still comes packaged only with v1.7. Likewise, the osgeo download site at http://download.osgeo.org/gdal/win32/ also lists GDAL versions up to v1.7.

So, if you want to get hold of the latest GDAL to import directly into SQL Server you’ll have to build it yourself from source, which can be downloaded from http://download.osgeo.org/gdal/gdal180.zip

Fortunately, the source has been very considerately packaged, and includes solution files that will build GDAL out-of-the-box in VS2005, 2008, and 2010. Simply load the .sln file, click build, and wait a few minutes:

Building GDAL 1.8 in VS2010

Then, if you look in the output directory (\warmerda\bld\bin, by default) you should see a lovely collection of utilities for working with spatial data – GDAL (for working with raster data), and OGR (for its vector sibling).

Here’s the output of calling ogr2ogr –formats, which retrieves the list of supported vector spatial formats – note the MSSQLSpatial format supported for both read/write:

image 

Example usage to load a shapefile to SQL Server as follows:

ogr2ogr -overwrite -f MSSQLSpatial "MSSQL:server=.\MSSQLSERVER2008;database=spatial;trusted_connection=yes" "TG20.shp"

And here’s some lovely OS vector shape data in SQL Server:

image

Follow

Get every new post delivered to your Inbox.

Join 53 other followers