# Closest grid square to a point in spherical coordinates

I am programming an algorithm where I have broken up the surface of a sphere into grid points (for simplicity I have the grid lines are parallel and perpendicular to the meridians). Given a point A on the sphere, I would like to be able to efficiently take any grid "square" and determine the point B in the square with the least spherical coordinate distance AB.

• Near the poles, the "squares" degenerate into "triangles".
• As noted by Trithematician, this is a special case of to finding the shortest distance from a point to an arc on a sphere. For the longitudinal lines, these are arcs of great circles, but for the latitudinal they are not.
• TonyK provides a method below that solves the longitudinal case.
• Minimising the distance in 3d also minimises the distance on the sphere

I actually asked this question here: http://stackoverflow.com/questions/1463606/closest-grid-square-to-a-point-in-spherical-coordinates and decided to use an approximation instead, but I am still curious as to whether there is an exact solution. I thought that since this is really more of a maths question than a programming question that I would repost it here.

-
A couple of questions to clarify: 1. Is the point A on the sphere or not? 3. What do you mean by "spherical coordinate distance"? Is it the Euclidean distance when the points are expressed in spherical coordinates, or something else? –  Elisha Peterson Oct 17 '09 at 11:55
A is on the sphere - now clarified in the question. The spherical coordinate distance is the distance along the spheres surface. However, this is largely irrelevant as minimising the distance along the sphere also minimises the distance in Euclidean space –  Casebash Oct 18 '09 at 12:52
Sorry that I wasn't the best question asker. I should have posted the progress I did make on this problem as well as the question. –  Casebash Oct 23 '09 at 3:33

The tricky part is to find the nearest point on a meridian to a given point P. Let's fix the meridian M at φ=0, to keep the algebra simple.
Suppose P has spherical coordinates (θ, φ), with 0 ≤ θ ≤ π and -π < φ ≤ π. Let C be the great circle through P which is perpendicular to M; then we are looking for an intersection of C and M. (There are two such intersections; in the end we will simply choose the one that is on the same side of the equator as P).
The normal of C meets the sphere on the meridian M, at A = (ρ, 0), say. P is on C, so OP is perpendicular to OA, where O is the centre of the sphere. In cartesian coordinates, with r = 1 for simplicity:

A = (sin ρ, 0, cos ρ)
P = (cos φ sin θ, sin φ sin θ, cos θ)

The scalar product is

A.P = sin ρ cos φ sin θ + cos ρ cos θ

This must be zero, giving

tan ρ = -1/(cos φ tan θ)

The required point, Q say, is on M and perpendicular to A, so if Q = (σ, 0) in spherical coordinates, we have

tan σ = -1/tan ρ = cos φ tan θ

So

Q = (arctan (cos φ tan θ), 0)

If this Q is on the wrong side of the equator, take

Q = (π - arctan (cos φ tan θ), π)

-
Again, the intersection points may not be on the arc –  Casebash Oct 20 '09 at 7:42
Well, then, that's easy - you just take the corner of the square that's nearest to P. –  TonyK Oct 20 '09 at 7:54
That solves the meridians (longitudinal lines), but what about the lines of latitude? –  Casebash Oct 20 '09 at 20:25
"The normal of C". C is a circle? Do you mean the normal through the center of C? –  Casebash Oct 20 '09 at 20:44
Yep, that works given the above understanding of a normal for a circle, that works for the longitudinals –  Casebash Oct 20 '09 at 21:01

So the question essentially is: given a point A on the sphere, and a curve, what is the closest point on the curve. Standard optimization guarantees this occurs at a point D where the curve is perpendicular to the geodesic AD. (And geodesics are pieces of great circles.)

In your case, the curves are limited to theta=constant and phi=constant (parts of the spherical rectangles).

If phi=constant, the perpendicular geodesics are all great circles through the north pole, i.e. theta=constant, so the closest point has the same value of theta as A.

If the curve is theta=constant, there are a couple of ways to go. You might try Lagrange multipliers or some spherical trigonometry.

-
Won't a great circle through A be perpendicular to exactly two points on the extended arc (unless the curve is equatorial). If these points are not included in the arc, then the optimisation fails –  Casebash Oct 20 '09 at 7:40
Actually, then you just take the closest point on the circle (prove using the spherical law of sin. –  Casebash Oct 22 '09 at 2:14

Here is a meta-answer, without going into detail on the underlying optimization problems.

From the title of your question, and from your description, it sounds like you have two different problems:

Start with a sphere subdivided into rectilinear regions in spherical coordinates, for arguments sake, let's say each region is closed, and an arbitrary point A.

Problem #1. Determine which region(s) contains point A. The point A can only be contained in more than one region if it lies on the boundary of each region. Let's assume you have already solved this one. Problem #2. Given an arbitrary region R, determine B, the closest point of R to A (in spherical distance or Euclidean distance, doesn't matter which). Solution: Two cases: Case (1) A is in R, in which case B=A. Case (2) A is not in R, in which case B lies on the boundary of R. The boundary of R in general [2.0] consists of 4 arcs: two meridians of longitude, which are arcs of great circles and two parallels of latitude, which are arcs of small circles. There are also special and degenerate subcases: [2.1] R is a sector of a north or south polar cap, in which case the boundary consists only of 3 arcs: two meridians meeting in a pole, and one parallel; [2.2] R is a whole polar cap, with only one parallel as a boundary; [2.3] R is the difference between two caps, with a boundary consisting of two parallels; [2.4] R is a sector of the whole sphere from one pole to the other, with two meridians as a boundary; [2.5] R is a hemisphere with both poles on the boundary, where the two meridians degenerate into a single great circle; [2.6] R is the whole sphere - this is excluded anyway, since in this case, A lies in R.

In the subcase [2.0], B may lie on any of the 4 arcs or on the 4 vertices which are the intersections of pairs of arcs. In any case, the small circle with centre A and radius AB either intersects R only at B, in which case B is unique, or all along one of the arcs, or at two vertices, in which cases B is not unique. Intersection along all of a parallel arc occurs only when A is one of the poles. Intersection along all of a meridian arc occurs only when A is on the equator at a longitude 90 degrees away from the meridian arc. Intersection at two vertices can occur when the longitude of A is halfway between the two meridian arcs (or 180 degrees away from halfway in between, far enough away that the radius AB is greater than the radius of the smaller of the two radii of the parallel arcs, but close enough that the larger radius parallel arc is not closer to A than the smaller radius arc.

As per the other answers, you can treat [2.0] as 4 constrained optimization problems: find the closest point to A on each of the arcs. Once you have all 4 solutions, pick the closest one.

(Still in [2.0]) To get a (very rough) approximation, let's assume you can easily find the "centre" C of the region R in spherical coordinates. Now try to find the intersection (if any) between the arc AC and each of the boundary arcs of R. This will result in a point D which is either on one of the boundary arcs, or is one of the 4 vertices.

-

I added another comment to my own answer, but it was the sixth, so it has been hidden by the display software; perhaps you didn't notice it, Casebash:

And the latitudinals (did I get that right?) are trivial, as Triathematician explained.

Paraphrasing Triathematician's answer: the nearest point from P = (θ, φ) to the line of latitude given by θ = θ0 is simply the point (θ0, φ).

-
Well, in the case where that point isn't inside the segment, we have to just take the closest point instead. I only just figured out that it can be proved with the sin law (in case you can't tell I am totally new to spherical geometry). –  Casebash Oct 22 '09 at 2:22