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.

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

Please, if any one knows a good solution help.

Thanks in response :)

Imagine this on a globe or sphere

share|improve this question
    
Do you have a platform in mind? C library, Python module, spatial database? –  DavidF Oct 16 '12 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%. –  whuber Oct 16 '12 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. –  Batjaa Oct 17 '12 at 6:33
    
Related question gis.stackexchange.com/questions/14909/… –  Kirk Kuykendall Oct 17 '12 at 15:16
add comment

1 Answer

up vote 2 down vote accepted

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.

share|improve this answer
1  
thanks whuber! just implemented it in java and working fine. I truly appreciate your help :) –  Batjaa Oct 18 '12 at 8:27
    
@Batjaa Could i see your java implementation? –  sabotero Apr 9 at 18:06
add comment

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.