How does ArcGIS compute the distance between two points with a non-equidistant projection?

Question

This is a follow-up question to my previous one, Can you suggest some well-written introductory texts about coordinate system projections?

---

Let's assume I'm working with the CH1903 map projection, which for all I know is conformal, but not equidistant. Meaning, angles (shape) have been preserved, but not areas, distances, or scale. (At least these have not been preserved exactly). So far so good.

I'm wondering what kind of calculation ArcGIS performs when I now want to calculate the distance between two points. In ArcObjects, I could use the IProximityOperator interface as follows:

IPoint a = ...,
       b = ...;

double distance = ((IProximityOperator)a).ReturnDistance(b);

Question: When I'm working with a reference system that does not accurately preserve distances, what would ArcGIS do when I query it for the distance between two points (as shown above)?

• Does it simply do some Pythagorean maths (a² + b² = c²) to get the distance, meaning the returned distance will only be as accurate as the projection allows?

• Or will it do something more complicated, like some form of re-projection, to get a more accurate distance?

(The same question, but more generally: Once that geometries have been projected, does ArcGIS perform all calculations simply in Euclidean space, or does the used map projection still influence calculations of distances, angles, areas, etc.?)

Answer 1

If you want a stable method of computing geodesic distances, I recommend Richie Carmichael's wrapper for ESRI's Projection Engine.

Update: I just tried Richie's code with ArcGIS 10.0 on Vista64 and get an exception after calling LoadLibrary. I'll look into that more later.

For now though, here is some code in response to questions in the comments of another answer.

The code compares IProximityOperator for points with and without spatial references. Then it shows how to use an azimuthal equidistant projection (with first point being the point of tangency) to find the great circle distance.

private void Test()
{
    IPoint p1 = new PointClass();
    p1.PutCoords(-98.0, 28.0);
    
    IPoint p2 = new PointClass();
    p2.PutCoords(-78.0, 28.0);

    Debug.Print("Euclidian Distance {0}", EuclidianDistance(p1, p2));
    Debug.Print("Distance with no spatialref {0}", GetDistance(p1, p2));

    ISpatialReferenceFactory srf = new SpatialReferenceEnvironmentClass();
    IGeographicCoordinateSystem gcs =
    srf.CreateGeographicCoordinateSystem((int)esriSRGeoCSType.esriSRGeoCS_WGS1984);

    p1.SpatialReference = gcs;
    p2.SpatialReference = gcs;

    Debug.Print("Distance with spatialref {0}", GetDistance(p1, p2));
    Debug.Print("Great Circle Distance {0}", GreatCircleDist(p1, p2));

}
private double GetDistance(IPoint p1, IPoint p2)
{
    return ((IProximityOperator)p1).ReturnDistance(p2);
}

private double EuclidianDistance(IPoint p1, IPoint p2)
{
    return Math.Sqrt(Math.Pow((p2.X - p1.X),2.0) + Math.Pow((p2.Y - p1.Y), 2.0));
}

private double GreatCircleDist(IPoint p1, IPoint p2)
{
    ISpatialReferenceFactory srf = new SpatialReferenceEnvironmentClass();
    IProjectedCoordinateSystem pcs =
    srf.CreateProjectedCoordinateSystem((int)esriSRProjCSType.esriSRProjCS_WGS1984N_PoleAziEqui);
    pcs.set_CentralMeridian(true, p1.X);
    ((IProjectedCoordinateSystem2)pcs).LatitudeOfOrigin = p1.Y;
    p1.SpatialReference = pcs.GeographicCoordinateSystem;
    p1.Project(pcs);
    p2.SpatialReference = pcs.GeographicCoordinateSystem;
    p2.Project(pcs);
    return EuclidianDistance(p1, p2);
}

Here's the output:

Euclidian Distance 20
Distance with no spatialref 20
Distance with spatialref 20
Great Circle Distance 1965015.61318737

I think it would be interesting to test this against the projection engine dll (pe.dll). Will post results if I ever get Richie's code to work.

Update: Once I changed Richies code to compile for x86, I got it to run. Interesting ... the great circle distance it give me is 1960273.80162999 - a significant difference from that returned from the azimuthal equidistant method above.

Answer 2

In ArcGIS 10, check out IGeometryServer2 which now has GetDistanceGeodesic (geodesic distance between two geometries), GetLengthsGeodesic (return the geodesic lengths of each polyline), and DensifyGeodesic (densify a polyline by plotting points along the geodesic lines connecting the vertices, uses IPolycurve4::GeodesicDensify) methods.

As mentioned in the other answers, ArcGIS still uses mostly planar calculations.

Melita Kennedy

---

Some commentary on the other answers (not enough rep yet to comment directly!).

Esri's azimuthal equidistant projection supports ellipsoids. The GreatCircleDist code is creating a PCS that uses an ellipsoid/spheroid-based GCS, thus distances from the center/origin point will be geodesic distances, not great circle distances. It could also be simplified. We know the projected coords of the first point because it's the center of the projection: 0,0. So only the 2nd point needs to be projected. A simplified EuclidianDistance function could then be used.

I checked the results against the pe.dll's geodesic functions and it matched. It looks like Richie's app is using a sphere, so it is returning great circle distances/coordinates in its test application. That's why the results don't match. I didn't recognize the radius values; I think I need to talk to him about it!

Answer 3

The accuracy of any answer about ArcGIS is subject to change at any time--for all we know, new procedures will be introduced in the very next service pack without warning or documentation. That being said, ESRI software has for a long time used Euclidean calculations (e.g., the Pythagorean formula for distances) whenever projected coordinates are used. Often, in calculations like the one you illustrate, the software doesn't even have access to the projection information, so what else can be done?

Your question itself seems to suggest that Euclidean distance calculations for an equidistant projection are correct. Nothing could be further from the truth. For a one-point equidistant projection, the Euclidean distance to the base point is guaranteed to equal the geodesic distance; for a two-point equidistant projection, the Euclidean distance to either base point is guaranteed to equal the geodesic distance. In return for those guarantees, the metric distortion between all other pairs of points is typically greatly increased compared to other projections one might choose.