Como posso usar o ArcGIS 10.1 para encontrar um ponto equidistante geodésico definido por três pontos?

Por exemplo, tenho coordenadas para três pontos de base em uma costa e preciso encontrar as coordenadas do ponto ao largo da costa, que são equidistantes das três. É um exercício simples de geometria, mas todas as medidas devem levar em consideração a geodésia.

Se eu estivesse abordando isso de uma maneira euclidiana, poderia medir os caminhos geodésicos que conectam os pontos de base, encontrar os pontos médios dos lados do triângulo resultante e criar ortodromos perpendiculares a cada um desses caminhos. Os três loxodromos provavelmente convergiriam no ponto equidistante. Se este for o método correto, deve haver uma maneira mais fácil de fazer isso no Arc.

Esta resposta está dividida em várias seções:

• Análise e Redução do Problema , mostrando como encontrar o ponto desejado com rotinas "enlatadas".

• Ilustração: um Protótipo de Trabalho , fornecendo código de trabalho.

• Exemplo , mostrando exemplos das soluções.

• Armadilhas , discutindo problemas em potencial e como lidar com eles.

• Implementação do ArcGIS , comentários sobre a criação de uma ferramenta personalizada do ArcGIS e onde obter as rotinas necessárias.

Análise e Redução do Problema

Vamos começar por observar que, no (perfeitamente redonda) modelo esférico haverá sempre ser uma solução --em verdade, exatamente duas soluções. Dados os pontos base A, B e C, cada par determina sua "bissetor perpendicular", que é o conjunto de pontos equidistantes dos dois pontos dados. Essa bissetriz é uma geodésica (grande círculo). A geometria esférica é elíptica : quaisquer duas geodésicas se cruzam (em dois pontos únicos). Assim, os pontos de interseção da bissetriz de AB e da bissetriz de BC são - por definição - equidistantes de A, B e C, resolvendo assim o problema. (Veja a primeira figura abaixo.)

As coisas parecem mais complicadas em um elipsóide, mas , como é uma pequena perturbação da esfera, podemos esperar um comportamento semelhante. (A análise disso nos levaria muito longe.) As fórmulas complicadas usadas (internamente em um SIG) para calcular distâncias precisas em um elipsóide não são uma complicação conceitual: o problema é basicamente o mesmo. Para ver como o problema é realmente simples, vamos declarar algo abstrato. Nesta declaração, "d (U, V)" refere-se à distância verdadeira e totalmente precisa entre os pontos U e V.

Dados três pontos A, B, C (como pares lat-lon) em um elipsóide, encontre um ponto X para o qual (1) d (X, A) = d (X, B) = d (X, C) e ( 2) essa distância comum é a menor possível.

Todas essas três distâncias dependem do X desconhecido . Assim, as diferenças nas distâncias u (X) = d (X, A) - d (X, B) e v (X) = d (X, B) - d (X, C) são funções com valor real de X. Novamente, de maneira um pouco abstrata, podemos reunir essas diferenças em um par ordenado. Também usaremos (lat, lon) como coordenadas para X, permitindo-nos considerá-lo também como um par ordenado, digamos X = (phi, lambda). Nesta configuração, a função

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

é uma função de uma parte de um espaço bidimensional que recebe valores no espaço bidimensional e nosso problema se reduz a

Encontre todo o possível (phi, lambda) para o qual F (phi, lambda) = (0,0).

Aqui é onde a abstração compensa: existe um grande número de ótimos softwares para resolver esse problema (busca de raiz multidimensional puramente numérica). O modo como funciona é que você escreve uma rotina para calcular F e depois a repassa para o software, juntamente com qualquer informação sobre restrições em sua entrada ( phi deve estar entre -90 e 90 graus e lambda deve estar entre -180 e 180 graus). Ele dispara por uma fração de segundo e retorna (normalmente) apenas um valor de ( phi , lambda ), se puder encontrar um.

Há detalhes para lidar, porque há uma arte nisso: existem vários métodos de solução para escolher, dependendo de como F "se comporta"; ajuda a "orientar" o software, fornecendo a ele um ponto de partida razoável para sua pesquisa (essa é uma maneira de obtermos a solução mais próxima , em vez de qualquer outra); e você geralmente precisa especificar o quão preciso você gostaria que a solução fosse (para que ele saiba quando interromper a pesquisa). (Para saber mais sobre o que os analistas de GIS precisam saber sobre esses detalhes, que surgem muito nos problemas de SIG, visite os tópicos de recomendação a serem incluídos em um curso de Ciência da computação para tecnologias geoespaciais e consulte a seção "Miscelânea" no final. )

Ilustração: um protótipo funcional

A análise mostra que precisamos programar duas coisas: uma estimativa inicial bruta da solução e o cálculo do próprio F.

A estimativa inicial pode ser feita por uma "média esférica" ​​dos três pontos base. Isso é obtido representando-as em coordenadas cartesianas geocêntricas (x, y, z), calculando a média dessas coordenadas e projetando essa média de volta para a esfera e reexpressando-a em latitude e longitude. O tamanho da esfera é imaterial e, portanto, os cálculos são feitos de maneira direta: como esse é apenas um ponto de partida, não precisamos de cálculos elipsoidais.

Para este protótipo de trabalho, usei o 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
  ]

(A Ifcondição final testa se a média pode falhar claramente ao indicar uma longitude; nesse caso, ela volta a uma média aritmética reta das latitudes e longitudes de sua entrada - talvez não seja uma ótima opção, mas pelo menos válida. Para aqueles que usam este código para orientação de implementação, observe que os argumentos do Mathematica ArcTan são revertidos em comparação com a maioria das outras implementações: seu primeiro argumento é a coordenada x, seu segundo é a coordenada y e retorna o ângulo feito pelo vetor ( x, y).)

Quanto à segunda parte, porque o Mathematica - como o ArcGIS e quase todos os outros GIS - contém código para calcular distâncias precisas no elipsóide, não há quase nada para escrever. Apenas chamamos a rotina de busca de raiz:

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
   ];

O aspecto mais notável dessa implementação é como ela evita a necessidade de restringir a latitude ( f) e a longitude ( q), sempre computando-as nos módulos 180 e 360 ​​graus, respectivamente. Isso evita ter que restringir o problema (que geralmente cria complicações). Os parâmetros de controle MaxIterationsetc. são ajustados para tornar esse código a maior precisão possível.

Para vê-lo em ação, vamos aplicá-lo aos três pontos base fornecidos em uma pergunta relacionada :

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

{-6,29692, 106,907}

As distâncias calculadas entre esta solução e os três pontos são

{1450.23206979, 1450.23206979, 1450.23206978}

(estes são metros). Eles concordam com o décimo primeiro dígito significativo (que é muito preciso, na verdade, uma vez que as distâncias raramente são precisas ou melhores que um milímetro). Aqui está uma figura desses três pontos (preto), seus três bissetores mútuos e a solução (vermelho):

Exemplo

Para testar esta implementação e entender melhor como o problema se comporta, aqui está um gráfico de contorno da discrepância quadrática média da raiz nas distâncias para três pontos base amplamente espaçados. (A discrepância do RMS é obtida calculando as três diferenças d (X, A) -d (X, B), d (X, B) -d (X, C) ed (X, C) -d (X , A) calculando a média de seus quadrados e obtendo a raiz quadrada.É igual a zero quando X resolve o problema e aumenta à medida que X se afasta de uma solução, e mede o quão "próximos" estamos de ser uma solução em qualquer local. )

Os pontos base (60, -120), (10, -40) e (45,10) são mostrados em vermelho nesta projeção do Plate Carree; a solução (49.2644488, -49.9052992) - que levou 0,03 segundos para calcular - está em amarelo. Sua discrepância RMS é inferior a três nanômetros , apesar de todas as distâncias relevantes serem milhares de quilômetros. As áreas escuras mostram valores pequenos do RMS e as áreas claras mostram valores altos.

Este mapa mostra claramente que outra solução está próxima (-49.2018206, 130.0297177) (calculada em um RMS de dois nanômetros, definindo o valor da pesquisa inicial diametralmente oposto à primeira solução).

Armadilhas

Instabilidade numérica

Quando os pontos de base são quase colineares e próximos, todas as soluções estarão a quase meio mundo de distância e extremamente difíceis de definir com precisão. A razão é que pequenas mudanças em um local em todo o mundo - movendo-o para ou longe dos pontos de base - induzirão apenas mudanças incrivelmente pequenas nas diferenças de distâncias. Apenas não há precisão e precisão suficientes incorporadas no cálculo usual das distâncias geodésicas para determinar o resultado.

Por exemplo, começando com os pontos base em (45.001, 0), (45, 0) e (44.999,0), que são separados ao longo do Meridiano Prime por apenas 111 metros entre cada par, triobtém a solução (11.8213, 77.745 ) As distâncias dele até os pontos de base são 8.127.964,998 77; 8.127.964,998 41; e 8.127.964,998 65 metros, respectivamente. Eles concordam com o milímetro mais próximo! Não tenho certeza de quão preciso seja esse resultado, mas não ficaria nem um pouco surpreso se outras implementações retornassem locais distantes deste, mostrando quase a mesma igualdade das três distâncias.

Tempo de computação

Esses cálculos, porque envolvem uma pesquisa considerável usando cálculos de distância complicados, não são rápidos, geralmente exigindo uma fração perceptível de segundo. Os aplicativos em tempo real precisam estar cientes disso.

Implementação do ArcGIS

Python é o ambiente de script preferido para o ArcGIS (começando com a versão 9). O pacote scipy.optimize possui um rootfinder multivariado rootque deve fazer o que FindRootfaz no código do Mathematica . É claro que o próprio ArcGIS oferece cálculos precisos de distância elipsoidal. O resto, então, são todos os detalhes da implementação: decida como as coordenadas do ponto base serão obtidas (de uma camada digitada pelo usuário? De um arquivo de texto? Do mouse?) E como a saída será apresentada (como coordenadas exibido na tela - como um ponto gráfico - como um novo objeto de ponto em uma camada?), escreva essa interface, transporta o código do Mathematica mostrado aqui (direto) e você estará pronto.

Como você observa, esse problema surge na determinação das fronteiras marítimas; geralmente é chamado de problema do "ponto triplo", e você pode pesquisar no Google e encontrar vários documentos sobre isso. Um desses documentos é meu (!) E ofereço uma solução precisa e rapidamente convergente. Consulte a Seção 14 de http://arxiv.org/abs/1102.1215

O método consiste nas seguintes etapas:

1. adivinhe um O de três pontos
2. use O como o centro de uma projeção equidistante azimutal
3. projeto A, B, C, para esta projeção
4. encontre o tri-ponto nesta projeção, O '
5. use O 'como o novo centro de projeção
6. repita até O 'e O coincidirem

A fórmula necessária para a solução de três pontos na projeção é apresentada no artigo. Contanto que você esteja usando uma projeção equidistante azimutal precisa, a resposta será exata. Convergência é um significado quadrático de que apenas algumas iterações são necessárias. Isso quase certamente superará os métodos gerais de busca de raiz sugeridos pelo @whuber.

Não posso ajudá-lo diretamente com o ArcGIS. Você pode pegar meu pacote python para fazer cálculos geodésicos em https://pypi.python.org/pypi/geographiclib e codificar a projeção com base nisso é simples.

Editar

O problema de encontrar o ponto triplo no caso degenerado do whuber (45 + eps, 0) (45,0) (45-eps, 0) foi considerado por Cayley em Nas linhas geodésicas de um esferóide oblato , Phil. Mag. (1870), http://books.google.com/books?id=4XGIOoCMYYAC&pg=PA15

Nesse caso, o tri-ponto é obtido seguindo um geodésico de (45,0) com o azimute 90 e encontrando o ponto em que a escala geodésica desaparece. Para o elipsóide WGS84, esse ponto é (-0.10690908732248, 89.89291072793167). A distância deste ponto para cada um de (45.001,0), (45,0), (44.999) é 10010287.665788943 m (dentro de um nanômetro ou mais). Isso é cerca de 1882 km a mais do que a estimativa da whuber (que apenas mostra o quão instável é este caso). Para uma terra esférica, o tri-ponto seria (0,90) ou (0, -90), é claro.

ADENDA: Aqui está uma implementação do método equidistante azimutal usando o 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

Testando isso usando o Octave, recebo

oitava: 1> formato longo
oitava: 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 o ponto triplo para Nova York, Tóquio e Wellington.

Este método é impreciso para os pontos colineares vizinhos, por exemplo, [45.001,0], [45,0], [44.999,0]. Nesse caso, você deve resolver M ₁₂ = 0 em uma geodésica que emana de [45,0] no azimute 90. A função a seguir faz o truque (usando o 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

Por exemplo, isso fornece:

[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

Eu estava curioso para ver com que rapidez a abordagem do @ cffk converge para uma solução, então escrevi um teste usando arcobjects, que produziu essa saída. As distâncias estão em 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

Aqui está o código fonte. (Editar) Alterou o FindCircleCenter para manipular interseções (pontos centrais) que caem da borda da projeção 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;
}

Há também uma abordagem alternativa na edição de junho de 2013 da Revista MSDN, Otimização do método Amoeba usando C # .

Editar

O código publicado anteriormente convergiu para o antípode em alguns casos. Alterei o código para produzir essa saída para os pontos de teste do @ cffk.

Aqui está a saída que agora produz:

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

Aqui está o 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

Aqui estão os resultados que eu recebo com 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