Currently the VE Imagery service has the ability to display 26 different icons on the map. Oddly enough the default pushpin icon from the JavaScript version of Virtual Earth is not one of these icons. Unfortunately there is currently no built in method to display custom icons on the map image. There is also currently a limitation of 10 pushpins being passed to the imagery service. This article shows how to add any number of custom icons to the map image in the correct location. The default pushpin icon from the JavaScript version of Virtual Earth will be used as the custom icon. The high level steps involved to accomplish this are as follows:

1. Calculate the best map view for an array of Location objects.
2. Request map image from the imagery service that matches the best map view calculated in step 1.
3. Calculate the pixel coordinate of each location and display custom icon on map image.

Calculate the Best Map View

The best map view is the zoom level and center point of a map image that contains all locations with the closest zoom level. The best map view will be stored using a structure.

/// <summary>
/// Structure to define a BestView object which defines a map view (centerpoint and zoom level)
/// </summary>
struct BestView
{
    public Location center;
    public int zoom;

    //Constructor
    public BestView(Location _center, int _zoom)
    {
        center = _center;
        zoom = _zoom;
    }
}

Listing 1Structure to define a BestView object

The following method calculates the best map view for an array of locations. This results in a map image similar to what is displayed if the VEMap.SetMapView function is used in the JavaScript version of Virtual Earth.

/// <summary>
/// Calculates best map view for an array of points. Similar to VEMap.SetMapView method
/// </summary>
/// <param name="points">Array of Location objects</param>
/// <returns></returns>
private BestView CalculateMapView(Location[] points)
{
    //Calculate bounding rectangle
    double maxLat = -90, minLat = 90, maxLon = -180, minLon = 180;

    //default zoom scales in km/pixel from http://msdn2.microsoft.com/en-us/library/aa940990.aspx
    double[] defaultScales = { 78.27152, 39.13576, 19.56788, 9.78394, 4.89197, 2.44598, 1.22299, 
0.61150, 0.30575, 0.15287, .07644, 0.03822, 0.01911, 0.00955, 0.00478, 0.00239, 0.00119, 0.0006, 0.0003 };

    //Calculate bounding box for array of locations
    for (int i = 0; i < points.Length; i++)
    {
        if (points[i].Latitude > maxLat)
            maxLat = points[i].Latitude;

        if (points[i].Latitude < minLat)
            minLat = points[i].Latitude;

        if (points[i].Longitude > maxLon)
            maxLon = points[i].Longitude;

        if (points[i].Longitude < minLon)
            minLon = points[i].Longitude;
    }

    //calculate center coordinate of bounding box
    double centerLat = (maxLat + minLat) / 2;
    double centerLong = (maxLon + minLon) / 2;

    //create a Location object for the center point
    Location centerPoint = new Location();
    centerPoint.Latitude = centerLat;
    centerPoint.Longitude = centerLong;

    //want to calculate the distance in km along the center latitude between the two longitudes
    double meanDistanceX = HaversineDistance(centerLat, minLon, centerLat, maxLon);

    //want to calculate the distance in km along the center longitude between the two latitudes
    double meanDistanceY = HaversineDistance(maxLat, centerLong, minLat, centerLong) * 2;

    //calculates the x and y scales
    var meanScaleValueX = meanDistanceX / mapWidth;
    var meanScaleValueY = meanDistanceY / mapHeight;

    double meanScale;

    //gets the largest scale value to work with
    if (meanScaleValueX > meanScaleValueY)
        meanScale = meanScaleValueX;
    else
        meanScale = meanScaleValueY;

    //intialize zoom level variable
    int zoom = 1;

    //calculate zoom level
    for (var i = 1; i < 19; i++)
    {
        if (meanScale >= defaultScales[i])
        {
            zoom = i;
            break;
        }
    }

    //return a BestView object with the center point and zoom level to use
    return new BestView(centerPoint, zoom);
}

Listing 2 Method to calculate the best map view for an array of Locations

The method that calculates the best map view uses a method called HaversineDistance. This method calculates the distance in kilometers between two coordinates using the Haversine formula. The code for this method is as follows:

/// <summary>
/// Calculate the distance in kilometers between two coordinates
/// </summary>
/// <param name="lat1"></param>
/// <param name="lon1"></param>
/// <param name="lat2"></param>
/// <param name="lon2"></param>
/// <returns></returns>
private double HaversineDistance(double lat1, double lon1, double lat2, double lon2)
{
    double earthRadius = 6371;
    double factor = Math.PI / 180;
    double dLat = (lat2 - lat1) * factor;
    double dLon = (lon2 - lon1) * factor;
    double a = Math.Sin(dLat / 2) * Math.Sin(dLat / 2) + Math.Cos(lat1 * factor) * Math.Cos(lat2 * factor) 
* Math.Sin(dLon / 2) * Math.Sin(dLon / 2);

    double c = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));

    return earthRadius * c;
}

Listing 3 Haversine distance method

Request Map Image from the Imagery Service

By specifying a specific zoom level and center point for a map, this ensures that we have at least one point of reference between pixel coordinates and a latitude/longitude coordinate. The following method retrieves an image URL of a map for a specific map view. The image is then retrieves using a stream.

/// <summary>
/// Gets map from VE Imagery Web Service
/// </summary>
/// <param name="view">Map view</param>
/// <returns></returns>
private Image GetMapImage(BestView view)
{
    //Map uri request object is created
    MapUriRequest mapUriRequest = new MapUriRequest();

    // Credentials are set using a valid Virtual Earth Token
    mapUriRequest.Credentials = new VEImageryService.Credentials();
    mapUriRequest.Credentials.Token = clientToken;

    //Pass centerpoint of map
    mapUriRequest.Center = view.center;

    // Map style and zoom level are set
    MapUriOptions mapUriOptions = new MapUriOptions();
    mapUriOptions.Style = MapStyle.Road;

    //Set zoom level of map
    mapUriOptions.ZoomLevel = view.zoom;

    // Size of the requested image in pixels is set
    mapUriOptions.ImageSize = new VEImageryService.SizeOfint();
    mapUriOptions.ImageSize.Height = mapHeight;
    mapUriOptions.ImageSize.Width = mapWidth;

    //Map options are added to the request
    mapUriRequest.Options = mapUriOptions;

    //ImageryServiceClient object created
    ImageryServiceClient imageryService = new ImageryServiceClient();

    //A request for a map uri is made
    MapUriResponse mapUriResponse = imageryService.GetMapUri(mapUriRequest);

    //Get image from URI
    System.Net.WebRequest request = 
System.Net.WebRequest.Create(mapUriResponse.Uri.Replace("{token}", clientToken));

    System.Net.WebResponse response = request.GetResponse();
    System.IO.Stream responseStream = response.GetResponseStream();

    //return image
    return Image.FromStream(responseStream);
}

Listing 4 Method to retrieve map image from the VE Web Services for best map view

Calculate Pixel Coordinate for a Location

Virtual Earth tiles are positioned using quad key tile positioning. In Virtual Earth a map tile is 256 pixels by 256 pixels. For zoom level one, there are 4 map tiles displayed, making the total map width and height to display the whole world 512 pixels by 512 pixels. If the zoom level is 2, there are 16 map tiles being displayed, making the total map width and height to display the whole world 1024 pixels by 1024 pixels. Additional information the Virtual Earth tiling system can be found here: http://msdn.microsoft.com/en-us/library/bb259689.aspx

With the Virtual Earth tiling system in mind, all pixel coordinates for a set of latitude/longitude coordinates can be calculated. Using the center point of the map image the top left hand corner of the map image that is returned by the Virtual Earth imagery service can be calculated as a pixel location on a map that displays the whole world. Relative pixel coordinates can be then calculated for all location latitude/longitude coordinates.

The following method takes in a Location object and a BestView object and returns a Point object which refers to the pixel location of a latitude/longitude coordinate on the map image.

/// <summary>
/// Converts a latlitude longitude coordinate to a pixel coordinate of a map based on the map view
/// </summary>
/// <param name="latlong"></param>
/// <param name="view"></param>
/// <returns></returns>
private Point LatLongToPixel(Location latlong, BestView view)
{
    //Formulas based on following article:
    //http://msdn.microsoft.com/en-us/library/bb259689.aspx

    //calcuate pixel coordinates of center point of map
    double sinLatitudeCenter = Math.Sin(view.center.Latitude * Math.PI / 180);
    double pixelXCenter = ((view.center.Longitude + 180) / 360) * 256 * Math.Pow(2, view.zoom);
    double pixelYCenter = (0.5 - Math.Log((1 + sinLatitudeCenter) / (1 - sinLatitudeCenter)) / (4 * Math.PI))
 * 256 * Math.Pow(2, view.zoom);
    
    //calculate pixel coordinate of location
    double sinLatitude = Math.Sin(latlong.Latitude * Math.PI / 180);
    double pixelX = ((latlong.Longitude + 180) / 360) * 256 * Math.Pow(2, view.zoom);
    double pixelY = (0.5 - Math.Log((1 + sinLatitude) / (1 - sinLatitude)) / (4 * Math.PI)) 
* 256 * Math.Pow(2, view.zoom);

    //calculate top left corner pixel coordiates of map image
    double topLeftPixelX = pixelXCenter - (mapWidth / 2);
    double topLeftPixelY = pixelYCenter - (mapHeight / 2);

    //calculate relative pixel cooridnates of location
    double x = pixelX - topLeftPixelX;
    double y = pixelY - topLeftPixelY;

    return new Point((int)Math.Floor(x), (int)Math.Floor(y));
}

Listing 5 LatLong to Pixel coordinate method

Putting it all together

In order for these methods to work in a WinForm application the following references will have to be made:

using System.Drawing.Design;
using System.Drawing.Drawing2D;
using System.IO;

Listing 6 Required References

A reference to the Virtual Earth Imagery Service and the Virtual Earth token service will also be required.

The following global variables will also be needed:

string clientToken;
int mapWidth = 600;
int mapHeight = 400;

Listing 7 Global Variables

The following method will calculate the best map view for an array of locations, retrieve a map image from the imagery service, add custom icons to the map image and then display the map image in a picture box.

/// <summary>
/// Generate map image with custom icons
/// </summary>
/// <param name="points">Array of Location objects</param>
private void GenerateImage(Location[] points)
{
    //Get Best Map View
    BestView mapView = CalculateMapView(points);

    //Generate default map for best view
    Image map = GetMapImage(mapView);

    //add custom icons to map image
    AddCustomIcon(mapView, map, points);

    //display map
    MapImageBox.Image = map;
}

Listing 8 Method to Generate map image with custom icons

The following method can be used to add the custom icons to a map image:

/// <summary>
/// Adds custom icon image to map image
/// </summary>
/// <param name="view"></param>
/// <param name="map"></param>
/// <param name="points"></param>
private void AddCustomIcon(BestView view, Image map, Location[] points)
{
    //define the map as a gaphics object
    Graphics graphics = Graphics.FromImage(map);

    //Define custom icon
    Image customIcon = Image.FromFile("../../VEDefaultPin.gif");

    //custom icon offset with respect to the top left corner of icon
    Point imageOffset = new Point(12, 12);

    //Add custom icons to map image
    for (int i = 0; i < points.Length; i++)
    {
        Point point = LatLongToPixel(points[i], view);
        int x = point.X - imageOffset.X;
        int y = point.Y - imageOffset.Y;
        graphics.DrawImageUnscaled(customIcon, x, y);
    }
}

Listing 9 Method to add custom icons to a map image

This method can add a custom icon image and specify an image offset similar to the VECustomIconSpecification.

Here is an example of five locations plotted onto a map using Virtual Earth:

Figure 2 Points plotted using the JavaScript version of Virtual Earth

Here is an example of the same points being plotted onto a map image that is generated from the Virtual Earth Imagery Service using the methods in this article.

Figure 3 Custom icons on a VE Imagery Service image

Conclusion

Combine these methods with a valid client token and you will be able to add your custom icons to a map image. Complete sample code can be found here:
http://cid-e7dba9a4bfd458c5.skydrive.live.com/self.aspx/VE%20Sample%20code/CustomIconVEService.zip

This article was written by Richard Brundritt. Richard works for Infusion Development.

Similar to the method described in the post "Determine if Two Polygons Overlap" (http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!475.entry), we will have to iterate through all the line segments that make up one polygon and see if they overlap with any of the line segments that make up the second polygon. The difference will be that we will create an array of all the points of intersections. The following function takes in two arrays of coordinates that make up two polygons. An array of VELatLong objects will be returned.

//poly1 and poly2 are arrays of VELatlongs that represent polygons
function PolygonIntersections(poly1, poly2)
{
    var intersections = new Array();
    if(poly1.length >= 3 && poly2.length >= 3)
    {
        //close polygons
        poly1.push(poly1[0]);
        poly2.push(poly2[0]);
        
        for(var i = 0; i < poly1.length-1;i++)
        {
            for(var k = 0; k < poly2.length-1; k++)
            {
                var intersect = SimplePolylineIntersection(poly1[i],poly1[i+1],poly2[k],poly2[k+1]);
                
                if(intersect!=null)
                    intersections.push(intersect);
            }
        }
        
        return intersections;
    }
    
    return null;
}

The code for the function "SimplePolylineIntersection" can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!474.entry

Here is a screen shot of this function being used and the points being plotted on the map:

In order to determine if two polygons overlap we need to iterate through the line segments that make up each polygon and determine if any of the line segments intersect. Once we find a single point where the lines intersect we can then stop the iterations. The following code can be used to determine if two polygons overlap:

//poly1 and poly2 are arrays of VELatlongs that represent polygons
function ArePolygonsOverlapped(poly1, poly2)
{
    if(poly1.length >= 3 && poly2.length >= 3)
    {
        //close polygons
        poly1.push(poly1[0]);
        poly2.push(poly2[0]);
        
        for(var i = 0; i < poly1.length-1;i++)
        {
            for(var k = 0; k < poly2.length-1; k++)
            {
    if(SimplePolylineIntersection(poly1[i],poly1[i+1],poly2[k],poly2[k+1])!=null)
                    return true;
            }
        }
        
        return false;
    }
    
    return null;
}

The code for the function "SimplePolylineIntersection" can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!474.entry

The theory for calculating the point of intersection of two line segments is similar to the theory described in the post "Approximate Point of Intersection of Two Planes" found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!473.entry

The difference is that in order to be a point of intersection between two line segments, the calculated point must be within both of the bounding boxes in which the line segments are contained in. When can thus use the following two functions to determine the poin tof intersection of two polylines.

function SimplePolylineIntersection(latlong1,latlong2,latlong3,latlong4)
{
    //Line segment 1 (p1, p2)
    var A1 = latlong2.Latitude - latlong1.Latitude;
    var B1 = latlong1.Longitude - latlong2.Longitude;
    var C1 = A1*latlong1.Longitude + B1*latlong1.Latitude;
    
    //Line segment 2 (p3,  p4)
    var A2 = latlong4.Latitude - latlong3.Latitude;
    var B2 = latlong3.Longitude - latlong4.Longitude;
    var C2 = A2*latlong3.Longitude + B2*latlong3.Latitude;

    var determinate = A1*B2 - A2*B1;

    var intersection;
    if(determinate != 0)
    {
        var x = (B2*C1 - B1*C2)/determinate;
        var y = (A1*C2 - A2*C1)/determinate;
        
        var intersect = new VELatLong(y,x);
        
        if(inBoundedBox(latlong1, latlong2, intersect) && 
inBoundedBox(latlong3, latlong4, intersect))
            intersection = intersect;
        else
            intersection = null;
    }
    else //lines are parrallel
        intersection = null; 
        
    return intersection;
}

//latlong1 and latlong2 represent two coordinates that make up the bounded box
//latlong3 is a point that we are checking to see is inside the box
function inBoundedBox(latlong1, latlong2, latlong3)
{
    var betweenLats;
    var betweenLons;
    
    if(latlong1.Latitude < latlong2.Latitude)
        betweenLats = (latlong1.Latitude <= latlong3.Latitude && 
latlong2.Latitude >= latlong3.Latitude);
    else
        betweenLats = (latlong1.Latitude >= latlong3.Latitude && 
latlong2.Latitude <= latlong3.Latitude);
        
    if(latlong1.Longitude < latlong2.Longitude)
        betweenLons = (latlong1.Longitude <= latlong3.Longitude && 
latlong2.Longitude >= latlong3.Longitude);
    else
        betweenLons = (latlong1.Longitude >= latlong3.Longitude && 
latlong2.Longitude <= latlong3.Longitude);
    
    return (betweenLats && betweenLons);
}

A more complex method for calculating the point of intersection of two line segments in 3D space can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!439.entry

The point of intersection of two planes in 3D space can be calculated using spherical mathematics, as described here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!438.entry Taking this route requires a significant number of calculation. If we wanted to calculate the points of intersection of a large number of planes this can cause performance issues. The following is a method that can be used to calculate a 2D approximation for the point of intersection.

Theory

A 2D plane, commonly known as a line, can be defined using a formula of the following format:   Ax + By = C
Where x is in the same plane as the longitude coordinates and y is in the same plane as the latitude.
The components A, B, and C can be calculated using the following formula's:

A = y₂-y₁
B = x₁-x₂
C = A*x₁+B*y₁

If we were to have two planes defined by the following equations:

A₁x + B₁y = C₁
A₂x + B₂y = C₂

We can then calculate the determinate of these two equations using the following formula:

determinate = A1*B2 - A2*B1

If the determinate is 0, this would indicate that the planes are parallel to each other. If the determinate is 1 then this would indicate that the planes are perpendicular to each other.

Now that we have the determinate of these equations we can now solve for the x and y components of the point of intersection.

x = (B2*C1 - B1*C2)/determinate
y = (A1*C2 - A2*C1)/determinate

Application

The following code takes in four coordinates, where two coordinates represents the endpoints of a line segment. This code will return either a VELatLong representing the point of intersection or will return null which would indicate that the planes are parallel.

function SimplePlaneIntersection(latlong1,latlong2,latlong3,latlong4)
{
    //Line segment 1 (p1, p2)
    var A1 = latlong2.Latitude - latlong1.Latitude;
    var B1 = latlong1.Longitude - latlong2.Longitude;
    var C1 = A1*latlong1.Longitude + B1*latlong1.Latitude;
    
    //Line segment 2 (p3,  p4)
    var A2 = latlong4.Latitude - latlong3.Latitude;
    var B2 = latlong3.Longitude - latlong4.Longitude;
    var C2 = A2*latlong3.Longitude + B2*latlong3.Latitude;

    var determinate = A1*B2 - A2*B1;

    var intersection;
    if(determinate != 0)
    {
        var x = (B2*C1 - B1*C2)/determinate;
        var y = (A1*C2 - A2*C1)/determinate;
        intersection = new VELatLong(y,x);
    }
    else //lines are parrallel
        intersection = null; 
        
    return intersection;
}

By default Virtual Earth does not draw polylines across the international dateline (-180/180 degrees longitude) in 2D mode. The following code is an example of how to over come this issue. This code draws a polyline from Sydney Australia to Los Angeles California, which is a common airplane route. What this code does is determine the point of intersection of the path with the international date line. It then creates two polylines to represent this path. The algorithms use din this code are documented in the VE Math section of this blog.

<!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 src="http://dev.virtualearth.net/mapcontrol/mapcontrol.ashx?v=6.2"></script>
<script>
var earthRadius = 6367; //radius in km

var map;
var start = new VELatLong(-33.93979,151.17625);
var end = new VELatLong(33.94051,-118.40738);
var correctedPolyline1;
var correctedPolyline2;
var standardPolyline;

function GetMap()
{
    map = new VEMap('mymap');
    map.LoadMap(new VELatLong(0,0),1);
    
    var pin1 = new VEShape(VEShapeType.Pushpin,start);
    pin1.SetTitle("start");
    
    var pin2 = new VEShape(VEShapeType.Pushpin,end);
    pin2.SetTitle("end");

    standardPolyline = new VEShape(VEShapeType.Polyline,[start,end]);
    standardPolyline.HideIcon();
    
    var points = crossDateLine(start,end);
    
    correctedPolyline1 = new VEShape(VEShapeType.Polyline,[points[0],points[1]]);
    correctedPolyline2 = new VEShape(VEShapeType.Polyline,[points[2],points[3]]);
    
    correctedPolyline1.HideIcon();
    correctedPolyline2.HideIcon();
    
    map.AddShape([pin1,pin2,standardPolyline,correctedPolyline1,correctedPolyline2]);
    showLines();
    map.AttachEvent("onchangeview",mapModeChanged);
}

function mapModeChanged(e)
{
    if(e!=null)
    {
        showLines();
    }
}

function showLines()
{
    if(map.GetMapMode()==VEMapMode.Mode2D)
    {
        correctedPolyline1.Show();
        correctedPolyline2.Show();
        standardPolyline.Hide();
    }
    else
    {
        correctedPolyline1.Hide();
        correctedPolyline2.Hide();
        standardPolyline.Show();
    }
}

function DegtoRad(x)
{
    return x*Math.PI/180;
}

function RadtoDeg(x)
{
    return x*180/Math.PI;
}

function Cartesian(x,y,z)
{
    this.X = x;
    this.Y = y;
    this.Z = z;
}

function convertSphericalToCartesian(latlong)
{
    var lat = DegtoRad (latlong.Latitude);
    var lon = DegtoRad (latlong.Longitude);
    var x = earthRadius * Math.cos(lat)*Math.cos(lon);
    var y = earthRadius * Math.cos(lat)*Math.sin(lon);
    var z = earthRadius * Math.sin(lat);
            
    return new Cartesian(x,y,z);
}

function convertCartesianToSpherical(cartesian)
{
    var r = Math.sqrt(cartesian.X*cartesian.X + cartesian.Y*cartesian.Y+ cartesian.Z*cartesian.Z); 
    var lat = RadtoDeg(Math.asin(cartesian.Z/r));
    var lon = lon = RadtoDeg(Math.atan2(cartesian.Y, cartesian.X));
            
    return new VELatLong(lat,lon);
}

function crossProduct(p1,p2)
{
    var x = p1.Y*p2.Z - p1.Z*p2.Y;
    var y = p1.Z*p2.X - p1.X*p2.Z;
    var z = p1.X*p2.Y - p1.Y*p2.X;
    
    return new Cartesian(x,y,z);
}

function vectorLength(p1)
{
    return Math.sqrt(p1.X*p1.X+p1.Y*p1.Y+p1.Z*p1.Z);
}

function unitVector(p1)
{
    var length = vectorLength(p1);
    
    return new Cartesian(p1.X/length,p1.Y/length,p1.Z/length);
}

function intersectionOfPlanes(latlong1,latlong2,latlong3,latlong4)
{
    var point1 = convertSphericalToCartesian(latlong1);
    var point2 = convertSphericalToCartesian(latlong2);
    var point3 = convertSphericalToCartesian(latlong3);
    var point4 = convertSphericalToCartesian(latlong4);
    
    var vector1 = crossProduct(point1,point2);
    var vector2 = crossProduct(point3,point4);
    
    var unitVector1 = unitVector(vector1);
    var unitVector2 = unitVector(vector2);
    
    if((unitVector1.X != unitVector2.X )&& (unitVector1.Y != unitVector2.Y) && (unitVector1.Z != unitVector2.Z))
    {
        var vectorDirector = crossProduct(unitVector1,unitVector2);
        
        var intersection1 = unitVector(vectorDirector);
        var intersection2 = inverseCartesian(intersection1);    
        
        return new Array(convertCartesianToSpherical(intersection1),convertCartesianToSpherical(intersection2));
    }
    
    return -1;
}

function inverseCartesian(p1)
{
    return new Cartesian(-1*p1.X,-1*p1.Y,-1*p1.Z);
}

function crossDateLine(startLatLong,endLatLong)
{
    var points = intersectionOfPlanes(startLatLong,endLatLong,new VELatLong(90,-180),new VELatLong(-90,-180));
    
    var intersection1;
    
    if(points[0].Longitude==-180)
        intersection1 = points[0];
    else
        intersection1 = points[1];
        
    return [startLatLong,new VELatLong(intersection1.Latitude,180),intersection1,endLatLong];
    
}
</script>
</head>
   <body onload="GetMap();">
   <div id='mymap' style="position:relative; height:600px; width:800px;"></div>
   </body>
</html>

Sometimes it might be important to calculate the midpoint along a polyline. this can be done by first calculating the total length of a polyline and then to iterate through the points of the polyline, calculating the distance between points until  we find the points where the midpoint falls between. From here we can calculate the coordinate of the midpoint. The following functions can calculate the midpoint coordinate of a polyline.

function PolylineCentroid(points)
{
    var totalDistance = PolylineLength(points);
    
    var midpoint = FindMidPoint(points,totalDistance);
    return midpoint;
}

function PolylineLength(points)
{
    var distance = 0;
    
    for(var i=0;i<points.length-1;i++)
    {
        distance += haversineDistance(points[i],points[i+1]);
    }
    
    return distance;
}

function FindMidPoint(points,totalDistance)
{
    var midDistance = totalDistance/2;
    
    var distance = 0;
    var subDistance=0;
    var i;
    
    for(i=0;i<points.length-1;i++)
    {
        subDistance = haversineDistance(points[i],points[i+1]);
        
        if((subDistance+distance)<midDistance)
            distance += subDistance;
        else
            break;
    }
    
    subDistance = midDistance - distance;
    var bearing  = calculateBearing(points[i],points[i+1]);
    return calculateCoord(points[i], bearing, subDistance);
}

Information on the haversineDistance method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!317.entry

Information on the calculateBearing method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!393.entry

Information on the calculateCoord method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!400.entry

Theory

The gravitational acceleration that an object experiences decreases the further it gets from the center of the Earth. The gravitational acceleration at sea level at any point on the Earth’s surface can be calculated using the International Gravity formula:

When at sea level the gravitational acceleration only changes with respect to the latitude and not the longitude. The Earth is not a smooth surface and as such we need to be able to calculate this acceleration at altitudes with respect to sea level. We can modify the previous formula so that an altitude (h) in meters can be used to calculate the correct gravitational acceleration.

Application

This function will take in a latitude and an altitude in meters. It will return a number representing the gravitational accelerations whose units are in

function gravitationalAcceleration(latitude, altitude)
{
    var lat = DegToRad(latitude);
    
    return 9.780327*(1+0.0053024*Math.pow(Math.sin(lat),2)-0.0000058*Math.pow(Math.sin(lat),2))
-0.000003086*altitude;
}

Listing 1 Gravitational Acceleration function
The following post has additional information on the DegToRad method: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!257.entry

Theory

Triangulation is a common method for finding a coordinate based on the direction the point in question is in relative to two other points and the distance between these two other points. Triangulation is commonly used in surveying, navigation, wireless communications, astronomy, and just about anywhere it is necessary to calculate the location of an object from a distance.

Take for example two points A and B, each having a heading pointing towards some distant point C. Using triangulation we can calculate the coordinate of C. As you can see in the following figure this creates a triangle.

Figure 1 Triangulation

Using the law of sines we can make the following connection between the three points.

Using the formulas that have already been derived in this article we can calculate the distance between point A and B, we can also calculate the bearing from A to B. Once we know this we can calculate the angles and. Once we have these two angles we can calculate the angle using the following formula:

From here we can calculate the distances AC and BC using the following formulas:

We only need to calculate one of the distances. Once we have that we can use the distance, and the point is based at either A of B, and the associated heading to calculate the coordinate of C.

Application

This function for calculating the coordinate of a point using triangulation will take in two VELatLong objects and two bearings that are between 0 and 180 degrees. This function will return a VELatLong object for the point in which the provided information triangulates to. When calculating the angle between two headings we are only interested in the inner angle because we are working with a triangle all the angles will be under 180˚.

function triangulate(latlongA,latlongB,bearingAC,bearingBC)
{
    var AB = haversineDistance(latlongA,latlongB);
    var bearingAB = calculateBearing(latlongA,latlongB);
    var angleBAC = angleBetweenHeadings(bearingAB,bearingAC,"inner");
    var angleABC = angleBetweenHeadings(180-bearingAB,bearingBC,"inner");
    var angleBCA = 180 - angleBAC - angleABC;
    
    var AC = AB*Math.sin(DegToRad(angleABC))/Math.sin(DegToRad(angleBCA));
    
    return calculateCoord(latlongA, bearingAC, AC);
}

Listing 1 Triangulation Function Information on the calculateBearing method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!393.entry
Information on the haversineDistance method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!317.entry
Information on the calculateCoord method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!400.entry
The following post has additional information on the RadToDeg and DegToRad methods: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!257.entry
Information on the angleBetweenHeading method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!444.entry

Theory

Often it will be useful to know the angle between to polylines. By using the connecting point as our origin and calculating the bearings of the two polylines we can calculate the difference in the bearings to find out the angle between the polylines. There are two angles that we can calculate, the inner angle and the outer angle. The inner angle is less than 180˚.

Figure 1 Angle between polylines

We can calculate the angle using the following formula:

Application

A line segment consists to two end points. Since we have two line segments that are connected we have three points. Since there are two possible angles that can be calculated we will make the function so that the type of angle to be returned can be specified. The two angle types that can be used in this function is “inner” and “outer”. This function will return a number between 0 and 360 which will be the angle between the two polylines.

function angleBetweenPoints(latlongA,latlongB,latlongC,angleType)
{
    var headingBA = calculateBearing(latlongB,latlongA);
    var headingBC = calculateBearing(latlongB,latlongC);

    return angleBetweenHeadings(headingBA,headingBC,angleType);
}

Listing 1 Calculate angle between three points

We’ve broken the function to calculate the angle between three points up into two parts. By doing this we increase the usability of our code. This second function allows us to calculate the angle between two headings. This function will take in two headings between 0 and 360 degrees and will take in the type of angle that should be returned. A number representing the angle in degrees will be returned,

function angleBetweenHeadings(headingBA,headingBC,angleType)
{
    var angle = ((headingBA-headingBC)+360)%360;
        
    if(angleType=="inner")
        if(angle>180)
            return 360-angle;
        else
            return angle;
        
    if(angleType=="outer")
        if(angle<180)
            return 360-angle;
        else
            return angle;
}

Listing 2 Angle between two headings

Information on the calculateBearing method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!393.entry

Theory

The point of intersection is a location where two line segments overlap. We can use the function that calculates the intersection of two planes to find the two possible points of intersections. Once we have these two points of intersections we can then calculate the length of the two polylines and the length from each vertex of the polylines to a point of intersection. The sum of the lengths of the vertices should be the length of the polyline in which the vertices belong to if the point of intersection is on the polyline. If the first point of intersection is not on one of the polylines we can then do the same calculations for the second possible point of intersection.

Application

This function will take in the four VELatLongs that define the end points of the two polyline line segments. We will first calculate the two possible points of intersections. We then will calculate the lengths of the two polylines and the distances from the vertices to the first point of intersection. We will then take the length of each polyline and subtract it by the lengths of its vertices to the point of intersection. In theory this should be 0 if the point of intersection is on the line. Since there is a 1% accuracy we will then check to see if the absolute value of this calculation is within 1% of the length of the longest polyline. If the point is not on one of the polylines then the calculation is repeated for the second point of intersection. This function will either return a VELatLong object for -1 if the lines do not intersect.

function pointOfIntersection(latlong1,latlong2,latlong3,latlong4)
{
    var poi = intersectionOfPlanes(latlong1,latlong2,latlong3,latlong4); 

    if(poi!=-1)
    {
        var lengthP1P2 = haversineDistance(latlong1,latlong2);
        var lengthP1I1 = haversineDistance(latlong1,poi[0]);
        var lengthP2I1 = haversineDistance(latlong2,poi[0]);
        var lengthP3P4 = haversineDistance(latlong3,latlong4);
        var lengthP3I1 = haversineDistance(latlong3,poi[0]);
        var lengthP4I1 = haversineDistance(latlong4,poi[0]);
        var max = Math.max(lengthP1P2,lengthP3P4);
        if((Math.abs(lengthP1P2-lengthP1I1-lengthP2I1) < max*0.01) && 
(Math.abs(lengthP3P4-lengthP3I1-lengthP4I1) < max*0.01))
        {
            return poi[0];
        }
        else
        {
            var lengthP1I2 = haversineDistance(latlong1,poi[1]);
            var lengthP2I2 = haversineDistance(latlong2,poi[1]);
            var lengthP3I2 = haversineDistance(latlong3,poi[1]);    
            var lengthP4I2 = haversineDistance(latlong4,poi[1]);
            if((Math.abs(lengthP1P2-lengthP1I2-lengthP2I2) < max*0.01) && 
(Math.abs(lengthP3P4-lengthP3I2-lengthP4I2) < max*0.01))
            {
                return poi[1];
            }
        }
    }
    return -1;
}

Listing 1 Point of Intersection between two Polylines

Information on the haversineDistance method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!317.entry

Information on the intersectionOfPlanes method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!438.entry

Theory

We can think of each polyline as a line segment that represents a plane defined by the end points of the line segment and the center of the Earth. The point of intersection would be the point where two planes intersects which would create a line that passes through the center of the Earth, thus representing a point on the surface of the Earth.

Figure 1 Intersection of two planes defined by line segments

We will look at the solution when Cartesian coordinates are used. A line segment can be defined by two points P₁ and P₂. Each point having an x, y and z component. We can define a plane that contains P₁, P₂ and the center of the Earth (P₀) taking the dot product of P₀ and the cross product of P₁ and P₂ which will equal zero.

We can evaluate this as follows:

We can represent a second line segment the same way which consists of points P₃, and P₄.

We can then solve for x and Y in terms of Z as follows:

The point of intersection with this line and the sphere of radius r has z such that the distance from the center of the Earth is r. Thus;

From this we can solve for z.

Knowing z we can find x and y, however the equations are starting to get fairly large. We can significantly reduce the number of calculations that need to performed to solve x, y and z by first solving sub components. For instance, let us define the following subcomponents.

We can then write simplified versions of x, y, and z:

We could then convert this to work directly with spherical coordinates; however the equations grow in size fairly quickly. It is better to convert the coordinates to Cartesian coordinates and then perform our calculations and then convert the points back into spherical coordinates.

Since the Earth is Spherical there will be two points on the surface of the Earth where these planes will intersect. This second point is the inverse coordinate.

Application

To calculate the intersection of two planes we have to define the planes with line segments. A line segment will be defined by two VELatLong objects. Since we want to find to intersection of two planes we need to pass four VELatLong objects into our function to represent the two line segments. We then have to convert these VELatLong objects into Cartesian coordinates. From there we can calculate two vectors that represent the two planes. We can then calculate the unit vectors of each plane and check that the planes are not equal. If the planes are equal we will return -1. If the planes are not equal we can calculate the vector director. The unit vector of the vector director is the first point of intersection. We can take the first intersection point and calculate its inverse coordinate to calculate the second point of intersection. We then have to convert our x, y, and z parameters into a VELatLong objects. We can then return an array containing these two VELatLongs. Note that this algorithm will calculate the point of intersections to within a distance of 1% of the largest polyline length used to define a plane.

function intersectionOfPlanes(latlong1,latlong2,latlong3,latlong4)
{
    var point1 = convertSphericalToCartesian(latlong1);
    var point2 = convertSphericalToCartesian(latlong2);
    var point3 = convertSphericalToCartesian(latlong3);
    var point4 = convertSphericalToCartesian(latlong4);
    
    var vector1 = crossProduct(point1,point2);
    var vector2 = crossProduct(point3,point4);
    
    var unitVector1 = unitVector(vector1);
    var unitVector2 = unitVector(vector2);
    
    if((unitVector1.X != unitVector2.X )&& (unitVector1.Y != unitVector2.Y) && 
(unitVector1.Z != unitVector2.Z))
    {
        var vectorDirector = crossProduct(unitVector1,unitVector2);
        
        var intersection1 = unitVector(vectorDirector);
        var intersection2 = inverseCartesian(intersection1);    
        
        return new Array(convertCartesianToSpherical(intersection1),
convertCartesianToSpherical(intersection2));
    }
    
    return -1
}

Listing 1 Intersection Point of two Planes

Information on the convertSphericalToCartesian and convertCartesianToSpherical methods can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!280.entry

Information on the crossProduct method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!365.entry

Information on the unitVector method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!374.entry

Information on the inverseCartesian method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!342.entry

Theory

The midpoint of a line segment is the point on the line segment that is equidistance from each end point. By using functions which we have already created we can easily calculate the midpoint of a line segment. The first step is to find the length of the line segment. This can be found be sending the end points into the Haversine formula. Our mid point is half this distance away from either of our end points. We then have to calculate the bearing from one end point to the other. Which every end point we choose as the origin of our bearing is the point we will pass into the function that calculates a coordinate based on a point of origin, a bearing and a distance.

Application

To calculate the midpoint of a line segment we can use the Haversine function to find the distance between the end points. We can then calculate the bearing between the points. Once we have the distance to the midpoint and the bearing we can then calculate the midpoint coordinate. This function will take in two VELatLong objects and will return one VELatLong object.

function midpoint(latlong1, latlong2)
{
    var arcLength = haversineDistance(latlong1,latlong2);
    var brng = calculateBearing(latlong1,latlong2);
    
    return calculateCoord(latlong1, brng, arcLength/2);
}

Listing 1 Midpoint between two coordinates
Information on the haversineDistance method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!317.entry

Information on the calculateBearing method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!393.entry

Information on the calculateCoord method can be found here: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!400.entry

Theory

Sometimes it’s useful to be able to calculate the coordinates of a destination based off a point of origin, a bearing (β) and a distance. One popular activity that uses this is orienteering. Orienteering is a competition which is a timed race where participants use a map, and a compass to navigate between check points. At each check point the participant will receive the bearing and distance to the next check point based on the current location of the check point. The distance that is used is the arc length (L) between the two points. Knowing this we can derive the central angle (δ) between the two points based at the center of the Earth.

Using the central angle and trigonometry we can calculate the latitude of the destination point using the following formula:

Now that we know the latitude of the destination point we can calculate the longitude by using the following formula:

Application

The function to calculate the destination coordinate will take in three parameters; a point of origin, a bearing, and a distance. With these parameters we can calculate the latitude and longitude coordinates of the destination point. This function returns a VELatLong object.

var earthRadius = 6367; //radius in km
function calculateCoord(origin, brng, arcLength)
{
    var lat1 = DegtoRad(origin.Latitude);
    var lon1 = DegtoRad(origin.Longitude);
    var centralAngle = arcLength /earthRadius;

    var lat2 = Math.asin( Math.sin(lat1)*Math.cos(centralAngle) + Math.cos(lat1)*Math.sin(centralAngle)*
Math.cos(DegtoRad(brng)));
    var lon2 = lon1+ Math.atan2(Math.sin(DegToRad(brng))*Math.sin(centralAngle)*Math.cos(lat1),
Math.cos(centralAngle)-Math.sin(lat1)*Math.sin(lat2));
                
    return new VELatLong(RadtoDeg(lat2),RadtoDeg(lon2));
}

Listing 1 Function to calculate a destination coordinate

The following post has additional information on the RadToDeg and DegToRad methods: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!257.entry

Theory

Bearing is the direction in which an object is pointing or traveling. This is also commonly known as a heading. In Virtual Earth polygons and polylines consist of an array of coordinates that represent endpoints of line segments. It’s useful to think of these line segments as vectors. Each vector has a length and a bearing. Taking this into consideration will allow us to come up with some very useful algorithms.

Virtual Earth uses bearings in the 3D map mode in the form of the heading which is the direction in which the camera is pointing. This heading uses 0 degrees as true North, 90 degrees as East, 180 degrees as South, and 270 degrees as West. These same bearings can be used in 2D map mode when representing the line segments of a polygon or polyline as vectors. This calculation requires us to use complex variable calculations to calculate the arctangent of a point in space. This calculation is often done using the atan2 function. The reason for this is that it can handle situations when the denominator is equal to 0. This complex function has the following complex argument:

For

For y = 0

Where φ is an angle in the range such that:

The sgn function, known as the sign function and signum function, is another complex function that extracts the sign of a real number. This function has the following characteristics:

The following formula can be used to calculate the bearing between two coordinates. Note that this formula takes in latitude and longitudes as radian values and returns the bearing as a radian value in the following range:

We can convert this bearing to degrees in the range [0,360] by using the following formula:

The reason for adding 360 to the previous formula, then calculate the remainder when divided by 360, is to ensure that the bearing is a positive angle between 0 and 360.

Application

The function we are going to create will take in two a VELatLong objects and return the bearing in degrees. We will have to take the latitudes and longitudes of each VELatLong and convert them to their radian value. This will be done using the DegtoRad function we created earlier. We will then calculate the two points we need to put into the arctangent function. We will then calculate the bearing in degrees.

function calculateBearing(latlong1,latlong2)
{
    var lat1 = DegtoRad(latlong1.Latitude);
    var lon1 = latlong1.Longitude;
    var lat2 = DegtoRad(latlong2.Latitude);
    var lon2 = latlong2.Longitude;
    var dLon = DegtoRad(lon2-lon1);
    var y = Math.sin(dLon) * Math.cos(lat2);
    var x = Math.cos(lat1)*Math.sin(lat2) - Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon); 

    var brng = (RadtoDeg(Math.atan2(y, x))+360)%360;
    return brng;
}

Listing 1 Bearing calculation function

The following post has additional information on the RadToDeg and DegToRad methods: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!257.entry

Theory

A unit vector is vector that has a length of one. Unit vectors are used to indicate directions. Any vector can be expressed as a product of a magnitude and a direction. The magnitude of a vector is its length. The following formula can be used to calculate the magnitude of a vector :

The unit vector of can be calculated using the following formula:

Application

The following algorithm calculates the length of vector. This algorithm takes in a Cartesian coordinate and returns a float representing the magnitude of the Cartesian coordinate.

function vectorLength(p1)
{
     return Math.sqrt(p1.X* p1.X+ p1.Y* p1.Y+ p1.Z* p1.Z);
}

Listing 1 Vector Length

The following algorithm calculates the unit of vector a Cartesian coordinate. This algorithm takes in a Cartesian coordinate and returns a Cartesian coordinate representing the unit vector.

function unitVector(p1)
{
     var length = vectorLength(p1);

     return new Cartesian(p1.X/length,p1.Y/length,p1.Z/length);
}

Listing 2 Unit Vector
The following post has additional information on the Cartesian object: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!280.entry

Theory

When working with spherical coordinates it is often useful to be able to perform cross and dot product calculations. The cross product of two Cartesian coordinates, and is defined as a Cartesian coordinates that is perpendicular to the plane of and .

Where

The dot product of two Cartesian coordinates, and is defined as a length of the projection of onto the vector .

Application

The following algorithm calculates the cross product of two Cartesian coordinates. This algorithm takes in two Cartesian coordinates and returns a Cartesian coordinate representing the cross product vector.

function crossProduct(p1,p2)
{
    var x = p1.Y*p2.Z - p1.Z*p2.Y;
    var y = p1.Z*p2.X - p1.X*p2.Z;
    var z = p1.X*p2.Y - p1.Y*p2.X;
    
    return new Cartesian(x,y,z);
}

Listing 1 Cross Product function

The following algorithm calculates the dot product of two Cartesian coordinates. This algorithm takes in two Cartesian coordinates and returns a float representing the dot product.

function dotProduct(p1,p2)
{
      return p1.X*p2.X+p1.Y*p2.Y+p1.Z*p2.Z;
}

Listing 2 Dot Production function

The following post has additional information on the Cartesian object: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!280.entry

Theory

The inverse of a spherical coordinate can be calculated by taking the inverse of the latitude (θ) and offsetting the longitude (Φ) by 180˚. The following formulas can be used to calculate these values.

The sgn function, known as the sign function and signum function, is a complex function that extracts the sign of a real number. This function has the following characteristics:

The inverse of a Cartesian coordinate can be calculated by taking the x, y, and z components of the coordinate and multiply them by -1.

Application

The function to calculate the inverse coordinate will take in a VELatLong object and will return a VELatLong object. There is not a sgn function in the JavaScript Math class but we can do this using simple logic.

function inverseSpherical(latlong)
{
     var lat = -1*latlong.Latitude;
     var lon = 180-Math.abs(latlong.Longitude);

     if(latlong.Longitude>0)
         lon *= -1;

     return new VELatLong(lat,lon);
}

Listing 1 Inverse Spherical Coordinate

The following algorithm calculates the inverse coordinate of a Cartesian coordinate. This algorithm takes in a Cartesian coordinate and returns a Cartesian coordinate.

function inverseCartesian(p1)
{
     return new Cartesian(-1*p1.X,-1*p1.Y,-1*p1.Z);
}

Listing 2 Inverse Cartesian Coordinate

The following post has additional information on the Cartesian object: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!280.entry

Theory

A large majority of these calculations are based on the Earth being a perfect sphere with evenly geocentric latitudes. Making this assumption allows us to use these formulas in any part of the world with minimal error while at the same time keeping the mathematics simple enough to understand. In this article the radius of the Earth is used for several calculations. The radius of the Earth can be represented using any form of measurement, kilometers (km), miles (mi), etc. For these formulas we can define the radius in any units as long as we also define any distance variables in the same units. For most cases using the average radius of the Earth, 6367 km, will result in accurate enough results.

In the real world we know that the Earth is an ellipsoidal in shape and that the circles of latitude are further apart near the equator and closer together near the poles. Ellipses have a minor axis and a major axis. For Earth, the minor axis passes through the poles and has an approximate polar radius of 6,356.7523142 km (b). The major axis is equivalent to the radius at the equator with an approximate radius of 6,378.137 km (a). When Earth is assumed to be a sphere the latitude is the angle between the equatorial plane and a line from the center of Earth (r). In the real world the angle of latitude is the angle between the radius of curvature in the minor axis and the equatorial plane (r_N).

Figure 1 Geodetic Latitude in Ellipsoidal Earth

The formula to calculate the radius of curvature in the minor axis is as follows:

The formula for calculating the radius of the Earth at a specific geodetic latitude with respect to the center of the Earth is as follows:

Application

All the formulas in this article will implement a global variable for the radius of the Earth as defined below. We will also make the radius of Earth at poles and equator global variables. Rather than using a and b as variable names we will us equatorRadius (a) and polarRadius (b).

var earthRadius = 6367; //radius in km
var equatorRadius = 6378.137; //equatorial radius in km
var polarRadius = 6356.7523142; //polar radius in km 

Listing 1 Global Variable: Radius of the Earth

For calculating the radius of the Earth at a specific geodetic latitude we will create two functions. The first function will take in a VELatLong and return the radius and the second function will take a latitude in degrees.

function getEarthRadius(latlong)
{
      return getEarthRadiusAtLatitude(latlong.Latitude);
}

function getEarthRadiusAtLatitude(latitude)
{
     var lat = DegtoRad(latitude);
     return equatorRadius*Math.sqrt(Math.pow(polarRadius,4)/Math.pow(equatorRadius,4)*Math.pow((Math.sin(lat)),2)+
Math.pow(Math.cos(lat),2))/Math.sqrt(1-(1-(polarRadius*polarRadius)/( equatorRadius*equatorRadius))*
Math.pow(Math.sin(lat),2));
}

Listing 2 Radius of Earth at a specific geodetic Latitude

The following post has additional information on the DegToRad method: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!257.entry

There are several different ways to calculate the distance between two coordinates, P₁(r, θ, Φ) and P₂(r, θ, Φ). The first method consists of converting the points to Cartesian coordinates, and then calculating the cord length (d) between the points (straight line distance through the Earth). From there we can calculate the central angle between the two points (δ). Once we have the central angle between the two points we can make a simple calculation to figure out the arc length (L) between the two points.

Figure 1 Arc Length calculation

To calculate the cord length between P₁ and P₂ we can use Pythagorean Theorem as follows:

Now that we have the cord length we can calculate the central angle between the two points. Using simple trigonometry we can calculate from r and as follows:

We want to be in radians. If we calculate using javascript the angle will be in radians. If we calculate using a calculator we will have to convert the angle from degrees to radians. We can then calculate the arc length using the following formula:

The second method for calculating the distance between two points is the Haversine method. The basic principal is similar to the last method of calculating the distance between two points except that it does not solve for x, y, and z, but uses the formulas of the components. The cord length between the two points can be calculated by using Pythagorean Theorem against the displacement in the x, y, and z axis.

By using the following trigonometric identities we can derive a simpler formula for the cord length between two points.

Half-angle formula:

Double-Angle formula:

Subtraction formula:

We can then use the atan2 function to derive the central angle between the two points:

We can then calculate the distance between the two points using the following formula:

Both of these methods produce fairly accurate results. If more precision is required the Vincenty formula can be used which calculates the distance between two points on an ellipsoid. The basic principals are similar to the other two formulas except that the radius varies between the two points.

Application

For the function that calculates the arc length between two points it will take in to VELatLong objects and return a number. This function has to convert the VELatLongs into Cartesian coordinates. When we use the convertSphericalToCartesian function the returned value is a Cartesian coordinate object with X, Y, and Z values.

var earthRadius = 6367; //radius in km

function arcLength(latlong1, latlong2)
{
   var point1 = convertSphericalToCartesian(latlong1);
   var point2 = convertSphericalToCartesian(latlong2);
   var cordLength = Math.sqrt(Math.pow(point1.X-point2.X,2)+Math.pow(point1.Y-point2.Y,2)+
Math.pow(point1.Z-point2.Z,2));
   var centralAngle = 2*Math.asin(cordLength/2/earthRadius);
   return earthRadius*centralAngle;
}

Listing 1 Arc Length function

The Haversine function will take in to VELatLong objects and return a number. In this function we will need to convert the latitude and longitude values from degrees to radians using the DegtoRad function.

function haversineDistance(latlong1,latlong2)
{
    var lat1 = DegtoRad(latlong1.Latitude);
    var lon1 = DegtoRad(latlong1.Longitude);
    var lat2 = DegtoRad(latlong2.Latitude);
    var lon2 = DegtoRad(latlong2.Longitude);
    var dLat = lat2-lat1;
    var dLon = lon2-lon1; 
    var cordLength = Math.pow(Math.sin(dLat/2),2)+Math.cos(lat1)*Math.cos(lat2)*Math.pow(Math.sin(dLon/2),2); 
    var centralAngle = 2 * Math.atan2(Math.sqrt(cordLength), Math.sqrt(1-cordLength)); 
    return earthRadius * centralAngle;
}

Listing 2 Haversine function

The following post has additional information on the RadToDeg and DegToRad methods: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!257.entry

The following post has additional information on the convertSphericalToCartesian method: http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!280.entry