¿Cómo puedo utilizar ArcGIS 10.1 para encontrar un punto geodésico equidistante definido por tres puntos?

Question

Por ejemplo, tengo las coordenadas de tres puntos base en una costa y necesito hallar las coordenadas del punto de la costa que equidista de los tres. Es un simple ejercicio de geometría, pero todas las mediciones deben tener en cuenta la geodesia.

Si lo enfocara de manera euclidiana, podría medir las trayectorias geodésicas que conectan los puntos base, hallar los puntos medios de los lados del triángulo resultante y crear ortodrómicas perpendiculares a cada una de esas trayectorias. Los tres loxodrómicas convergerían presumiblemente en el punto equidistante. Si este es el método correcto, tiene que haber una manera más fácil de hacerlo en Arc.

Answer 1

Esta respuesta está dividida en varias secciones:

• Análisis y reducción del problema , mostrando cómo encontrar el punto deseado con rutinas "enlatadas".

• Ilustración: un prototipo funcional , dando código de trabajo.

• Ejemplo con ejemplos de soluciones.

• Escollos En la reunión se debaten los posibles problemas y la forma de resolverlos.

• Implantación de ArcGIS , comentarios sobre la creación de una herramienta ArcGIS personalizada y dónde obtener las rutinas necesarias.

---

Análisis y reducción del problema

Comencemos observando que en el modelo esférico (perfectamente redondo) habrá siempre sea una solución --de hecho, exactamente dos soluciones. Dados los puntos base A, B y C, cada par determina su "bisectriz perpendicular", que es el conjunto de puntos equidistantes de los dos puntos dados. Esta bisectriz es una geodésica (gran círculo). La geometría esférica es elíptica dos geodésicas cualesquiera se intersecan (en dos únicos puntos). Así, los puntos de intersección de la bisectriz de AB y la bisectriz de BC son -por definición- equidistantes de A, B y C, con lo que se resuelve el problema. (Véase la primera figura a continuación).

Las cosas parecen más complicadas en un elipsoide, pero como se trata de una pequeña perturbación de la esfera, podemos esperar un comportamiento similar. (Sin embargo, las complicadas fórmulas utilizadas (internamente en un SIG) para calcular distancias precisas en un elipsoide no son una complicación conceptual: el problema es básicamente el mismo. Para ver lo sencillo que es el problema en realidad, Digámoslo de forma algo abstracta. En este enunciado, "d(U,V)" se refiere a la distancia verdadera, totalmente exacta, entre los puntos U y V.

Dados tres puntos A, B, C (como pares latitud-longitud) en un elipsoide, encuentre un punto X para el que (1) d(X,A) = d(X,B) = d(X,C) y (2) esta distancia común sea lo más pequeña posible.

Estas tres distancias dependen de la incógnita X . Así, el diferencias en distancias u(X) = d(X,A) - d(X,B) y v(X) = d(X,B) - d(X,C) son funciones de valor real de X. De nuevo, de forma algo abstracta, podemos ensamblar estas diferencias en un par ordenado. También utilizaremos (lat, lon) como coordenadas para X, lo que nos permite considerarlo también como un par ordenado, digamos X = (phi, lambda). En esta configuración, la función

F(phi, lambda) = (u(X), v(X))

es una función de una porción de un espacio bidimensional que toma valores en el espacio bidimensional y nuestro problema se reduce a

Encuentra todos los posibles (phi, lambda) para los que F(phi, lambda) = (0,0).

Aquí es donde la abstracción merece la pena: Existe un gran número de programas informáticos para resolver este problema (puramente numérico y multidimensional de búsqueda de raíces). El funcionamiento consiste en escribir una rutina para calcular F y, a continuación, lo transmite al programa junto con la información sobre las restricciones de entrada ( phi debe estar entre -90 y 90 grados y lambda debe estar entre -180 y 180 grados). Se aleja durante una fracción de segundo y vuelve (normalmente) justo un valor de ( phi , lambda ), si lo encuentra.

Hay detalles de los que ocuparse, porque hay un arte en esto: hay varios métodos de solución para elegir, dependiendo de cómo F se "comporta"; ayuda a "dirigir" el software dándole un punto de partida razonable para su búsqueda (es una forma de obtener el más cercano en lugar de cualquier otra); y normalmente hay que especificar la precisión que se desea que tenga la solución (para que sepa cuándo detener la búsqueda). (Para saber más sobre lo que necesitan saber los analistas de SIG acerca de estos detalles, que aparecen mucho en los problemas de SIG, visite Recomendar temas para un curso de Informática para las Tecnologías Geoespaciales y busque en la sección "Miscelánea" cerca del final).

---

Ilustración: un prototipo funcional

El análisis muestra que necesitamos programar dos cosas: una estimación inicial aproximada de la solución y el cálculo de F sí mismo.

La estimación inicial puede realizarse mediante una "media esférica" de los tres puntos base. Se obtiene representándolos en coordenadas cartesianas geocéntricas (x,y,z), haciendo la media de esas coordenadas y proyectando esa media de nuevo a la esfera y volviéndola a expresar en latitud y longitud. El tamaño de la esfera es irrelevante, por lo que los cálculos resultan más sencillos: al tratarse sólo de un punto de partida, no necesitamos cálculos elipsoidales.

Para este prototipo de trabajo he utilizado Mathematica 8.

sphericalMean[points_] := Module[{sToC, cToS, cMean},
  sToC[{f_, l_}] := {Cos[f] Cos[l], Cos[f] Sin[l], Sin[f]};
  cToS[{x_, y_, z_}] := {ArcTan[x, y], ArcTan[Norm[{x, y}], z]};
  cMean = Mean[sToC /@ (points Degree)];
  If[Norm[Most@cMean] < 10^(-8), Mean[points], cToS[cMean]] / Degree
  ]

(El final If comprueba si la media puede no indicar claramente una longitud; si es así, vuelve a una media aritmética directa de las latitudes y longitudes de su entrada; puede que no sea una gran elección, pero al menos es válida. Para quienes utilicen este código como guía de implementación, tenga en cuenta que los argumentos de Mathematica ArcTan se invierten en comparación con la mayoría de las demás implementaciones: su primer argumento es la coordenada x, el segundo es la coordenada y, y devuelve el ángulo formado por el vector (x,y)).

En cuanto a la segunda parte, porque Mathematica --como ArcGIS y casi todos los demás SIG-- contiene código para calcular distancias precisas en el elipsoide, no hay casi nada que escribir. Nos limitamos a llamar a la rutina de búsqueda de raíces:

tri[a_, b_, c_] := Block[{d = sphericalMean[{a, b, c}], sol, f, q},
   sol = FindRoot[{GeoDistance[{Mod[f, 180, -90], Mod[q, 360, -180]}, a] == 
                   GeoDistance[{Mod[f, 180, -90], Mod[q, 360, -180]}, b] ==
                   GeoDistance[{Mod[f, 180, -90], Mod[q, 360, -180]}, c]}, 
           {{f, d[[1]]}, {q, d[[2]]}}, 
           MaxIterations -> 1000, AccuracyGoal -> Infinity, PrecisionGoal -> 8];
   {Mod[f, 180, -90], Mod[q, 360, -180]} /. sol
   ];

El aspecto más destacable de esta implementación es cómo evita la necesidad de restringir la latitud ( f ) y la longitud ( q ) calculándolos siempre en módulo de 180 y 360 grados, respectivamente. Esto evita tener que restringir el problema (lo que a menudo crea complicaciones). Los parámetros de control MaxIterations etc. se ajustan para que este código proporcione la mayor precisión posible.

Para verlo en acción, apliquémoslo a los tres puntos base dados en un pregunta relacionada :

sol = tri @@ (bases = {{-6.28530175, 106.9004975375}, {-6.28955287, 106.89573839}, {-6.28388865789474, 106.908087643421}})

{-6.29692, 106.907}

Las distancias calculadas entre esta solución y los tres puntos son

{1450.23206979, 1450.23206979, 1450.23206978}

(son contadores). Coinciden hasta el undécimo dígito significativo (lo cual es demasiado preciso, en realidad, ya que las distancias rara vez tienen una precisión superior a un milímetro, más o menos). Aquí tienes una imagen de estos tres puntos (en negro), sus tres bisectrices mutuas y la solución (en rojo):

---

Ejemplo

Para probar esta implementación y comprender mejor cómo se comporta el problema, he aquí un gráfico de contorno de la discrepancia cuadrática media de las distancias para tres puntos base muy espaciados. (La discrepancia RMS se obtiene calculando las tres diferencias d(X,A)-d(X,B), d(X,B)-d(X,C), y d(X,C)-d(X,A), promediando sus cuadrados y sacando la raíz cuadrada. Es igual a cero cuando X resuelve el problema y, en caso contrario, aumenta a medida que X se aleja de una solución, por lo que mide lo "cerca" que estamos de ser una solución en cualquier punto).

Los puntos base (60,-120), (10,-40) y (45,10) se muestran en rojo en esta proyección de Plate Carree; la solución (49,2644488, -49,9052992), que se calculó en 0,03 segundos, aparece en amarillo. Su discrepancia RMS es inferior a tres nanómetros a pesar de que todas las distancias relevantes son de miles de kilómetros. Las zonas oscuras muestran valores pequeños de la RMS y las claras, valores altos.

Este mapa muestra claramente que otra solución se encuentra cerca de (-49.2018206, 130.0297177) (calculada a un RMS de dos nanómetros estableciendo el valor de búsqueda inicial diametralmente opuesto a la primera solución).

---

Escollos

Inestabilidad numérica

Cuando los puntos base son casi colineales y están próximos entre sí, todas las soluciones estarán a casi medio mundo de distancia y serán extremadamente difíciles de precisar. La razón es que pequeños cambios en una ubicación a lo largo del globo -alejándola o acercándola a los puntos base- sólo inducirán cambios increíblemente diminutos en las diferencias de distancias. El cálculo habitual de las distancias geodésicas carece de la exactitud y precisión necesarias para precisar el resultado.

Por ejemplo, comenzando con los puntos base en (45,001,0), (45,0) y (44,999,0), que están separados a lo largo del Primer Meridiano por sólo 111 metros entre cada par, tri obtiene la solución (11,8213, 77,745). Las distancias desde ella a los puntos base son 8.127.964,998 77; 8.127.964,998 41; y 8.127.964,998 65 metros, respectivamente. Coinciden al milímetro. No estoy seguro de la exactitud de este resultado, pero no me extrañaría lo más mínimo que otras implementaciones devolvieran ubicaciones muy alejadas de ésta, mostrando una igualdad casi tan buena de las tres distancias.

Tiempo de cálculo

Estos cálculos, al implicar una búsqueda considerable mediante complicados cálculos de distancias, no son rápidos y suelen requerir una fracción de segundo considerable. Las aplicaciones en tiempo real deben ser conscientes de ello.

---

Implantación de ArcGIS

Python es el entorno de scripting preferido para ArcGIS (a partir de la versión 9). El sitio paquete scipy.optimize tiene un buscador de raíces multivariante root que debe hacer lo que FindRoot hace en el Mathematica código. Por supuesto, el propio ArcGIS ofrece cálculos precisos de distancias elipsoidales. El resto, por tanto, son detalles de implementación: decidir cómo se obtendrán las coordenadas del punto base (¿de una capa? ¿introducidas por el usuario? ¿de un archivo de texto? ¿del ratón?) y cómo se presentará la salida (¿como coordenadas mostradas en pantalla? ¿como un punto gráfico? ¿como un nuevo objeto punto en una capa?), escribir esa interfaz, portar el programa a la interfaz de ArcGIS y ejecutar el programa. Mathematica código que se muestra aquí (sencillo), y ya estará todo listo.

Answer 2

Como usted señala, este problema se plantea a la hora de determinar las fronteras marítimas; a menudo se denomina el problema de los "tres puntos" y se puede buscar en Google y encontrar varios artículos que lo tratan. Uno de ellos es mío y ofrezco una solución precisa y de rápida convergencia. Véase Sección 14 de http://arxiv.org/abs/1102.1215

El método consta de los siguientes pasos:

1. adivinar un tri-punto O
2. utilizar O como centro de una proyección equidistante acimutal
3. proyectar A, B, C, a esta proyección
4. encontrar el tri-punto en esta proyección, O'
5. utilizar O' como nuevo centro de proyección
6. repetir hasta que O' y O coincidan

La fórmula necesaria para la solución de tres puntos en la proyección es en el documento. Si se utiliza una proyección acimutal equidistante equidistante, la respuesta será exacta. La convergencia es cuadrática, lo que significa que sólo se necesitan unas pocas iteraciones. Esto casi seguro que superan los métodos generales de búsqueda de raíces sugerido por @whuber.

No puedo ayudarte directamente con ArcGIS. Usted puede agarrar mi paquete python para hacer cálculos geodésicos de https://pypi.python.org/pypi/geographiclib y codificar la proyección basándose en esto es sencillo.

---

Editar

El problema de encontrar el tri-punto en el caso degenerado de @whuber (45+eps,0) (45,0) (45-eps,0) fue considerado por Cayley en En el líneas geodésicas en un esferoide oblato Phil. Mag. (1870), http://books.google.com/books?id=4XGIOoCMYYAC&pg=PA15

En este caso, el punto triple se obtiene siguiendo una geodésica desde (45,0) con acimut 90 y encontrando el punto en el que la escala geodésica se desvanece. Para el elipsoide WGS84, este punto es (-0,106908732248, 89.89291072793167). La distancia desde este punto a cada uno de (45,001,0), (45,0), (44,999) es de 10010287,665788943 m (dentro de un nanómetro más o menos). Esto es alrededor de 1882 km más que la estimación de whuber (que sólo va a lo inestable que es este caso). Para una Tierra esférica, el punto triple sería (0,90) o (0,-90), por supuesto.

ADDENDUM: He aquí una implementación del método equidistante azimutal utilizando Matlab

function [lat, lon] = tripoint(lat1, lon1, lat2, lon2, lat3, lon3)
% Compute point equidistant from arguments
% Requires:
%   http://www.mathworks.com/matlabcentral/fileexchange/39108
%   http://www.mathworks.com/matlabcentral/fileexchange/39366
  lats = [lat1, lat2, lat3];
  lons = [lon1, lon2, lon3];
  lat0 = lat1;  lon0 = lon1; % feeble guess for tri point
  for i = 1:6
    [x, y] = eqdazim_fwd(lat0, lon0, lats, lons);
    a = [x(1), y(1), 0];
    b = [x(2), y(2), 0];
    c = [x(3), y(3), 0];
    z = [0, 0, 1];
    % Eq. (97) of http://arxiv.org/abs/1102.1215
    o = cross((a*a') * (b - c) + (b*b') * (c - a) + (c*c') * (a - b), z) ...
        / (2 * dot(cross(a - b, b - c), z));
    [lat0, lon0] = eqdazim_inv(lat0, lon0, o(1), o(2))
  end
  % optional check
  s12 = geoddistance(lat0, lon0, lats, lons); ds12 = max(s12) - min(s12)
  lat = lat0; lon = lon0;
end

Probando esto con Octave obtengo

octave:1> format long
octave:2> \[lat0,lon0\] = tripoint(41,-74,36,140,-41,175)
lat0 =  15.4151378380375
lon0 = -162.479314381144
lat0 =  15.9969703299812
lon0 = -147.046790722192
lat0 =  16.2232960167545
lon0 = -147.157646039471
lat0 =  16.2233394851560
lon0 = -147.157748279290
lat0 =  16.2233394851809
lon0 = -147.157748279312
lat0 =  16.2233394851809
lon0 = -147.157748279312
ds12 =  3.72529029846191e-09
lat0 =  16.2233394851809
lon0 = -147.157748279312

como punto de enlace con Nueva York, Tokio y Wellington.

Este método es inexacto para puntos colineales vecinos, por ejemplo, [45,001,0], [45,0], [44,999,0]. En ese caso, debes resolver para M ₁₂ \= 0 en una geodésica que emana de [45,0] en el acimut 90. La siguiente función lo consigue (utilizando el método de Newton):

function [lat2,lon2] = semiconj(lat1, lon1, azi1)
% Find the point where neighboring parallel geodesics emanating from
% close to [lat1, lon1] with azimuth azi1 intersect.

  % First guess is 90 deg on aux sphere
  [lat2, lon2, ~, ~, m12, M12, M21, s12] = ...
      geodreckon(lat1, lon1, 90, azi1, defaultellipsoid(), true);
  M12
  % dM12/ds2 = - (1 - M12*M21/m12)
  for i = 1:3
    s12 = s12 - M12 / ( -(1 - M12*M21)/m12 ); % Newton
    [lat2, lon2, ~, ~, m12, M12, M21] = geodreckon(lat1, lon1, s12, azi1);
    M12
  end
end

Para el ejemplo, esto da:

\[lat2,lon2\] = semiconj(45, 0, 90)
M12 =  0.00262997817649321
M12 = -6.08402492665097e-09
M12 =  4.38017677684144e-17
M12 =  4.38017677684144e-17
lat2 = -0.106909087322479
lon2 =  89.8929107279317

Answer 3

Tenía curiosidad por ver lo rápido que el enfoque de @cffk converge en una solución, así que escribí una prueba utilizando arcobjects, que produjo este resultado. Las distancias están en metros:

0 longitude: 0 latitude: 90
    Distances: 3134.05443974188 2844.67237777542 3234.33025754997
    Diffs: 289.382061966458 -389.657879774548 -100.27581780809
1 longitude: 106.906152157596 latitude: -6.31307123035178
    Distances: 1450.23208989615 1450.23208089398 1450.23209429293
    Diffs: 9.00216559784894E-06 -1.33989510686661E-05 -4.39678547081712E-06
2 longitude: 106.906583669013 latitude: -6.29691590176649
    Distances: 1450.23206976414 1450.23206976408 1450.23206976433
    Diffs: 6.18456397205591E-11 -2.47382558882236E-10 -1.85536919161677E-10
3 longitude: 106.906583669041 latitude: -6.29691590154641
    Distances: 1450.23206976438 1450.23206976423 1450.23206976459
    Diffs: 1.47565515362658E-10 -3.61751517630182E-10 -2.14186002267525E-10
4 longitude: 106.906583669041 latitude: -6.29691590154641
    Distances: 1450.23206976438 1450.23206976423 1450.23206976459
    Diffs: 1.47565515362658E-10 -3.61751517630182E-10 -2.14186002267525E-10

Aquí está el código fuente. (Edit) Cambiado el FindCircleCenter para manejar las intersecciones (puntos centrales) que caen fuera del borde de la proyección azimutal:

public static void Test()
{
    var t = Type.GetTypeFromProgID("esriGeometry.SpatialReferenceEnvironment");
    var srf = Activator.CreateInstance(t) as ISpatialReferenceFactory2;
    var pcs = srf.CreateProjectedCoordinateSystem((int)esriSRProjCSType.esriSRProjCS_WGS1984N_PoleAziEqui)
        as IProjectedCoordinateSystem2;

    var pntA = MakePoint(106.9004975375, -6.28530175, pcs.GeographicCoordinateSystem);
    var pntB = MakePoint(106.89573839, -6.28955287, pcs.GeographicCoordinateSystem);
    var pntC = MakePoint(106.908087643421, -6.28388865789474, pcs.GeographicCoordinateSystem);

    int maxIter = 5;
    for (int i = 0; i < maxIter; i++)
    {
        var msg = string.Format("{0} longitude: {1} latitude: {2}", i, pcs.get_CentralMeridian(true), pcs.LatitudeOfOrigin);
        Debug.Print(msg);
        var newCenter = FindCircleCenter(ProjectClone(pntA, pcs), ProjectClone(pntB, pcs), ProjectClone(pntC, pcs));
        newCenter.Project(pcs.GeographicCoordinateSystem); // unproject
        var distA = GetGeodesicDistance(newCenter, pntA);
        var distB = GetGeodesicDistance(newCenter, pntB);
        var distC = GetGeodesicDistance(newCenter, pntC);
        Debug.Print("\tDistances: {0} {1} {2}", distA, distB, distC);
        var diffAB = distA - distB;
        var diffBC = distB - distC;
        var diffAC = distA - distC;
        Debug.Print("\tDiffs: {0} {1} {2}", diffAB, diffBC, diffAC);

        pcs.set_CentralMeridian(true, newCenter.X);
        pcs.LatitudeOfOrigin = newCenter.Y;
    }
}
public static IPoint FindCircleCenter(IPoint a, IPoint b, IPoint c)
{
    // from http://blog.csharphelper.com/2011/11/08/draw-a-circle-through-three-points-in-c.aspx
    // Get the perpendicular bisector of (x1, y1) and (x2, y2).
    var x1 = (b.X + a.X) / 2;
    var y1 = (b.Y + a.Y) / 2;
    var dy1 = b.X - a.X;
    var dx1 = -(b.Y - a.Y);

    // Get the perpendicular bisector of (x2, y2) and (x3, y3).
    var x2 = (c.X + b.X) / 2;
    var y2 = (c.Y + b.Y) / 2;
    var dy2 = c.X - b.X;
    var dx2 = -(c.Y - b.Y);

    // See where the lines intersect.
    var cx = (y1 * dx1 * dx2 + x2 * dx1 * dy2 - x1 * dy1 * dx2 - y2 * dx1 * dx2)
        / (dx1 * dy2 - dy1 * dx2);
    var cy = (cx - x1) * dy1 / dx1 + y1;

    // make sure the intersection point falls
    // within the projection.
    var earthRadius = ((IProjectedCoordinateSystem)a.SpatialReference).GeographicCoordinateSystem.Datum.Spheroid.SemiMinorAxis;

    // distance is from center of projection
    var dist = Math.Sqrt((cx * cx) + (cy * cy));
    double factor = 1.0;
    if (dist > earthRadius * Math.PI)
    {
        // apply a factor so we don't fall off the edge
        // of the projection
        factor = earthRadius / dist;
    }
    var outPoint = new PointClass() as IPoint;
    outPoint.PutCoords(cx * factor, cy* factor);
    outPoint.SpatialReference = a.SpatialReference;
    return outPoint;
}

public static double GetGeodesicDistance(IPoint pnt1, IPoint pnt2)
{
    var pc = new PolylineClass() as IPointCollection;
    var gcs = pnt1.SpatialReference as IGeographicCoordinateSystem;
    if (gcs == null)
        throw new Exception("point does not have a gcs");
    ((IGeometry)pc).SpatialReference = gcs;
    pc.AddPoint(pnt1);
    pc.AddPoint(pnt2);
    var t = Type.GetTypeFromProgID("esriGeometry.SpatialReferenceEnvironment");
    var srf = Activator.CreateInstance(t) as ISpatialReferenceFactory2;
    var unit = srf.CreateUnit((int)esriSRUnitType.esriSRUnit_Meter) as ILinearUnit;
    var pcGeodetic = pc as IPolycurveGeodetic;
    return pcGeodetic.get_LengthGeodetic(esriGeodeticType.esriGeodeticTypeGeodesic, unit);
}

public static IPoint ProjectClone(IPoint pnt, ISpatialReference sr)
{
    var clone = ((IClone)pnt).Clone() as IPoint;
    clone.Project(sr);
    return clone;
}

public static IPoint MakePoint(double longitude, double latitude, ISpatialReference sr)
{
    var pnt = new PointClass() as IPoint;
    pnt.PutCoords(longitude, latitude);
    pnt.SpatialReference = sr;
    return pnt;
}

También hay un enfoque alternativo en la edición de junio de 2013 de MSDN Magazine, Optimización del método de la ameba con C# .

---

Editar

El código publicado anteriormente convergía a la antípoda en algunos casos. He modificado el código para que produzca este resultado para los puntos de prueba de @cffk.

Esta es la salida que produce ahora:

0 0
0 longitude: 0 latitude: 0
    MaxDiff: 1859074.90170379 Distances: 13541157.6493561 11682082.7476523 11863320.2116807
1 longitude: 43.5318402621384 latitude: -17.1167429904981
    MaxDiff: 21796.9793742411 Distances: 12584188.7592282 12588146.4851222 12566349.505748
2 longitude: 32.8331167578493 latitude: -16.2707976739314
    MaxDiff: 6.05585224926472 Distances: 12577536.3369782 12577541.3560203 12577542.3928305
3 longitude: 32.8623898057665 latitude: -16.1374156408507
    MaxDiff: 5.58793544769287E-07 Distances: 12577539.6118671 12577539.6118666 12577539.6118669
4 longitude: -147.137582018133 latitude: 16.1374288796667
    MaxDiff: 1.12284109462053 Distances: 7441375.08265703 7441376.12671342 7441376.20549812
5 longitude: -147.157742373074 latitude: 16.2233413614432
    MaxDiff: 7.45058059692383E-09 Distances: 7441375.70752843 7441375.70752842 7441375.70752842
5 longitude: -147.157742373074 latitude: 16.2233413614432 Distance 7441375.70752843
iterations: 5

Aquí está el código revisado:

class Program
{
    private static LicenseInitializer m_AOLicenseInitializer = new tripoint.LicenseInitializer();

    [STAThread()]
    static void Main(string[] args)
    {
        //ESRI License Initializer generated code.
        m_AOLicenseInitializer.InitializeApplication(new esriLicenseProductCode[] { esriLicenseProductCode.esriLicenseProductCodeStandard },
        new esriLicenseExtensionCode[] { });
        try
        {
            var t = Type.GetTypeFromProgID("esriGeometry.SpatialReferenceEnvironment");
            var srf = Activator.CreateInstance(t) as ISpatialReferenceFactory2;
            var pcs = srf.CreateProjectedCoordinateSystem((int)esriSRProjCSType.esriSRProjCS_World_AzimuthalEquidistant)
                as IProjectedCoordinateSystem2;
            Debug.Print("{0} {1}", pcs.get_CentralMeridian(true), pcs.LatitudeOfOrigin);
            int max = int.MinValue;
            for (int i = 0; i < 1; i++)
            {
                var iterations = Test(pcs);
                max = Math.Max(max, iterations);
                Debug.Print("iterations: {0}", iterations);
            }
            Debug.Print("max number of iterations: {0}", max);
        }
        catch (Exception ex)
        {
            Debug.Print(ex.Message);
            Debug.Print(ex.StackTrace);
        }
        //ESRI License Initializer generated code.
        //Do not make any call to ArcObjects after ShutDownApplication()
        m_AOLicenseInitializer.ShutdownApplication();
    }
    public static int Test(IProjectedCoordinateSystem2 pcs)
    {
        var pntA = MakePoint(-74.0, 41.0, pcs.GeographicCoordinateSystem);
        var pntB = MakePoint(140.0, 36.0, pcs.GeographicCoordinateSystem);
        var pntC = MakePoint(175.0, -41.0, pcs.GeographicCoordinateSystem);

        //var r = new Random();
        //var pntA = MakeRandomPoint(r, pcs.GeographicCoordinateSystem);
        //var pntB = MakeRandomPoint(r, pcs.GeographicCoordinateSystem);
        //var pntC = MakeRandomPoint(r, pcs.GeographicCoordinateSystem);

        int maxIterations = 100;
        for (int i = 0; i < maxIterations; i++)
        {
            var msg = string.Format("{0} longitude: {1} latitude: {2}", i, pcs.get_CentralMeridian(true), pcs.LatitudeOfOrigin);
            Debug.Print(msg);
            var newCenter = FindCircleCenter(ProjectClone(pntA, pcs), ProjectClone(pntB, pcs), ProjectClone(pntC, pcs));
            var c = ((IClone)newCenter).Clone() as IPoint;
            newCenter.Project(pcs.GeographicCoordinateSystem); // unproject
            //newCenter = MakePoint(-147.1577482, 16.2233394, pcs.GeographicCoordinateSystem);
            var distA = GetGeodesicDistance(newCenter, pntA);
            var distB = GetGeodesicDistance(newCenter, pntB);
            var distC = GetGeodesicDistance(newCenter, pntC);
            var diffAB = Math.Abs(distA - distB);
            var diffBC = Math.Abs(distB - distC);
            var diffAC = Math.Abs(distA - distC);
            var maxDiff = GetMax(new double[] {diffAB,diffAC,diffBC});
            Debug.Print("\tMaxDiff: {0} Distances: {1} {2} {3}",maxDiff, distA, distB, distC);
            if (maxDiff < 0.000001)
            {
                var earthRadius = pcs.GeographicCoordinateSystem.Datum.Spheroid.SemiMinorAxis;
                if (distA > earthRadius * Math.PI / 2.0)
                {
                    newCenter = AntiPode(newCenter);
                }
                else
                {
                    Debug.Print("{0} longitude: {1} latitude: {2} Distance {3}", i, pcs.get_CentralMeridian(true), pcs.LatitudeOfOrigin, distA);
                    return i;
                }
            }
            //Debug.Print("\tDiffs: {0} {1} {2}", diffAB, diffBC, diffAC);

            pcs.set_CentralMeridian(true, newCenter.X);
            pcs.LatitudeOfOrigin = newCenter.Y;
        }
        return maxIterations;
    }

    public static IPoint FindCircleCenter(IPoint a, IPoint b, IPoint c)
    {
        // from http://blog.csharphelper.com/2011/11/08/draw-a-circle-through-three-points-in-c.aspx
        // Get the perpendicular bisector of (x1, y1) and (x2, y2).
        var x1 = (b.X + a.X) / 2;
        var y1 = (b.Y + a.Y) / 2;
        var dy1 = b.X - a.X;
        var dx1 = -(b.Y - a.Y);

        // Get the perpendicular bisector of (x2, y2) and (x3, y3).
        var x2 = (c.X + b.X) / 2;
        var y2 = (c.Y + b.Y) / 2;
        var dy2 = c.X - b.X;
        var dx2 = -(c.Y - b.Y);

        // See where the lines intersect.
        var cx = (y1 * dx1 * dx2 + x2 * dx1 * dy2 - x1 * dy1 * dx2 - y2 * dx1 * dx2)
            / (dx1 * dy2 - dy1 * dx2);
        var cy = (cx - x1) * dy1 / dx1 + y1;

        // make sure the intersection point falls
        // within the projection.
        var earthRadius = ((IProjectedCoordinateSystem)a.SpatialReference).GeographicCoordinateSystem.Datum.Spheroid.SemiMinorAxis;

        // distance is from center of projection
        var dist = Math.Sqrt((cx * cx) + (cy * cy));
        double factor = 1.0;
        if (dist > earthRadius * Math.PI)
        {
            // apply a factor so we don't fall off the edge
            // of the projection
            factor = earthRadius / dist;
        }
        var outPoint = new PointClass() as IPoint;
        outPoint.PutCoords(cx * factor, cy* factor);
        outPoint.SpatialReference = a.SpatialReference;
        return outPoint;
    }

    public static IPoint AntiPode(IPoint pnt)
    {
        if (!(pnt.SpatialReference is IGeographicCoordinateSystem))
            throw new Exception("antipode of non-gcs projection not supported");
        var outPnt = new PointClass() as IPoint;
        outPnt.SpatialReference = pnt.SpatialReference;
        if (pnt.X < 0.0)
            outPnt.X = 180.0 + pnt.X;
        else
            outPnt.X = pnt.X - 180.0;
        outPnt.Y = -pnt.Y;
        return outPnt;
    }

    public static IPoint MakeRandomPoint(Random r, IGeographicCoordinateSystem gcs)
    {
        var latitude = (r.NextDouble() - 0.5) * 180.0;
        var longitude = (r.NextDouble() - 0.5) * 360.0;
        //Debug.Print("{0} {1}", latitude, longitude);
        return MakePoint(longitude, latitude, gcs);
    }
    public static double GetMax(double[] dbls)
    {
        var max = double.MinValue;
        foreach (var d in dbls)
        {
            if (d > max)
                max = d;
        }
        return max;
    }
    public static IPoint MakePoint(IPoint[] pnts)
    {
        double sumx = 0.0;
        double sumy = 0.0;
        foreach (var pnt in pnts)
        {
            sumx += pnt.X;
            sumy += pnt.Y;
        }
        return MakePoint(sumx / pnts.Length, sumy / pnts.Length, pnts[0].SpatialReference);
    }
    public static double GetGeodesicDistance(IPoint pnt1, IPoint pnt2)
    {
        var pc = new PolylineClass() as IPointCollection;
        var gcs = pnt1.SpatialReference as IGeographicCoordinateSystem;
        if (gcs == null)
            throw new Exception("point does not have a gcs");
        ((IGeometry)pc).SpatialReference = gcs;
        pc.AddPoint(pnt1);
        pc.AddPoint(pnt2);
        var t = Type.GetTypeFromProgID("esriGeometry.SpatialReferenceEnvironment");
        var srf = Activator.CreateInstance(t) as ISpatialReferenceFactory2;
        var unit = srf.CreateUnit((int)esriSRUnitType.esriSRUnit_Meter) as ILinearUnit;
        var pcGeodetic = pc as IPolycurveGeodetic;
        return pcGeodetic.get_LengthGeodetic(esriGeodeticType.esriGeodeticTypeGeodesic, unit);
    }

    public static IPoint ProjectClone(IPoint pnt, ISpatialReference sr)
    {
        var clone = ((IClone)pnt).Clone() as IPoint;
        clone.Project(sr);
        return clone;
    }

    public static IPoint MakePoint(double longitude, double latitude, ISpatialReference sr)
    {
        var pnt = new PointClass() as IPoint;
        pnt.PutCoords(longitude, latitude);
        pnt.SpatialReference = sr;
        return pnt;
    }
}

---

Editar

Estos son los resultados que obtengo con esriSRProjCS_WGS1984N_PoleAziEqui

0 90
0 longitude: 0 latitude: 90
    MaxDiff: 1275775.91880553 Distances: 8003451.67666723 7797996.2370572 6727675.7578617
1 longitude: -148.003774863594 latitude: 9.20238223616225
    MaxDiff: 14487.6784785809 Distances: 7439006.46128994 7432752.45732905 7447240.13580763
2 longitude: -147.197808459106 latitude: 16.3073233548167
    MaxDiff: 2.32572609744966 Distances: 7441374.94409209 7441377.26981819 7441375.90768183
3 longitude: -147.157734641831 latitude: 16.2233338760411
    MaxDiff: 7.72997736930847E-08 Distances: 7441375.70752842 7441375.70752848 7441375.7075284
3 longitude: -147.157734641831 latitude: 16.2233338760411 Distance 7441375.70752842