Thursday, 15 October 2015

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


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.?)



Answer



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.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...