Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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 (a2 + b2 = c2) 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.?)

share|improve this question
2  
Please create a new question instead of modifying your original. Otherwise, you subvert all the mechanisms on this site: what do ratings mean when two or more questions are in play in one thread? What would it mean to mark one answer as correct? Etc. –  whuber Aug 23 '10 at 15:15
1  
@whuber: Although all that's been written in this thread is still on-topic WRT the original question asked, I agree that there are now really two questions asked here. It's too late to change that now, but will keep your advice in mind for a next time. –  stakx Aug 23 '10 at 21:06
add comment

3 Answers 3

up vote 9 down vote accepted

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.

share|improve this answer
    
The reason for the discrepancy is likely because the line segment (in the PCS) connecting the points is not a projected geodesic, which would be curvilinear when projected. Accordingly you get a smaller value than you should. A test of this theory is simple to do: take a simple geodesic (such as the equator) and compare two calculations of the distance between two widely separated points on the geodesic. One is the direct calculation, as in your code; the other breaks the geodesic into segments, directly calculates the segment lengths, and adds them up. The latter should be more accurate. –  whuber Aug 23 '10 at 15:13
add comment

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!

share|improve this answer
1  
Melita - Great to see you here! –  Kirk Kuykendall Oct 7 '10 at 3:00
1  
I agree, welcome aboard! –  matt wilkie Oct 13 '10 at 21:48
add comment

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.

share|improve this answer
    
@whuber: Thanks for answering. Concerning the first paragraph: I figured ArcGIS might see that the CH1903 map projection (which uses the Bessel 1841 ellipsoid) is used, and then project the points back onto that ellipsoid via the datum, and then do distance calculations on the ellipsoid. From your answer I take that ArcGIS won't do all that and will remain in Euclidean X-Y space to do calculations. (How about other GIS software?) -- 2nd paragraph: You're right of course, thanks for clarifying this point. –  stakx Aug 20 '10 at 18:32
    
A hidden reprojection mechanism is possible only if the point objects maintain references to a projection. I don't believe they do. –  whuber Aug 20 '10 at 18:37
    
@whuber: Would it be sufficient (for more accurate calculations) to know the ellipsoid used for the projection? AFAIK, ArcGIS stores a reference to the used projection with each feature class (data layer). –  stakx Aug 21 '10 at 16:00
1  
Actually IPoint, which derives from IGeometry, has SpatialReference as a property. help.arcgis.com/en/sdk/10.0/arcobjects_net/componenthelp/… However, I don't think ReturnDistance uses it. It might be worth testing to see if that's changed though. –  Kirk Kuykendall Aug 21 '10 at 20:49
1  
@stakx I've updated my answer to include code that shows that setting spatialref has no impact on ReturnDistance. –  Kirk Kuykendall Aug 23 '10 at 3:36
show 5 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.