Currently when storing data in the customer service site you have the ability to query your point location data by entity id, property, radius search, and find near route. There are two key search functionalities missing. The first is bounding box searching, and the second is point in polygon searches. A bounding box is simple 4 sided regular polygon and as such if we can perform a point in polygon search we can use the same algorithm to perform a bounding box search. This post is going to outline two methods that can be used to perform such a search against data that is stored in the Customer Service Site.

Method 1

Lets start off with a polygon (blue). If we find the maximum and minimum latitude and longitude coordinates that this polygon has we can create a bounding box (green) that encloses the polygon. We can then calculate the center point and radius from center point to a corner of this bounding box. Once we know this information we can then enclose this bounding box with a circle (red). The information used to create this circle can be used with the radius search tools that are currently available.

If you perform a radius search using the calculated information you will end up with a lot of extra data points. You can filter this data points by running them through a point in polygon algorithm like the one outlined here: http://msdn.microsoft.com/en-us/library/cc451895.aspx. You can then return the filtered data to the user. Using this method your polygon must fit inside a bounding box whose sides are no longer than 353.55 miles. The benefit of this approach is that there are only a few calculations required initially. One down side is that there is a lot of extra data that could be potentially returned.

Method 2

The second method is very similar to the first. The maximum and minimum latitude and longitude coordinates of the polygon are used to create a bounding box so that the center of the polygon can be calculated. However, instead of calculating the radius from the center point to a corner of the bounding box, we can loop through and calculate the radii’s from the center point to each point in the polygon. The largest radius can then be used to perform a radius search.

Using this approach should reduce the amount of data that is returned by radius search as it is able to enclose the polygon tighter than that previous method. The down side is the number of calculations that need to be performed initially. For complex polygons this can result in slower performance. Using this method  you are limited to a maximum radius of 250 miles.

Conclusion

Using either of these methods will allow you to perform both bounding box and point in polygon searches against your data in the customer service site. Which method should you use. If you are performing bounding box searches or are using complex, predefined polygons then the first method would be ideal. If youare letting the user create the polygon chances are there will only be a limited number of points and as such method 2 would be ideal.

I have put together an example that allows you to draw out polygons using the right mouse button (press the left mouse button when done). After you have drawn a polygon you can then use either method to query the FourthCoffeeShops data source. You can down load the sample code here: http://cid-e7dba9a4bfd458c5.skydrive.live.com/self.aspx/VE%20Sample%20code/CSSFindByPolygon.zip

Often when displaying polygons you only need or want to show the polygons that are in the viewable area of the map. By doing this an increase in performance should occur with the map control as there will be less data that will need to be handled by the control. One method to go about this is to store the bounding coordinates of your polygon compare then to the bounding box of the viewable map. This leaves you with the task of determining if two bounding boxes overlap which should be much easier than determining if a polygon with any number of sides is in the map view. What this blog is meant to do is to present a mathematical solution to this problem.

For this solution we need to break a bounding box into some key components. A bounding box has a center point and two important radii’s. A radius from the center point to and edge in the x direction, and a radius in from the center point to and edge in the y direction. The following diagram illustrates this:

 

By looking at this diagram we can derive two key logical statements about how to determine if two bounding boxes overlap. The first statement is that if the radius from the center of bounding box A to bounding box B is less than or equal to the sum of the radius of A in the x direction and the radius of B in the x direction ; then bounding boxes A and B are aligned such that they share common coordinates in the x plane. We can write this like so: . Using this same train of thought we can make a similar statement for the y plane; . If both of these statements are true then bound bounding boxes must overlap or share a common edge.

We can calculate the x coordinate of the center point by adding the top left x value to the radius on the bounding box in the x direction. We can then determine the distance between centers by taking the absolute value of the different between the x coordinates of the center points of bounding box A and B. This gives use the distance which can be represented as follows:

The radius can be calculated by taking the difference between the top left x coordinate and the bottom right x coordinate and dividing it by 2. By applying this to the above formula we get the following:

this can then be further reduced to the following:

Similarly we can can derive a formula like this for :

Now we can look at the other side of the equation where we add the radii's of bounding box A and B. Doing this we get the following for the x direction:

this can then be reduced to the following:

Doing the same for the y direction we get the following:

this can then be reduced to the following:

We can now take the formula and put it all together:

this can then be reduced to the following:

Lets call this condition 1.

Doing the same for the y direction we get the following:

Lets call this condition 2.

Now if both condition 1 and condition 2 are true then the bounding boxes A and B must overlap.

We can now create a simple function that takes in two VELatLongRectangle objects and returns boolean depending on if the two bounding boxes overlap.

function DoBoundingBoxesIntersect(bb1, bb2) {                 //First bounding box, top left corner, bottom right corner                 var ATLx = bb1.TopLeftLatLong.Longitude;                 var ATLy = bb1.TopLeftLatLong.Latitude;                 var ABRx = bb1.BottomRightLatLong.Longitude;                 var ABRy = bb1.BottomRightLatLong.Latitude;                 //Second bounding box, top left corner, bottom right corner                 var BTLx = bb2.TopLeftLatLong.Longitude;                 var BTLy = bb2.TopLeftLatLong.Latitude;                 var BBRx = bb2.BottomRightLatLong.Longitude;                 var BBRy = bb2.BottomRightLatLong.Latitude;                 var rabx = Math.abs(ATLx + ABRx - BTLx - BBRx);                 var raby = Math.abs(ATLy + ABRy - BTLy - BBRy);                 //rAx + rBx                 var raxPrbx = ABRx - ATLx + BBRx - BTLx;                 //rAy + rBy                 var rayPrby = ATLy - ABRy + BTLy - BBRy;                 if(rabx <= raxPrbx && raby <= rayPrby)                 {                                 return true;                 }                 return false; }

 

In the past it has been useful to be able to determine the best map view to display an array of locations. In the AJAX control there is a method call SetMapView which determines the best map view for an array of coordinates. The VEWS and Silverlight control do not have this functionality. A long time ago someone wanted to know how to determine the best map view before loading the map. This required us to calculate out the best map view ourselves using map scale information. I later implemented the same method when working with the VEWS so that I could geo-reference images from the Imagery service (http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!488.entry). The method worked fairly well but I knew there was a better way to do this. When working with the Silverlight control I came across a situation where I needed this functionality again so I decided to create the updated method for determining the best map view using similar math that is used for the tiling system. This eliminated the need to maintain a list of scales for each zoom level and also reduce the amount of calculations that needed to be performed. Below is the updated algorithm for determining the best map view for a list of locations:

/// <summary> /// Calculates the best map view for a list of locations for a map /// </summary> /// <param name="locations">List of location objects</param> /// <param name="mapWidth">Map width in pixels</param> /// <param name="mapHeight">Map height in pixels</param> /// <param name="buffer">Width in pixels to use to create a buffer around the map. This is to keep pushpins from being cut off on the edge</param> /// <returns>Returns a MapViewSpecification with the best map center point and zoom level for the given set of locations</returns> public static MapViewSpecification BestMapView(IList<Location> locations, double mapWidth, double mapHeight, int buffer) {     MapViewSpecification mapView;     Location center = new Location();     double zoomLevel = 0;     double maxLat = -85;     double minLat = 85;     double maxLon = -180;     double minLon = 180;     //calculate bounding rectangle     for (int i = 0; i < locations.Count; i++)     {         if (locations[i].Latitude > maxLat)         {             maxLat = locations[i].Latitude;         }         if (locations[i].Latitude < minLat)         {             minLat = locations[i].Latitude;         }         if (locations[i].Longitude > maxLon)         {             maxLon = locations[i].Longitude;         }         if (locations[i].Longitude < minLon)         {             minLon = locations[i].Longitude;         }     }     center.Latitude = (maxLat + minLat) / 2;     center.Longitude = (maxLon + minLon) / 2;     double zoom1=0, zoom2=0;     //Determine the best zoom level based on the map scale and bounding coordinate information     if (maxLon != minLon && maxLat != minLat)     {         //best zoom level based on map width         zoom1 = Math.Log(360.0 / 256.0 * (mapWidth - 2*buffer) / (maxLon - minLon)) / Math.Log(2);         //best zoom level based on map height         zoom2 = Math.Log(180.0 / 256.0 * (mapHeight - 2*buffer) / (maxLat - minLat)) / Math.Log(2);     }     //use the most zoomed out of the two zoom levels     zoomLevel = (zoom1 < zoom2) ? zoom1 : zoom2;     mapView = new MapViewSpecification(center, zoomLevel);     return mapView; }

Recently I was asked how to draw route lines in Birds eye so that they become dashed when behind a building. Needless to say this turned into a pretty good brain teaser as the dashed lines on the Birds eye tiles are part of the image. So after some thought the possible methods that I could think of included:

• Use the 3D API of Virtual earth to determine where buildings are and see if based on the Birds eye orientation if a building is between the camera and the road. This could be done ahead of time and the data stored. It could also be done programmatically if the user is using 3D. This would require a significant amount of development. This method would also only work where 3D models exist.
• Create regions where tall buildings exist and then if a polyline is drawn in this region the direction in which the line runs can be determined. If the road runs horizontal more horizontal than vertical to the current birdseye orientation then make the line dashed. This method would be easier to develop and could be used dynamically, but would not be perfect. However, looking at some Birds eye imagery I’ve noticed that the dashed lines do not always show up when they should.

I ended up putting together some sample code on how to do the second method. This is were the code for the article "Dashed Polylines in Virtual Earth" originated from.

What the code does is defines a set of circular regions where tall buildings are known to exist (user defined) and then calculates the heading of a section of a path (http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!393.entry). If heading is more perpendicular than parallel to the Birds eye orientation the line is drawn as a dashed line.

As an added bonus, polylines done appear to be rendered when added in Birds eye mode. So the map needs to be flipped to a different map style, the polyline added, and then flipped back to Birds eye. For simplicity I used an array of points for my route. An actual VE route can be used when using client tokens and retrieving the route geometry (http://rbrundritt.spaces.live.com/blog/cns!E7DBA9A4BFD458C5!782.entry).

The complete source code for this sample can be found here: http://cid-e7dba9a4bfd458c5.skydrive.live.com/self.aspx/VE%20Sample%20code/dashedRoutes.zip

Here is a screen shot of this code in action:

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