[Math] Closest grid square to a point in spherical coordinates

algorithmsmg.metric-geometry

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: https://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.

Best Answer

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 θ), π)

Related Question