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

Question

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.

Answer 1

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
radian = 1 / degree
radius <- 6371007.2 # meters

# 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.

Answer 2

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])});
    A[0] *= radius;
    A[1] *= radius;
    B[0] *= radius;
    B[1] *= radius;
    C[0] *= radius;
    C[1] *= radius;

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

}