Bagaimana saya bisa menggunakan ArcGIS 10.1 untuk menemukan titik berjarak sama geodesik yang didefinisikan oleh tiga titik?

Sebagai contoh, saya memiliki koordinat untuk tiga titik dasar pada garis pantai dan saya perlu menemukan koordinat titik di lepas pantai yang berjarak sama dari ketiga titik tersebut. Ini adalah latihan sederhana dalam geometri, tetapi semua pengukuran harus memperhitungkan geodesi.

Jika saya mendekati ini dengan cara Euclidian, saya bisa mengukur jalur geodesik yang menghubungkan titik dasar, menemukan titik tengah sisi segitiga yang dihasilkan dan membuat ortodrom tegak lurus untuk masing-masing jalur tersebut. Tiga loxodromes mungkin akan bertemu pada titik yang sama. Jika ini adalah metode yang benar, pasti ada cara yang lebih mudah untuk melakukannya di Arc.

Jawaban ini dibagi menjadi beberapa bagian:

• Analisis dan Pengurangan Masalah , menunjukkan cara menemukan titik yang diinginkan dengan rutinitas "kalengan".

• Ilustrasi: Prototipe Bekerja , memberikan kode kerja.

• Contoh , menunjukkan contoh solusi.

• Perangkap , membahas potensi masalah dan cara mengatasinya.

• Implementasi ArcGIS , komentar tentang membuat alat ArcGIS khusus dan tempat untuk mendapatkan rutinitas yang diperlukan.

Analisis dan Pengurangan Masalah

Mari kita mulai dengan mengamati bahwa dalam model bola (bulat sempurna) akan selalu ada solusi - pada kenyataannya, tepat dua solusi. Diberikan titik dasar A, B, dan C, masing-masing pasangan menentukan "garis-beratnya tegak lurus," yang merupakan himpunan titik-titik yang berjarak sama dari dua titik yang diberikan. Garis bagi ini adalah geodesik (lingkaran besar). Geometri bola berbentuk bulat panjang : setiap dua geodesik berpotongan (dalam dua titik unik). Dengan demikian, titik-titik persimpangan garis-bagi AB dan garis-garis BC adalah - menurut definisi - berjarak sama dari A, B, dan C, dengan demikian menyelesaikan masalah. (Lihat gambar pertama di bawah ini.)

Hal-hal terlihat lebih rumit pada ellipsoid, tetapi karena itu adalah gangguan kecil dari bola, kita dapat mengharapkan perilaku yang sama. (Analisis ini akan membawa kita terlalu jauh.) Rumus rumit yang digunakan (secara internal dalam GIS) untuk menghitung jarak akurat pada ellipsoid bukanlah komplikasi konseptual, meskipun: masalahnya pada dasarnya sama. Untuk melihat betapa sederhananya masalahnya, mari kita nyatakan dengan agak abstrak. Dalam pernyataan ini, "d (U, V)" mengacu pada jarak benar dan akurat antara titik U dan V.

Diberikan tiga titik A, B, C (sebagai pasangan lat-lon) pada ellipsoid, cari titik X yang (1) d (X, A) = d (X, B) = d (X, C) dan ( 2) jarak bersama ini sekecil mungkin.

Ketiga jarak ini semua bergantung pada X yang tidak diketahui . Jadi perbedaan jarak u (X) = d (X, A) - d (X, B) dan v (X) = d (X, B) - d (X, C) adalah fungsi X yang bernilai riil. Sekali lagi, agak abstrak, kita dapat menggabungkan perbedaan-perbedaan ini menjadi pasangan yang teratur. Kami juga akan menggunakan (lat, lon) sebagai koordinat untuk X, memungkinkan kami untuk mempertimbangkannya sebagai pasangan yang dipesan, juga, katakanlah X = (phi, lambda). Dalam pengaturan ini, fungsinya

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

adalah fungsi dari bagian dari ruang dua dimensi yang mengambil nilai dalam ruang dua dimensi dan masalah kita berkurang menjadi

Temukan semua yang mungkin (phi, lambda) untuk yang F (phi, lambda) = (0,0).

Di sinilah abstraksi terbayar: ada banyak perangkat lunak yang hebat untuk memecahkan masalah ini (pencarian akar multidimensi numerik murni). Cara kerjanya adalah bahwa Anda menulis rutin untuk menghitung F , kemudian Anda meneruskannya ke perangkat lunak bersama dengan informasi tentang pembatasan inputnya ( phi harus berada di antara -90 dan 90 derajat dan lambda harus berada di antara -180 dan 180 derajat). Itu engkol pergi untuk sepersekian detik dan mengembalikan (biasanya) hanya satu nilai ( phi , lambda ), jika dapat menemukan satu.

Ada detail untuk ditangani, karena ada seni untuk ini: ada berbagai metode solusi untuk dipilih, tergantung pada bagaimana F "berperilaku"; ini membantu untuk "mengarahkan" perangkat lunak dengan memberikannya titik awal yang masuk akal untuk pencariannya (ini adalah salah satu cara kita bisa mendapatkan solusi terdekat , daripada yang lain); dan Anda biasanya perlu menentukan seberapa akurat solusi yang Anda inginkan (sehingga ia tahu kapan harus menghentikan pencarian). (Untuk informasi lebih lanjut tentang apa yang perlu diketahui analis GIS tentang perincian seperti itu, yang banyak muncul dalam masalah GIS, silakan kunjungi Topik yang disarankan untuk dimasukkan dalam kursus Ilmu Komputer untuk Teknologi Geospasial dan lihat di bagian "Lain-lain" di bagian akhir. )

Ilustrasi: Prototipe yang Berfungsi

Analisis menunjukkan kita perlu memprogram dua hal: perkiraan awal kasar solusi dan perhitungan F itu sendiri.

Perkiraan awal dapat dibuat dengan "rata-rata bulat" dari tiga titik dasar. Ini diperoleh dengan mewakili mereka dalam koordinat geosentris Cartesian (x, y, z), rata-rata koordinat tersebut, dan memproyeksikan rata-rata itu kembali ke bola dan mengekspresikannya kembali dalam lintang dan bujur. Ukuran bola tidak material dan perhitungannya dibuat langsung: karena ini hanyalah titik awal, kita tidak perlu perhitungan ellipsoidal.

Untuk prototipe yang berfungsi ini saya menggunakan 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
  ]

(Kondisi terakhir Ifmenguji apakah rata-rata mungkin gagal untuk menunjukkan dengan jelas garis bujur; jika demikian, ia jatuh kembali ke rata-rata aritmatika lurus dari lintang dan bujur inputnya - mungkin bukan pilihan yang bagus, tetapi setidaknya yang valid. Bagi mereka yang menggunakan kode ini untuk panduan implementasi, perhatikan bahwa argumen Mathematica ArcTan dibalik dibandingkan dengan kebanyakan implementasi lainnya: argumen pertama adalah koordinat x, yang kedua adalah koordinat y, dan mengembalikan sudut yang dibuat oleh vektor ( x, y).)

Sejauh bagian kedua, karena Mathematica - seperti ArcGIS dan hampir semua GIS lainnya - berisi kode untuk menghitung jarak akurat pada ellipsoid, hampir tidak ada yang bisa ditulis. Kami hanya memanggil rutinitas pencarian root:

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

Aspek yang paling penting dari implementasi ini adalah bagaimana ia menghindari kebutuhan untuk membatasi lintang ( f) dan bujur ( q) dengan selalu menghitungnya masing-masing modulo 180 dan 360 derajat. Ini menghindari keharusan membatasi masalah (yang sering menimbulkan komplikasi). Parameter kontrol MaxIterationsdll. Di-tweak untuk membuat kode ini memberikan akurasi setinggi mungkin.

Untuk melihatnya dalam tindakan, mari kita terapkan pada tiga poin dasar yang diberikan dalam pertanyaan terkait :

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

{-6.29692, 106.907}

Jarak yang dihitung antara solusi ini dan tiga titik adalah

{1450.23206979, 1450.23206979, 1450.23206978}

(Ini adalah meter). Mereka setuju melalui angka signifikan kesebelas (yang sebenarnya terlalu tepat, karena jarak jarang akurat hingga lebih baik dari satu milimeter atau lebih). Berikut adalah gambar dari tiga poin ini (hitam), tiga garis timbal balik mereka, dan solusinya (merah):

Contoh

Untuk menguji implementasi ini dan mendapatkan pemahaman yang lebih baik tentang bagaimana masalah itu terjadi, berikut adalah plot kontur dari perbedaan rata-rata kuadrat akar dalam jarak untuk tiga titik dasar yang tersebar luas. (Perbedaan RMS diperoleh dengan menghitung ketiga perbedaan d (X, A) -d (X, B), d (X, B) -d (X, C), dan d (X, C) -d (X , A), rata-rata kuadrat mereka, dan mengambil akar kuadrat.Ini sama dengan nol ketika X memecahkan masalah dan sebaliknya meningkat ketika X menjauh dari solusi, dan dengan demikian mengukur seberapa "dekat" kita untuk menjadi solusi di lokasi mana pun. )

Titik dasar (60, -120), (10, -40), dan (45,10) ditunjukkan dengan warna merah dalam proyeksi Plate Carree ini; solusinya (49.2644488, -49.9052992) - yang membutuhkan 0,03 detik untuk menghitung - berwarna kuning. Perbedaan RMS-nya kurang dari tiga nanometer , meskipun semua jarak yang relevan adalah ribuan kilometer. Area gelap menunjukkan nilai kecil RMS dan area terang menunjukkan nilai tinggi.

Peta ini jelas menunjukkan solusi lain yang terletak dekat (-49.2018206, 130.0297177) (dihitung dengan RMS dua nanometer dengan menetapkan nilai pencarian awal secara diametris berlawanan dengan solusi pertama.)

Perangkap

Ketidakstabilan angka

Ketika poin dasar hampir collinear dan berdekatan, semua solusi akan hampir setengah dunia dan sangat sulit untuk dijabarkan secara akurat. Alasannya adalah bahwa perubahan kecil di lokasi di seluruh dunia - bergerak menuju atau menjauh dari titik dasar - hanya akan menyebabkan perubahan kecil dalam perbedaan jarak. Hanya saja tidak ada cukup akurasi dan presisi yang dibangun dalam perhitungan jarak geodetik yang biasa untuk memberikan hasil.

Misalnya, dimulai dengan titik dasar di (45.001, 0), (45, 0), dan (44.999.0), yang dipisahkan sepanjang Prime Meridian dengan hanya 111 meter antara masing-masing pasangan, trimendapatkan solusi (11.8213, 77.745) ). Jarak dari itu ke titik dasar adalah 8.127.964.998 77; 8.127.964.998 41; dan 8.127.964.998 masing-masing 65 meter. Mereka setuju dengan milimeter terdekat! Saya tidak yakin seberapa akurat hasil ini, tetapi tidak akan sedikit terkejut jika implementasi lain mengembalikan lokasi yang jauh dari yang ini, menunjukkan kesetaraan hampir sama baiknya dari ketiga jarak.

Waktu perhitungan

Perhitungan ini, karena melibatkan pencarian yang cukup dengan menggunakan perhitungan jarak yang rumit, tidak cepat, biasanya membutuhkan sepersekian detik. Aplikasi waktu nyata perlu mengetahui hal ini.

Implementasi ArcGIS

Python adalah lingkungan skrip yang disukai untuk ArcGIS (dimulai dengan versi 9). The paket scipy.optimize memiliki rootfinder multivariat rootyang harus melakukan apa yang FindRootdilakukannya dalam Mathematica kode. Tentu saja ArcGIS sendiri menawarkan perhitungan jarak ellipsoidal yang akurat. Selebihnya, adalah semua detail implementasi: tentukan bagaimana koordinat titik dasar akan diperoleh (dari lapisan? Diketik oleh pengguna? Dari file teks? Dari mouse?) Dan bagaimana output akan disajikan (sebagai koordinat ditampilkan di layar "sebagai titik grafik" sebagai objek titik baru dalam sebuah lapisan?), tulis antarmuka itu, port kode Mathematica yang ditunjukkan di sini (langsung), dan Anda akan siap.

Seperti yang Anda perhatikan, masalah ini muncul dalam menentukan batas laut; sering disebut sebagai masalah "tri-point" dan Anda dapat Google ini dan menemukan beberapa makalah yang mengatasinya. Salah satu makalah ini adalah oleh saya (!) Dan saya menawarkan solusi konvergen yang akurat dan cepat. Lihat Bagian 14 dari http://arxiv.org/abs/1102.1215

Metode ini terdiri dari langkah-langkah berikut:

1. tebak tri-point O
2. gunakan O sebagai pusat proyeksi azimut yang berjarak sama
3. proyek A, B, C, untuk proyeksi ini
4. temukan tri-point dalam proyeksi ini, O '
5. gunakan O 'sebagai pusat proyeksi baru
6. ulangi sampai O 'dan O bertepatan

Formula yang diperlukan untuk solusi tri-point dalam proyeksi diberikan dalam makalah. Selama Anda menggunakan proyeksi azimuth akurat berjarak sama, jawabannya akan tepat. Konvergensi adalah makna kuadrat yang hanya diperlukan beberapa iterasi. Ini hampir pasti akan mengungguli metode pencarian root umum yang disarankan oleh @whuber.

Saya tidak bisa langsung membantu Anda dengan ArcGIS. Anda dapat mengambil paket python saya untuk melakukan perhitungan geodesik dari https://pypi.python.org/pypi/geographiclib dan mengkodekan proyeksi berdasarkan ini sederhana.

Edit

Masalah menemukan tri-point dalam kasus degenerate @ whuber (45 + eps, 0) (45.0) (45-eps, 0) dipertimbangkan oleh Cayley dalam On the geodesic lines pada sebuah spheroid oblate , Phil. Mag. (1870), http://books.google.com/books?id=4XGIOoCMYYAC&pg=PA15

Dalam hal ini, tri-point diperoleh dengan mengikuti geodesik dari (45,0) dengan azimuth 90 dan menemukan titik di mana skala geodesik menghilang. Untuk ellipsoid WGS84, titik ini adalah (-0.10690908732248, 89.89291072793167). Jarak dari titik ini ke masing-masing (45,001,0), (45,0), (44,999) adalah 10010287,665788943 m (dalam satu nanometer atau lebih). Ini sekitar 1882 km lebih dari perkiraan whuber (yang hanya menunjukkan betapa tidak stabilnya kasus ini). Untuk bumi berbentuk bola, tri-point akan menjadi (0,90) atau (0, -90), tentu saja.

TAMBAHKAN: Berikut adalah implementasi dari metode persamaan jarak azimut menggunakan 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

Menguji ini menggunakan Oktaf saya dapatkan

oktaf: 1> panjang format
oktaf: 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

sebagai titik tri untuk New York, Tokyo, dan Wellington.

Metode ini tidak akurat untuk titik kolinear tetangga, misalnya, [45,001,0], [45,0], [44,999,0]. Dalam hal ini, Anda harus memecahkan M ₁₂ = 0 pada berasal geodesik dari [45,0] di azimuth 90. Fungsi berikut melakukan trik (menggunakan metode 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

Sebagai contoh, ini memberi:

[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

Saya ingin tahu seberapa cepat pendekatan cffk menyatu pada solusi, jadi saya menulis tes menggunakan arcobjects, yang menghasilkan output ini. Jarak dalam meter:

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

Ini kode sumbernya. (Sunting) Mengubah FindCircleCenter untuk menangani persimpangan (titik tengah) yang jatuh dari tepi proyeksi azimut:

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

Ada juga pendekatan alternatif dalam edisi Juni 2013 dari Majalah MSDN, Amoeba Method Optimization menggunakan C # .

Edit

Kode yang sebelumnya diposting di konvergensi ke antipode dalam beberapa kasus. Saya telah mengubah kode sehingga menghasilkan output ini untuk titik uji @ cffk.

Inilah output yang sekarang dihasilkannya:

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

Berikut kode yang direvisi:

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

Edit

Inilah hasil yang saya dapatkan dengan 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