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

When representing the location of objects in three dimensions there are several different types of coordinate systems that can be used to represent the location with respect to some point of origin. Locations on a Virtual Earth map are represented by the spherical coordinate system. The spherical coordinate system uses three parameters to represent a point in space, a radial distance (r) from a point of origin to the point in space, a zenith angle (θ) from the positive z-axis, and azimuth angle (Φ) from the positive x-axis. In Virtual Earth the radial distance is the radius of the Earth, the zenith angle is the latitude, and the azimuth angle is the longitude. The origin is considered to be the center of the Earth. The Cartesian coordinate system is the standard x, y, and z coordinate system that is commonly used to represent a point in space. The following figure gives a graphical representation of the spherical coordinate system. It also points out the Cartesian coordinate system equivalent dimensions.

Figure 1 Spherical Coordinate System

By looking at the Spherical coordinate figure we can calculate the x, y, and z components of the Cartesian coordinate system. We will start with the z component because it is the easiest to calculate. Since the z component is in the same plane as r and the θ angle, and this plane is perpendicular to the xy plane which the Φ angle is in, we can see that the z component is equivalent to the opposite side of the triangle, and r is equivalent to the hypotenuse. Taking this information into consideration we can calculate the z component as follows:

To calculate the x and y components we need to first calculate the projection of the radius line in the xy plane. This projection would be equivalent to the hypotenuse of the two triangles in the xy plane. The following formula will give use the length of the projection of r:

The x component is equal to the projection of the r projection in the xz plane. We can calculate the x component using the following formula:

The y component is equal to the projection of the r projection in the yz plane. We can calculate the y component using the following formula:

Now that we can convert from Spherical to Cartesian coordinates, it would be useful to be able to do the reverse. Using Pythagorean Theorem we can calculate the length of r.

Using the formula for the z component we can derive the formula to calculate the angle θ.

By rearranging the formulas for the x and y components we can solve the value for the Φ angle.

Application

Before we convert the coordinates we need to create a definition for a Cartesian coordinate object. The following function is how we will define a Cartesian coordinate:

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

Listing 1 Cartesian Coordinate Object

For the conversion between Spherical and Cartesian coordinates we will take in a VELatLong object and use a constant value for the radius of the earth. If we want to perform more advance calculations we could extract the altitude value of the VELatLong object, if it’s relative to the WGS 84 ellipsoid, and add the radius of the earth to get the total distance from the center of the Earth to a point in space. We will return a Cartesian coordinate object that contains the x, y and z values. This will make it much easier to follow trace the code.

var earthRadius = 6367; //radius in km
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);
}

Listing 2 Spherical to Cartesian coordinate conversion

For the conversion from Cartesian coordinates to Spherical coordinates we will take in Cartesian coordinate object. This function will return a VELatLong that represents our spherical coordinates.

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 = RadtoDeg(Math.atan2(cartesian.Y, cartesian.X));
    return new VELatLong(lat,lon);
}

Listing 3 Cartesian to Spherical coordinate conversion

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

When working with trigonometric functions in JavaScript the units have to be in radians and not in degrees. Radians are a standard unit of angular measurement. A radian is the angle that is swept out when the radius of a circle and the swept out arc length are equal in length.

Figure 1 One Radian

A complete circle, 360º, is equivalent to 2π radians. Thus the conversion from degrees to radians is as follows:

The conversion from radians to degrees is as follows:

Application

The conversion between degrees and radians will be used often. Knowing this we will create two functions that will convert from degrees to radians and back again.

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

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

Listing 1 Conversion function for radians and degrees

Often it is useful to be able to calculate the center of a polygon. There are several ways to do this. One way is to calculate the average latitude and longitude points of the coordinates that make up the polygon. This method however can result in a point that is no where close to the center of the polygon if there is a large number of points grouped together. A second method is to calculate the center of mass of the polygon which gives us the center of a polygons centroid based on it's area. This can be done in 3-dimensions to account for the curvature of the earth, however assuming the earth is flat for small areas makes the math far easier. The standard formula for calculating the area of a polygon is as follows:

We can then use this area and calculate the x and y components of the centroid.

We can use these formulas to calculate an approximate centroid of a polygon.  The following function that takes in an array of coordinates that represent a polygon and returns a VELatLong object for the centroid.

function simplePolygonCentroid(points)
{
    var sumY = 0;
    var sumX = 0;
    var partialSum = 0;
    var sum = 0;
    
    //close polygon
    points.push(points[0]);
    
    var n = points.length;
    
    for(var i=0;i<n-1;i++)
    {
        partialSum = points[i].Longitude*points[i+1].Latitude - points[i+1].Longitude*points[i].Latitude;
        sum += partialSum;
        sumX += (points[i].Longitude+points[i+1].Longitude) * partialSum;
        sumY += (points[i].Latitude+points[i+1].Latitude) * partialSum;
    }
    
    var area = 0.5*sum;
    
    return new VELatLong(sumY/6/area,sumX/6/area);
}

Many times when performing calculations with polygons it is useful to remove any consecutively matching points and any null displacements that occur between three points in the array of points of a polygon. A null displacement is when you travel from point A to point B and then back to point A. To check if there are consecutively matching points we have to loop through the array of points that make up the polygon and compare consecutive points. If a match is found one of the points are removed from the data array. To check if there are any null displacements we will need to see if the first and third point of a series of three points is equal. If they are we have we will remove the second and third point from the data array.

Application

The algorithm to clean the data of a polygon will take in an array of VELatLongs and will return an array VELatLongs of the same size or less.

function cleanData(points)
{
    for(var i=0;i<points.length-1;i++)
    {
        if((points[i].Latitude==points[i+1].Latitude) && (points[i].Longitude==points[i+1].Longitude))
        {
            points.splice(i,1);
            i--;
        }
        
        if(i+2<points.length-1)
        {
            if((points[i].Latitude==points[i+2].Latitude) && (points[i].Longitude==points[i+2].Longitude))
            {
                points.splice(i,2);
                i-=2;
            }
        }
    }
    
    if((points[0].Latitude==points[points.length-1].Latitude)&&
(points[0].Longitude==points[points.length-1].Longitude))
        points.splice(points[points.length-1],1);
        
    return points;
}

Virtual Earth uses four different types of distance units; meters, feet, kilometers, and miles. Conversion between distance units consists of multiplying or dividing a number by a constant. There are several ways we can handle this situation, we can have twelve constants to represent every possible conversion between these four units, or we can make a function that uses three constants and calculates the conversion. We will implement the second method.

Theory

For the theory we will first go through the twelve different conversions and the connection between each unit.

Meters to Kilometers = Meters -> Kilometers =

Kilometers to Meters = Kilometers -> Meters =

Feet to Miles = Feet -> Miles =

Miles to Feet = Miles -> Feet =

Kilometers to Miles = Kilometers -> Miles =

Miles to Kilometers = Miles -> Kilometers =

By combining the above six conversions we can make all the rest of the possible conversions between the four distance units.

Meters to Miles = Meters -> Kilometers -> Miles =

Miles to Meters = Miles -> Kilometers -> Meters =

Feet to Kilometers = Feet -> Miles -> Kilometers =

Kilometers to Feet = Kilometers -> Miles -> Feet =

Feet to Meters = Feet -> Miles -> Kilometers -> Meters =

Meters to Feet = Meters -> Kilometers -> Miles -> Feet =

Application

To apply the conversions we will use if statements that will determine if the conversion that is being performed uses one of the six main conversion types. The function to convert between distance units will take in two parameters, a distance and the conversion type. The distance will be a number and the conversion type will be a string that represents the type of conversion that is to be made. The syntax for the conversion type parameter will consist of two units separated by the number two. For example, the conversion type for meters (m) to feet (ft) would be “m2ft”. The following four units are used for the conversion type parameter: meter (m), feet (ft), kilometers (km), and miles (mi).

function convertDistance(distance, conversionType)
{
    //meters to kilometers
    if(conversionType=="m2km"||conversionType=="m2mi"||conversionType=="m2ft") 
        distance /= 1000;
    
    //kilometers to meters
    if(conversionType=="ft2m"||conversionType=="mi2m"||conversionType=="km2m")
        distance *= 1000;        
    
    //feet to miles
    if(conversionType=="ft2mi"||conversionType=="ft2m"||conversionType=="ft2km")
        distance /= 5280;
    
    //miles to feet
    if(conversionType=="mi2ft"||conversionType=="m2ft"||conversionType=="km2ft")
        distance *= 5280;
        
    //kilometers to miles
    if(conversionType=="m2mi"||coversionType=="m2ft"||conversionType=="km2ft"||conversionType=="km2mi")
        distance /= 1.609344;
    
    //miles to kilometers
    if(conversionType=="mi2m"||conversionType=="ft2m"||conversionType=="ft2km"||conversionType=="mi2km")
        distance *= 1.609344;
        
    return distance;
}

This function takes in conversion types that consist of a from unit and to unit which is separated by the number two. The following are the unit short forms used:

• meters: m
• feet: ft
• miles: mi
• kilometers: km

By using the match regular expression method of JavaScript we can simplify the above algorithm by matching the starting and ending characters of the conversion type. This algorithm reduces to the following:

function convertDistance(distance, conversionType)
{
    //meters to kilometers
    if(conversionType.match("^m2")=="m2") 
        distance /= 1000;
    
    //kilometers to meters
    if(conversionType.match("2m$")=="2m")
        distance *= 1000;        
    
    //feet to miles
    if(conversionType.match("^ft2")=="ft2")
        distance /= 5280;
    
    //miles to feet
    if(conversionType.match("2ft$")=="2ft")
        distance *= 5280;
        
    //kilometers to miles
    if(conversionType=="m2mi"||coversionType=="m2ft"||conversionType=="km2ft"||conversionType=="km2mi")
        distance /= 1.609344;
    
    //miles to kilometers
    if(conversionType=="mi2m"||conversionType=="ft2m"||conversionType=="ft2km"||conversionType=="mi2km")
        distance *= 1.609344;
        
    return distance;
}

The following is a list of conversion types that can be passed:

• meters to kilometers: "m2km"
• meters to miles: "m2mi"
• meters to feet: "m2ft"
• feet to meters: "ft2m"
• feet to kilometers: "ft2km"
• feet to miles: "ft2mi"
• kilometers to miles: "km2mi"
• kilometers to feet: "km2ft"
• kilometers to meters: "km2m"
• miles to feet: "mi2ft"
• miles to meters: "mi2m"
• miles to kilometers: "mi2km"

Virtual Earth gives us the ability to get turn by turn directions which is great. However the directions may take longer to read through when you need to know which way to turn at the last minute. One way to aid this is to bold the key words in the directions so that they stand out. The following code is an simple example of how to display the plain text, un-bolded directions.

function OnGotRoute( route )
{
   var legs = route.RouteLegs;
   var leg = null;
   var directions = "";
   
   for (var i=0; i < legs.length; i++)
   {
      leg = legs[i];
      legItinerary = leg.Itinerary;
      
      for ( var j=0; j< legItinerary.Items.length; j++ )
      {
         item = legItinerary.Items[j];
         directions += item.Text+"<br/>";
      }
   }
   
   document.getElementById("directions").innerHTML = directions;
}

We can determine which words to bold by using regular expressions. Searching and bolding every keyword would require a list of all the possible key words which would result in a really long regular expression. However, there are a limited number of words that are commonly used by Virtual Earth in the directions which we don't want to bold. Knowing this we can bold the key words by bolding the complete directions instruction and un-bolding all the common words. The following code example is a modified version of the previous code that bolds the complete direction instruction and un-bolds all common words, thus bolding all keywords.

function OnGotRoute( route )
{
   var legs = route.RouteLegs;
   var leg = null;
   var directions = "";
   
   for (var i=0; i < legs.length; i++)
   {
      leg = legs[i];
      legItinerary = leg.Itinerary;
      
      for ( var j=0; j< legItinerary.Items.length; j++ )
      {
         item = legItinerary.Items[j];
         directions += "<b>"+item.Text.replace(/\b(Take| Turn|^Turn|Bear|Depart|onto|for|and|the|then|Keep
|follow|signs|toward|At|exit|ramp|stay|on|to|Make a|Road name changes|Arrive)\b/gi,"</b>$1<b>")+"</b><br/>";
      }
   }
   
   document.getElementById("directions").innerHTML = directions; 
}

The following is an example of set of bolded directions from zip code 12345 to zip code 04556:

Depart River Rd
Turn left onto Edison Ave
Turn right onto SR-5 / State St
Turn right onto Broadway
Take ramp left and follow signs for I-890 / Schenectady County Memorial Veteran Hwy East
Take ramp for I-90 East / Gov Thomas E Dewey Thruway East / New York State Thruway East toward Albany
Keep right onto I-87 South / Gov Thomas E Dewey Thruway South / New York State Thruway South
At exit 21A, take ramp right for Gov Thomas E Dewey Thruway East / New York State Thruway East toward Mass Tpke / Boston
Keep straight onto I-90 East / Gov Thomas E Dewey Thruway East / New York State Thruway East
At exit 11A, take ramp right for I-495 North toward Marlborough / N.H.-Maine
Keep straight onto I-95 North
At exit 44, take ramp right for I-295 North toward South Portland / Downtown Portland
At exit 28, take ramp right toward Coastal Route / Brunswick / Bath
Keep straight onto US-1 North / Blue Star Memorial Hwy
Keep left to stay on US-1 / Blue Star Memorial Hwy / Mill St
Keep left to stay on US-1
Make a U-turn at Cochran Rd
Turn left onto SR-27
Turn left onto Old County Rd
Road name changes to E Old County Rd
Turn left onto road
Arrive on the left

When version 6.0 of Virtual Earth was released it was announced that KML file types were supported. Unfortunately, it there were some delays and support was not added until after the release. It just happened that I needed to load some KML data into a Virtual Earth map so I threw together this program that translates KML files to GeoRSS. When KML support was added there were the following limitations:

• Maximum of 200 place markers
• Network links go only one level
• Icon, StyleUrl, IconStyle,LineStyle, and PolygonStyle are unsupported

By translating your KML files to GeoRSS format you gain the ability to add upwards of 10,000 points, and ability to follow network links to through all levels. This translator is setup to handle the following KML elements:

• Placemark
• Name
• Description
• Point
• Polygon (Only the outerBoundary)
• LineString
• MultiGeometry (Splits geometry up into individual shapes)
• Style
• StyleURL
• IconStyle
• NetworkLink (Note: this program will follow nested network links through all levels which could potentially result in an endless loop)

This program has the ability to either translate an individual KML or a go through a directory and translate all KML files within the directory.

You can download the complete source code here.

It has been a standard practice for MapPoint and Virtual Earth users to Geocode their location data using the Virtual Earth/MapPoint Customer Service site. Since the introduction of rooftop geocoding in Virtual Earth users have been asking for the ability to Geocode their location data with the same accuracy. This can now be done using the Virtual Earth geocoding service.

I have thrown together a simple program that can do this. You can download it here: http://cid-e7dba9a4bfd458c5.skydrive.live.com/self.aspx/VE%20Tools/BatchGeocodeWithVE.zip. This program takes in tab or pipe delimited data source files, similar to the Customer Service site. It fills in the latitude and longitude fields and adds three more fields; Confidence, MatchCode, and CalculationMethod.

Running a couple of addresses through this program you can see a difference in the accuracy of the coordinates provided by the VE Geocoding service compared to the coordinates provided by the Customer Service site.