The Generate Near Table tool (under Analysis Tools > Proximity) will get you started. You can give the tool a search radius around your reference points and tell it to find the 4 nearest other points, it will calculate the distance and angle to each of the 4 points from the reference point. Make sure the Angle box is checked, Find Only Closest Feature is unchecked, and you enter the appropriate number of nearest points to find:
Once that's done you will have angles to each point from the reference point, so you'll need to do a little math to figure out the angles in your screenshot. The Near tool calculates angles relative to the X axis so an angle of 0 would be due east, 90 would be north and -90 would be south.
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;
}
}
}
}
Best Answer
For accurate calculations, convert (lat, lon, elevation) directly to earth-centered (x,y,z). (If you don't do this, you need to retain additional information about the local normal ["up"] directions in order to compute angles accurately at nonzero elevations.)
Elevation
Given two points (x,y,z) and (x',y',z') in an earth-centered coordinate system, the vector from the first to the second is (dx,dy,dz) = (x'-x, y'-y, z'-z), whence the cosine of the angle made to the normal at (x,y,z) is the inner product of the unit length versions of those vectors:
Obtain its principal inverse cosine. Subtract this from 90 degrees if you want the angle of view relative to a nominal horizon. This is the "elevation."
Azimuth
A similar calculation obtains the local direction of view ("azimuth"). We need a level vector (u,v,w) that points due north. One such vector at the location (x,y,z) is (-zx, -zy, x^2+y^2). (The inner product of these two vectors is zero, proving it is level. Its projection onto the Equatorial plane is proportional to (-x,-y) which points directly inward, making it the projection of a north-pointing vector. These two calculations confirm that this is indeed the desired vector). Therefore
We also need the sine of the azimuth, which is similarly obtained once we know a vector pointing due East (locally). Such a vector is (-y, x, 0), because it clearly is perpendicular to (x,y,z) (the up direction) and the northern direction. Therefore
These values enable us to recover the azimuth as the inverse tangent of the cosine and sine.
Example
A pilot in an airplane flying west at 4000 meters, located at (lat, lon) = (39, -75), sees a jet far ahead located at (39, -76) and flying at 12000 meters. What is are the angles of view (relative to the level direction at the pilot's location)?
The XYZ coordinates of the airplanes are (x,y,z) = (1285410, -4797210, 3994830) and (x',y',z') = (1202990, -4824940, 3999870), respectively (in the ITRF00 datum, which uses the GRS80 ellipsoid). The pilot's view vector therefore is (dx,dy,dz) = (-82404.5, -27735.3, 5034.56). Applying the formula gives the cosine of the view angle as 0.0850706. Its inverse cosine is 85.1199 degrees, whence the elevation is 4.88009 degrees: the pilot is looking up by that much.
A north-pointing level vector is (-5.13499, 19.1641, 24.6655) (times 10^12) and an east-pointing level vector is (4.79721, 1.28541, 0) (times 10^6). Therefore, applying the last two formulas, the cosine of the azimuth equals 0.00575112 and its sine equals -0.996358. The ArcTangent function tells us the angle for the direction (0.00575112, -0.996358) is 270.331 degrees: almost due west. (It's not exactly west because the two planes lie on the same circle of latitude, which is curving toward the North pole: see Why is the 'straight line' path across continent so curved? for an extended explanation.)
By the way, this example confirms we got the orientation correct for the azimuth calculation: although it was clear the east-pointing vector was orthogonal to the other two vectors, until now it was not plain that it truly points east and not west.