Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

For example, I have coordinates for three base points on a coastline and I need to find the coordinates of the point off the coast which is equidistant from all three. It’s a simple exercise in geometry, but all the measurements must take geodesy into account.

If I was approaching this in a Euclidian manner, I could measure the geodesic paths connecting the base points, find the midpoints of the sides of the resulting triangle and create perpendicular loxodromes to each of those paths. The three loxodromes would presumably converge at the equidistant point. If this is the correct method, there’s got to be an easier way to do it in Arc.

Many thanks!

I need to find O

share|improve this question
    
Are there constraints on the relative positions of the 3 points? Picture east coast, middle point is farthest east. Your solution wouldn't work as the perpendiculars wouldn't converge offshore. I'm sure we can come up with other bad cases! –  mkennedy May 30 '13 at 23:46
    
I wonder if you could use a distance preserving projection and run the computation from there? progonos.com/furuti/MapProj/Normal/CartProp/DistPres/… Not sure of the algorithm to do it, there must be one... maybe it's the barycentre: en.wikipedia.org/wiki/Barycentric_coordinate_system –  alexgleith May 31 '13 at 4:45
    
For solutions to a closely related problem, search our site for "trilateration". Also, gis.stackexchange.com/questions/10332/… is a duplicate but has no adequate answers (most likely because the question was asked in a confused manner). –  whuber May 31 '13 at 8:27
    
@mkennedy There are, in principle, no bad cases, only numerically unstable ones. These occur when the three base points are collinear; the two solutions (on a spherical model) occur at the two poles of the common geodesic; in an ellipsoidal model they occur near where those poles would be expected. –  whuber May 31 '13 at 15:10
    
The use of loxodromes here would not be correct: they are not the perpendicular bisectors. On the sphere, these lines will be parts of great circles (geodesics), but on the ellipsoid, they will depart slightly from being geodesics. –  whuber May 31 '13 at 16:55
add comment

3 Answers

This answer is divided into multiple sections:

  • Analysis and Reduction of the Problem, showing how to find the desired point with "canned" routines.

  • Illustration: a Working Prototype, giving working code.

  • Example, showing examples of the solutions.

  • Pitfalls, discussing potential problems and how to cope with them.

  • ArcGIS implementation, comments about creating a custom ArcGIS tool and where to obtain the needed routines.


Analysis and Reduction of the Problem

Let's begin by observing that in the (perfectly round) spherical model there will always be a solution--in fact, exactly two solutions. Given base points A, B, and C, each pair determines its "perpendicular bisector," which is the set of points equidistant from the two given points. This bisector is a geodesic (great circle). Spherical geometry is elliptic: any two geodesics intersect (in two unique points). Thus, the intersection points of AB's bisector and BC's bisector are--by definition--equidistant from A, B, and C, thereby solving the problem. (See the first figure below.)

Things look more complicated on an ellipsoid, but because it is a small perturbation of the sphere, we can expect similar behavior. (The analysis of this would take us too far afield.) The complicated formulas used (internally within a GIS) to compute accurate distances on an ellipsoid are not a conceptual complication, though: the problem is basically the same. To see how simple the problem really is, let's state it somewhat abstractly. In this statement, "d(U,V)" refers to the true, fully accurate distance between points U and V.

Given three points A, B, C (as lat-lon pairs) on an ellipsoid, find a point X for which (1) d(X,A) = d(X,B) = d(X,C) and (2) this common distance is as small as possible.

These three distances all depend on the unknown X. Thus the differences in distances u(X) = d(X,A) - d(X,B) and v(X) = d(X,B) - d(X,C) are Real-valued functions of X. Again, somewhat abstractly, we may assemble these differences into an ordered pair. We will also use (lat, lon) as coordinates for X, allowing us to consider it as an ordered pair, too, say X = (phi, lambda). In this setup, the function

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

is a function from a portion of a two-dimensional space taking values in two-dimensional space and our problem reduces to

Find all possible (phi, lambda) for which F(phi, lambda) = (0,0).

Here is where the abstraction pays off: plenty of great software exists to solve this (purely numerical multidimensional root-finding) problem. The way it works is that you write a routine to compute F, then you pass it along to the software along with any information about restrictions on its input (phi must lie between -90 and 90 degrees and lambda must lie between -180 and 180 degrees). It cranks away for a fraction of a second and returns (typically) just one value of (phi, lambda), if it can find one.

There are details to handle, because there's an art to this: there are various solution methods to choose from, depending on how F "behaves"; it helps to "steer" the software by giving it a reasonable starting point for its search (this is one way we can obtain the nearest solution, rather than any other one); and you usually need to specify how accurate you would like the solution to be (so it knows when to stop the search). (For more on what GIS analysts need to know about such details, which come up a lot in GIS problems, please visit Recommend topics to be included in a Computer Science for Geospatial Technologies course and look in the "Miscellany" section near the end.)


Illustration: a Working Prototype

The analysis shows we need to program two things: a crude initial estimate of the solution and the calculation of F itself.

The initial estimate can be made by a "spherical average" of the three base points. This is obtained by representing them in geocentric Cartesian (x,y,z) coordinates, averaging those coordinates, and projecting that average back to the sphere and re-expressing it in latitude and longitude. The size of the sphere is immaterial and the calculations are thereby made straightforward: because this is just a starting point, we don't need ellipsoidal calculations.

For this working prototype I used 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
  ]

(The final If condition tests whether the average might fail clearly to indicate a longitude; if so, it falls back to a straight arithmetic mean of the latitudes and longitudes of its input--maybe not a great choice, but at least a valid one. For those using this code for implementation guidance, note that the arguments of Mathematica's ArcTan are reversed compared to most other implementations: its first argument is the x-coordinate, its second is the y-coordinate, and it returns the angle made by the vector (x,y).)

As far as the second part goes, because Mathematica--like ArcGIS and almost all other GISes--contains code to compute accurate distances on the ellipsoid, there's almost nothing to write. We just call the root-finding routine:

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

The most noteworthy aspect of this implementation is how it dodges the need to constrain the latitude (f) and longitude (q) by always computing them modulo 180 and 360 degrees, respectively. This avoids having to constrain the problem (which often creates complications). The control parameters MaxIterations etc. are tweaked to make this code provide the greatest possible accuracy it can.

To see it in action, let's apply it to the three base points given in a related question:

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

{-6.29692, 106.907}

The computed distances between this solution and the three points are

{1450.23206979, 1450.23206979, 1450.23206978}

(these are meters). They agree through the eleventh significant digit (which is too precise, actually, since distances are rarely accurate to better than a millimeter or so). Here's a picture of these three points (black), their three mutual bisectors, and the solution (red):

Figure 1


Example

To test this implementation and get a better understanding of how the problem behaves, here is a contour plot of the root mean square discrepancy in distances for three widely spaced base points. (The RMS discrepancy is obtained by computing all three differences d(X,A)-d(X,B), d(X,B)-d(X,C), and d(X,C)-d(X,A), averaging their squares, and taking the square root. It equals zero when X solves the problem and otherwise increases as X moves away from a solution, and so measures how "close" we are to being a solution at any location.)

Figure 2

The base points (60,-120), (10,-40), and (45,10) are shown in red in this Plate Carree projection; the solution (49.2644488, -49.9052992)--which required 0.03 seconds to compute--is in yellow. Its RMS discrepancy is less than three nanometers, despite all relevant distances being thousands of kilometers. The dark areas show small values of the RMS and the light areas show high values.

This map clearly shows another solution lies near (-49.2018206, 130.0297177) (computed to an RMS of two nanometers by setting the initial search value diametrically opposite the first solution.)


Pitfalls

Numerical instability

When the base points are nearly collinear and close together, all solutions will be nearly half a world away and extremely difficult to pin down accurately. The reason is that small changes in a location across the globe--moving it towards or away from the base points--will induce only incredibly tiny changes in the differences of distances. There's just not enough accuracy and precision built into the usual computation of geodetic distances to pin down the result.

For example, starting with the base points at (45.001, 0), (45, 0), and (44.999,0), which are separated along the Prime Meridian by only 111 meters between each pair, tri obtains the solution (11.8213, 77.745). The distances from it to the base points are 8,127,964.998 77; 8,127,964.998 41; and 8,127,964.998 65 meters, respectively. They agree to the nearest millimeter! I'm not sure how accurate this result may be, but would not be in the least surprised if other implementations returned locations far from this one, showing almost as good equality of the three distances.

Computation time

These calculations, because they involve considerable searching using complicated distance computations, are not fast, usually requiring a noticeable fraction of a second. Real-time applications need to be aware of this.


ArcGIS implementation

Python is the preferred scripting environment for ArcGIS (beginning with version 9). The scipy.optimize package has a multivariate rootfinder root which should do what FindRoot does in the Mathematica code. Of course ArcGIS itself offers accurate ellipsoidal distance calculations. The rest, then, is all implementation details: decide how the base point coordinates will be obtained (from a layer? typed in by the user? from a text file? from the mouse?) and how the output will be presented (as coordinates displayed on screen? as a graphic point? as a new point object in a layer?), write that interface, port the Mathematica code shown here (straightforward), and you will be all set.

share|improve this answer
3  
+1 Very thorough. I think you might have to start charging for this, @whuber. –  Nick Ochoski May 31 '13 at 16:56
2  
@Radar Thanks. I'm hoping people will buy my book when(ever) it eventually appears :-). –  whuber May 31 '13 at 16:57
1  
Will do Bill ... Send a draft!!! –  Dan Patterson May 31 '13 at 22:48
    
Excellent! Still, it seems like an analytical solution would be possible. By restating the problem into 3d cartesian space with 3 points A,B,C and E where E is the center of the earth. Next find two planes Plane1 and Plane2. Plane1 would be the plane that is normal to planeABE and passing through E,midpoint(A,B). Likewise, Plane2 would be the plane normal to planeACE and passing through E,midpoint(C,E). The lineO formed by intersection of Plane1 and Plane2 represents points equidistant to the 3 points. The closer of the two points to A(or B or C) where lineO intersects sphere is pointO. –  Kirk Kuykendall Jun 1 '13 at 2:20
    
That analytical solution, @Kirk, applies only to the sphere. (Intersections of planes with the ellipsoid are never perpendicular bisectors in the metric of the ellipsoid except for a few obvious exceptional cases: when they are meridians or the Equator.) –  whuber Jun 3 '13 at 12:40
show 2 more comments

As you note, this problem arises in determining maritime boundaries; it's often referred to as the "tri-point" problem and you can Google this and find several papers addressing it. One of these papers is by me(!) and I offer an accurate and rapidly convergent solution. See Section 14 of http://arxiv.org/abs/1102.1215

The method consists of the following steps:

  1. guess a tri-point O
  2. use O as the center of an azimuthal equidistant projection
  3. project A, B, C, to this projection
  4. find the tri-point in this projection, O'
  5. use O' as the new center of projection
  6. repeat until O' and O coincide

The necessary formula for the tri-point solution in the projection is given in the paper. As long as you're using an accurate azimuthal equidistant projection, the answer will be exact. Convergence is quadratic meaning that only a few iterations are needed. This will almost certainly outperform the general root-finding methods suggested by @whuber.

I can't directly help you with ArcGIS. You can grab my python package for doing geodesic calculations from https://pypi.python.org/pypi/geographiclib and coding up the projection based on this is simple.


Edit

The problem of finding the tri-point in @whuber's degenerate case (45+eps,0) (45,0) (45-eps,0) was considered by Cayley in On the geodesic lines on an oblate spheroid, Phil. Mag. (1870), http://books.google.com/books?id=4XGIOoCMYYAC&pg=PA15

In this case, the tri-point is obtained by following a geodesic from (45,0) with azimuth 90 and finding the point at which the geodesic scale vanishes. For the WGS84 ellipsoid, this point is (-0.10690908732248, 89.89291072793167). The distance from this point to each of (45.001,0), (45,0), (44.999) is 10010287.665788943 m (within a nanometer or so). This is about 1882 km more than whuber's estimate (which just goes to show how unstable this case is). For a spherical earth, the tri-point would be (0,90) or (0,-90), of course.

ADDENDUM: Here's an implementation of the azimuthal equidistant method using 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

Testing this out using Octave I get

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

as the tri point for New York, Tokyo, and Wellington.

This method is inaccurate for neighboring colinear points, e.g., [45.001,0], [45,0], [44.999,0]. In that case, you should solve for M12 = 0 on a geodesic emanating from [45,0] at azimuth 90. The following function does the trick (using Newton's method):

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

For the example, this gives:

[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
share|improve this answer
    
+1. However, it's unclear that a general root finder will perform any less well: the function is nicely behaved near its optimum and Newton's method, for example, will also converge quadratically. (Mathematica is typically taking about four steps to converge; each step requires four evaluations to estimate the Jacobian.) The real advantage I see to your method is that it can readily be scripted in a GIS without resorting to use of a root finder. –  whuber Jun 5 '13 at 21:13
    
I agree. My method is equivalent to Newton and so, in contrast to Mathematica's root finding method, there's no need to estimate the gradient by taking differences. –  cffk Jun 5 '13 at 21:32
    
Right--but you have to do the reprojection each time, which looks like it's about the same amount of work. I do appreciate the simplicity and elegance of your approach, though: it's immediately obvious that it should work and will converge quickly. –  whuber Jun 5 '13 at 22:24
    
I've posted results for the same test points in my answer. –  Kirk Kuykendall Jun 11 '13 at 14:02
add comment

I was curious to see how quickly @cffk's approach converges on a solution, so I wrote a test using arcobjects, which produced this output. Distances are in meters:

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

Here is the source code. (Edit) Changed the FindCircleCenter to handle intersections (centerpoints) that fall off the edge of the azimuthal projection:

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

There's also an alternative approach in June 2013 issue of MSDN Magazine, Amoeba Method Optimization using C#.


Edit

The previously posted code converged to the antipode in some cases. I've altered the code so that it produces this output for @cffk's test points.

Here's the output it now produces:

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

Here's the revised code:

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

Here's the results I get with 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
share|improve this answer
    
That's impressively quick convergence! (+1) –  whuber Jun 8 '13 at 19:03
    
You should be using an honest-to-goodness azimuthal equidistant projection centered on newCenter. Instead you're using the projection centered on the N pole and shifting the origin to newCenter. It may therefore be accidental that you get a decent solution in this case (perhaps because the points are close together?). It would be good to try it with 3 points thousands of km apart. An implementation of the azimuthal equidistant projection is given in mathworks.com/matlabcentral/fileexchange/… –  cffk Jun 11 '13 at 2:29
    
@cffk The only other way I see to create an azimuthal equidistant projection centered on a particular point is to use the same method but with esriSRProjCS_World_AzimuthalEquidistant instead of esriSRProjCS_WGS1984N_PoleAziEqui (or esriSRProjCS_WGS1984S_PoleAziEqui). The only difference though, is that it's centered on 0,0 instead of 0,90 (or 0,-90). Can you guide me in running a test with the mathworks to see if this produces different results from an "honest-to-goodness" projection? –  Kirk Kuykendall Jun 11 '13 at 3:12
    
@KirkKuykendall: see the addendum to my first answer. –  cffk Jun 11 '13 at 4:24
1  
@KirkKuykendall So maybe ESRI did implement an "honest-to-goodness" projection? The key property required for this algorithm to work is that distances measured from the "center point" are true. And this property is easy enough to check. (The azimuthal property relative to the center point is secondary and may only affect the rate of convergence.) –  cffk Jun 12 '13 at 1:09
show 3 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.