# Line intersection with circle on a sphere (globe or earth)

How can I find the intersection point of a line and a circle on a sphere (earth surface)

Please, if any one knows a good solution help.

• Do you have a platform in mind? C library, Python module, spatial database? Commented Oct 16, 2012 at 15:57
• How much accuracy is needed? Spherical trig will make short (and easy) work of this (using, for instance, the spherical Law of Cosines applied to triangle ACD), but it may have errors up to about 0.3%. Commented Oct 16, 2012 at 16:14
• I'm planning to implement it in java, the accuracy should be high. because mostly of the circles radius will be about 500 meters and some to few kilometers. I'm given the latitude and longitude of line start and destination point, and of course the radius of the circle in meters, and center againg in lat, lon. Commented Oct 17, 2012 at 6:33
• Related question gis.stackexchange.com/questions/14909/… Commented Oct 17, 2012 at 15:16
• Note that you're actually trying to find the intersection of two circular arcs on a sphere, rather than a circle and a line. (with one arc spanning 2π and the other spanning a much smaller angle), widening the scope of what you can search for. Commented Jan 27 at 17:11

For such short distances any reasonable low-distortion projection will work. Perhaps the simplest is to project (lat, lon) to (R*lat, R*cos(lat0) * lon) where `R` is the authalic earth radius of 6371007.2 meters and `lat0` is a typical latitude for the problem, such as the latitude of the circle's center. Solve the problem in this projection.

### Example

The input data are

``````c, a, b = (48.137024, 11.575249), (48.139115, 11.578081), (48.146303, 11.593102)
r       = 1000
``````

The locations are expressed as (lat, lon) in degrees and the distance is in meters. We can find the point(s) of intersection in the form

``````x = a + t(b-a)
``````

for some real-valued parameter t. Let X be the projection of x and C the projection of c. (Notice that, since the projection is linear, X = A + t(B-A) where A is the projection of a and B is the projection of b.) The defining equation is

``````r = Distance(X, C) = Norm(X - C) = Norm(A-C + t(B-A)).
``````

Squaring both sides and expanding,

``````r^2 = Norm(A-C + t(B-A))^2 = (A-C).(A-C) + 2(A-C).(B-A) t + (B-A).(B-A) t^2.
``````

(where `.`, as in `(A-C).(B-A)`, denotes the dot product). This quadratic equation is readily solved for up to two values of t which, when plugged in to the defining equation for x, give up to two solutions. If an intersection along arc ab is sought, then only values of t between 0 and 1 will be valid.

With the sample data, the solution is (48.142739, 11.585654). Using calculations for the ITRF00 ellipsoid, this point is actually 1001.74 meters from the center--within the 0.3% expected for the difference between a spherical calculation (on which the projection is based) and an ellipsoidal calculation.

Here is (working, tested) `R` code illustrating the full calculation.

``````# Constants
degree <- 2 * pi / 360

# Input
a0 <- c(48.139115, 11.578081) * degree # Point A in (lat, lon)
b0 <- c(48.146303, 11.593102) * degree # Point B
c0 <- c(48.137024, 11.575249) * degree # Center
r <- 1000.0                            # Radius (meters)

# Projection
A <- a0 * c(1, cos(c0[1])) * radius
B <- b0 * c(1, cos(c0[1])) * radius
C <- c0 * c(1, cos(c0[1])) * radius

# Compute coefficients of the quadratic equation
v <- A - C; u <- B - A
alpha <- sum(u * u)
beta <- sum(u * v)
gamma <- sum(v * v) - r^2

# Solve the equation.
`%pm%` <- function(x,y) c(x+y, x-y)
t <- (-beta %pm% sqrt(beta^2 - alpha * gamma)) / alpha
t <- t[0 <= t & t <= 1]             # Limit the solution to arc a0-b0.
x <- (a0 + (b0-a0) %o% t) * radian  # Columns are (lat, lon)
``````

The quadratic formula used for the solution might need to be replaced by something better in order to handle cases where `alpha` is close to zero (that is, when points `a` and `b` are extremely close). It also should be enhanced to detect the case `beta^2 - alpha * gamma < 0`, when no solutions exist. Actual details depend on the programming environment: many of them already supply good equation solvers.

If greater than 0.3% accuracy is needed, then a projection based on an ellipsoidal datum is required.

• thanks whuber! just implemented it in java and working fine. I truly appreciate your help :) Commented Oct 18, 2012 at 8:27
• @Batjaa please post or send your java implementation somewhere. Commented Apr 22, 2014 at 16:40
• If this becomes solving the quadratic equation, shouldn't that be `-β/2α +/- sqrt(β² - 4αγ)/2α`? (as written, it seems to solve `-β/α +/- sqrt(β² - αγ)/α` instead?) Commented Jan 27 at 0:51
• @Mike `beta` in the code is half of your beta. Commented Jan 27 at 18:54
• ahh, thank you. Commented Jan 27 at 19:11

Java implementation of solution above, in case somebody ever need this)

``````public class GISIntersectionPointsCalculator implements Serializable {
static final double degree = 2 * Math.PI / 360;
static final double radian = 1 / degree;
static final double radius = 6371007.2;

static double[][] calculateIntersectionPoints(PointLatLon a, PointLatLon b, PointLatLon center, double r){

double[] a0 = {a.getLatitude() * degree, a.getLongitude() * degree}; // Point A in (lat, lon)
double[] b0 = {b.getLatitude() * degree, b.getLongitude() * degree}; // Point B
double[] c0 = {center.getLatitude() * degree, center.getLongitude() * degree}; // Center

double[] A = multiplyArray(a0, new double[]{1, Math.cos(c0[0])});
double[] B = multiplyArray(b0, new double[]{1, Math.cos(c0[0])});
double[] C = multiplyArray(c0, new double[]{1, Math.cos(c0[0])});

// Compute coefficients of the quadratic equation
double[] v = subtractArrays(A, C);
double[] u = subtractArrays(B, A);
double alpha = dotProduct(u, u);
double beta = dotProduct(u, v);
double gamma = dotProduct(v, v) - Math.pow(r, 2);

// Solve the equation
double[] t = solveQuadraticEquation(alpha, beta, gamma);
t = filterSolution(t); // Limit the solution to arc a0-b0
double[][] x = calculateCoordinates(a0, b0, t); // Columns are (lat, lon)

return x;
}

private static double[] multiplyArray(double[] arr, double[] multiplier) {
return new double[]{arr[0] * multiplier[0], arr[1] * multiplier[1]};
}

private static double[] subtractArrays(double[] arr1, double[] arr2) {
return new double[]{arr1[0] - arr2[0], arr1[1] - arr2[1]};
}

private static double dotProduct(double[] arr1, double[] arr2) {
return arr1[0] * arr2[0] + arr1[1] * arr2[1];
}

private static double[] filterSolution(double[] t) {
return Arrays.stream(t).filter(value -> value >= 0 && value <= 1).toArray();
}

private static double[] solveQuadraticEquation(double alpha, double beta, double gamma) {
double discriminant = Math.pow(beta, 2) - alpha * gamma;
double[] solutions = new double[]{
(-beta + Math.sqrt(discriminant)) / alpha,
(-beta - Math.sqrt(discriminant)) / alpha
};
return solutions;
}

private static double[][] calculateCoordinates(double[] a0, double[] b0, double[] t) {
double[][] coordinates = new double[t.length][2];
for (int i = 0; i < t.length; i++) {
coordinates[i][0] = (a0[0] + (b0[0] - a0[0]) * t[i]) * radian;
coordinates[i][1] = (a0[1] + (b0[1] - a0[1]) * t[i]) * radian;
}
return coordinates;
}
``````

}