[GIS] Measuring Angle of two Points from North in ArcGIS 10.3 C# SDK

anglesarcgis-10.3arcobjectsazimuthc

What is the most reliable method to calculate angle of a line from north, defined by it’s from and to points? I used to use IMeasurementTool interface found in defense solutions namespace found in Military Analyst for this purpose when we were using ArcGIS 10.0.

Since Military Analyst is no longer available for version later that 10.1, I need to find a method to replace its functionality. I have tested several methods for this purpose. Example code that uses same methods in a similar fashion can be found bellow ;

By constructing a line and using its angle;

IPoint from_point;
IPoint to_point;

To and from points are available with two point coordinate chosen by the user.

ILine2 measurement_line = new LineClass();
measurement_line.SpatialReference = from_point.SpatialReference;
measurement_line.FromPoint = from_point;
measurement_line.ToPoint   = to_point;
Double angle_from_horizontal = (180.0 * measurement_line.Angle) / Math.PI;
Double angle_from_north      = 90.0 - angle_from_horizontal;

By using IConstructAngle interface and GeometryEnvironmentClass object.

IConstructAngle angle_measurement = new GeometryEnvironmentClass();
IPoint north_point = new PointClass();
north_point.SpatialReference = from_point.SpatialReference;
north_point.X = from_point.X;
north_point.Y = 90.0;
FLOAT64 angle_radian = aci_olcer.ConstructThreePoint(north_point, from_point, to_point);
Double angle_from_horizontal = (180.0 * angle_radian) / Math.PI;
Double angle_from_north      = angle_from_horizontal;

By constructing a geodetic line and QueryTangent method of ICurve interface in order to find the tangent line at the start point.

PolylineClass geodeticline = new PolylineClass();
geodeticline.SpatialReference = from_point.SpatialReference;
IConstructGeodetic geodetic_calculations = (IConstructGeodetic)geodeticline;  

geodetic_calculations.ConstructGeodeticLineFromPoints(esriGeodeticType.esriGeodeticTypeGeodesic,  from_point,  to_point,  linearUnitMeter, esriCurveDensifyMethod.esriCurveDensifyByDeviation,  0.01);

ICurve tangent_finder = (ICurve)geodeticline;
ILine tanjant_line = new LineClass();
tangent_finder.QueryTangent(esriSegmentExtension.esriNoExtension, 0,  true,  0.0000001, tanjant_line);
Double angle_from_horizontal  = (180.0 * tanjant_line.Angle) / Math.PI;
Double angle_from_north       = 90 - angle_from_horizontal;

By using new MensurationClass.

IDistanceMeasurement distanceResult = new DistanceMeasurementClass();
IAngularMeasurement angularResult = new AngularMeasurementClass();
double result;
double uncertainty;
double result_angluar = 0;
double uncertainty_angular;

if (measure.CanMeasure == true)
{
measure.GetDistance(test_1_pnt, test_2_pnt, out distanceResult);
result = distanceResult.DistanceMeasurement;
uncertainty = distanceResult.DistanceUncertainty;
measure.GetAzimuthAngle(test_1_pnt, test_2_pnt, out angularResult);
result_angluar = angularResult.AngleMeasurement;
uncertainty_angular = angularResult.AngleUncertainty;
}

I compared the results of the methods I enlisted with the result of IMeasurementTool method and none of the produced the results produced by the MeasurementTool. Moreover this angle is actually used for later calculations where several points constructing corners of a rectangle are found by using distance, angle and starting point as arguments. If the angle measurements are not reliable, results deviate from original points. For example, I found out the by using other methods, resultant rectangle is deviated at an angle of 4 to 5 degrees from the original line for lines having a length of 200 000 meters located at whereabouts of Scandinavia. I also tested this for a wide range of locations throughout the World Map and found that resultant deviation is depended on the latitude, as results on the equator has near zero deviation where deviations near the poles are huge, and also depended on the length of the original line.

This is somewhat understandable because I found out that, as far as I understand, all of the methods I enlisted apart from the geodesic line construction returns angles calculated by using latitude and longitude values as Cartesian coordinates. Geodesic line construction on the other hand uses QueryTangent method to calculate angle at from_point. Results of the geodetic line are closer to the results on of the MeasurementTool but not same again. I assume that this may be due to the QueryTangent method used to calculate the angle at from point if it uses Cartesian latitude and longitude rations for arctangent. Also since I couldn’t find any sources on the implementation details, these are all my guesses.

Since when we approach the poles, actual distance between observed distances on the map changes, it is understandable to have an increased amount of error near the poles. So I used another method where I calculate geodetic distances of the two sides of a triangle formed by to and from points and a point representing the corner of a right triangle. Later on I used geodesic distance values to calculate an arctangent.

IGeometryServer2 geometry_server = new GeometryServerImplClass();
Double lat_y_dist_f64;
Double long_x_dist_f64;
IPoint corner_point = new PointClass();
corner_point.SpatialReference = spatialReference;
corner_point.X = to_point.X;
corner_point.Y = from_point.Y;
lat_y_dist_f64 = geometry_server.GetDistanceGeodesic(spatialReference, corner_point, to_point,  linearUnitMeter);
long_x_dist_f64 = geometry_server.GetDistanceGeodesic(spatialReference, from_point,  corner_point, linearUnitMeter);
Double angle_radian = Math.Atan2(long_x_dist_f64, lat_y_dist_f64);
Double angle_from_north = (180.0 * angle_radian) / Math.PI;

Results were much closer to their results of the MeasurementTool where originally 5 degrees of deviation observed at whereabouts of 40 degree latitude reduced to 0.003 degree for 2000 meter lines. While results were almost identical o the results of the MeasurementTool, they were not exactly same and deviation again increased near the poles and for longer lines.

So my question is what am I missing here? What is the most reliable and accurate method for calculating angles from north for geometries defined by two points? I believe that since this question is a very basic and fundamental question, ArcGIS should have a far easier solution to this problem which I am simply overlooking.

Best Answer

I've used Richie Carmichael's code from his blog post. In a later post he said he fixed some memory leaks and uploaded to arcscripts, however that link is broken. I'll look for the code - I think I downloaded it years ago. Anyway, here is the code from his first post. Take a look at the azimuth outputs.

The code uses the "projection engine" dll. I've read this dll is used throughout numerous Esri toolsets.

public class ProjectionEngineTest{
    public void Main() {
        double distance = 0d;
        double azimuthForward = 0d;
        double azimuthBack = 0d;
        ProjectionEngine projectionEngine = ProjectionEngineExplorer.GetInstance();
        projectionEngine.GetGeodesicDistance(
            ProjectionEngine.GLOBERADIUS,
            0,
            ProjectionEngine.ToRadians(51.50d),
            ProjectionEngine.ToRadians(-0.12d),
            ProjectionEngine.ToRadians(48.86d),
            ProjectionEngine.ToRadians(2.35d),
            out distance,
            out azimuthForward,
            out azimuthBack);
        MessageBox.Show(string.Format("London to Paris is {0} meters", distance.ToString()));               
    }
}

using System;
using System.Collections.Generic;
using System.IO;
using System.Runtime.InteropServices;
using System.Text;
using System.Windows.Forms;
using Microsoft.Win32;

namespace ESRI.ArcGIS.Sample {
    /// <summary>
    /// Projection Engine Abstract Class
    /// </summary>
    public abstract class ProjectionEngine {
        public const double GLOBERADIUS = 6367444.65712259d;
        private GeodesicDistance m_gd = null;
        private GeodesicCoordinate m_gc = null;

        [DllImport("kernel32.dll")]
        private static extern IntPtr LoadLibrary(String dllname);

        [DllImport("kernel32.dll")]
        private static extern IntPtr GetProcAddress(IntPtr hModule, String procname);

        private delegate void GeodesicDistance(
            [In] double semiMajorAxis,
            [In] double eccentricity,
            [In] double longitude1,
            [In] double latitude1,
            [In] double longitude2,
            [In] double latitude2,
            [Out] out double distance,
            [Out] out double azimuthForward,
            [Out] out double azimuthBack);

        private delegate void GeodesicCoordinate(
            [In] double semiMajorAxis,
            [In] double eccentricity,
            [In] double longitude1,
            [In] double latitude1,
            [In] double distance,
            [In] double azimuth,
            [Out] out double longitude2,
            [Out] out double latitude2);
        //
        // CONSTRUCTOR
        //
        protected ProjectionEngine() {
            IntPtr pe = ProjectionEngine.LoadLibrary(this.DLLPath);
            IntPtr gd = ProjectionEngine.GetProcAddress(pe, "pe_geodesic_distance");
            IntPtr gc = ProjectionEngine.GetProcAddress(pe, "pe_geodesic_coordinate");

            this.m_gd = (GeodesicDistance)Marshal.GetDelegateForFunctionPointer(
                gd, typeof(GeodesicDistance));
            this.m_gc = (GeodesicCoordinate)Marshal.GetDelegateForFunctionPointer(
                gc, typeof(GeodesicCoordinate));
        }
        //
        // PUBLIC METHODS
        //
        public abstract string DLLPath { get;}
        //
        // PUBLIC METHODS
        //
        /// <summary>
        /// Returns the geodestic azimuth and distance between two geographic locations.
        /// </summary>
        /// <param name="semiMajorAxis">Semi Major Axis</param>
        /// <param name="eccentricity">Globe Eccentricity</param>
        /// <param name="longitude1">From Longitude (radians)</param>
        /// <param name="latitude1">From Latitude (radians)</param>
        /// <param name="longitude2">To Longitude (radians)</param>
        /// <param name="latitude2">To Latitude (radians)</param>
        /// <param name="distance">Returned Geodetic Distance</param>
        /// <param name="azimuthForward">Returned Forward Azimuth (radians)</param>
        /// <param name="azimuthBack">Returned Reverse Azimuth (radians)</param>
        public void GetGeodesicDistance(
            double semiMajorAxis,
            double eccentricity,
            double longitude1,
            double latitude1,
            double longitude2,
            double latitude2,
            out double distance,
            out double azimuthForward,
            out double azimuthBack) {
            this.m_gd(
                semiMajorAxis,
                eccentricity,
                longitude1,
                latitude1,
                longitude2,
                latitude2,
                out distance,
                out azimuthForward,
                out azimuthBack);
        }
        /// <summary>
        /// Returns the geographic location based on the azimuth and distance from another geographic location
        /// </summary>
        /// <param name="semiMajorAxis">Semi Major Axis</param>
        /// <param name="eccentricity">Globe Eccentricity</param>
        /// <param name="longitude1">From Longitude (radians)</param>
        /// <param name="latitude1">From Latitude (radians)</param>
        /// <param name="distance">Distance from "From Location"</param>
        /// <param name="azimuth">Azimuth from "From Location"</param>
        /// <param name="longitude2">Out Logitude (in radians)</param>
        /// <param name="latitude2">Out Latitude (in radians)</param>
        public void GetGeodesicCoordinate(
            double semiMajorAxis,
            double eccentricity,
            double longitude1,
            double latitude1,
            double distance,
            double azimuth,
            out double longitude2,
            out double latitude2) {
            this.m_gc(
                semiMajorAxis,
                eccentricity,
                longitude1,
                latitude1,
                distance,
                azimuth,
                out longitude2,
                out latitude2);
        }
        /// <summary>
        /// Converts Radians to Decimal Degrees.
        /// </summary>
        /// <param name="degrees">In angle in decimal degrees.</param>
        /// <returns>Returns angle in radians.</returns>
        public static double ToRadians(double degrees) { return degrees * (Math.PI / 180d); }
        /// <summary>
        /// Converts Radians to Decimal Degrees.
        /// </summary>
        /// <param name="radians">In angle in radians.</param>
        /// <returns>Returns angle in decimal degrees.</returns>
        public static double ToDegrees(double radians) { return radians * (180d / Math.PI); }
    }
    /// <summary>
    /// Projection Engine Singleton for ArcGIS Explorer
    /// </summary>
    public sealed class ProjectionEngineExplorer : ProjectionEngine {
        private static ProjectionEngineExplorer projectionEngine;
        //
        // PUBLIC METHODS
        //
        public static ProjectionEngineExplorer GetInstance() {
            if (projectionEngine == null) {
                projectionEngine = new ProjectionEngineExplorer();
            }
            return projectionEngine;
        }
        public override string DLLPath {
            get {
                RegistryKey coreRuntime = Registry.LocalMachine.OpenSubKey(
                    @"SOFTWARE\ESRI\E2\CoreRuntime", false);
                object installDir = coreRuntime.GetValue("InstallDir", null);
                coreRuntime.Close();
                string folder = installDir.ToString();
                string bin = Path.Combine(folder, "bin");
                string pedll = Path.Combine(bin, "pe.dll");
                return pedll;
            }
        }
    }
    /// <summary>
    /// Projection Engine Singleton for ArcGIS Desktop
    /// </summary>
    public sealed class ProjectionEngineDesktop : ProjectionEngine {
        private static ProjectionEngineDesktop projectionEngine;
        //
        // PUBLIC METHODS
        //
        public static ProjectionEngineDesktop GetInstance() {
            if (projectionEngine == null) {
                projectionEngine = new ProjectionEngineDesktop();
            }
            return projectionEngine;
        }
        public override string DLLPath {
            get {
                RegistryKey coreRuntime = Registry.LocalMachine.OpenSubKey(
                    @"SOFTWARE\ESRI\CoreRuntime", false);
                object installDir = coreRuntime.GetValue("InstallDir", null);
                coreRuntime.Close();
                string folder = installDir.ToString();
                string bin = Path.Combine(folder, "bin");
                string pedll = Path.Combine(bin, "pe.dll");
                return pedll;
            }
        }
    }
}
Related Question