Looking closely at your data it looks like in the distance calculations you have used different coordinates from those in the *activity_objects* table.
In example on the first query you wrote:
[...] ST_GeographyFromText('SRID=4326;POINT(-70.01 15.01)'));
where it should have been:
[...] ST_GeographyFromText('SRID=4326;POINT(-73.01 15.01)'));
Notice -70.01 being used in place of -73.01.
hence the discrepancies in all your subsequent calculations.
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
If the data you will be using is on the scale of the given example, you can treat the problem and Earth as flat without significant accuracy loss. If that does not violate your use case, then calculate the metric shift for both dimensions like you had a triangle and convert it back to degrees. The only tricky part is that the longitude degrees are not of constant length (narrowing from the equator to the poles) and you have to multiply them by the cosine of the (freshly computed) latitude to compensate before using them as a conversion factor. You end up with the coordinate shift in degrees, which you then just add to the first point's coordinates.