[GIS] Getting ST_ClosestPoint for Geography type

great circlepostgis

Using PostGIS 2.0.0 I'd like to find the point on a LINESTRING that is closest to a given point. The LINESTRING represents a great circle line (ie. geography type). ST_ClosestPoint appears to do exactly what I want, however, I find that the returned point does not lie on the great circle line:

WITH points AS (
    SELECT ST_GeomFromEWKT('SRID=4326;POINT(41 12)') AS point,
        ST_GeomFromEWKT('SRID=4326;LINESTRING(40 5,41 15)') AS border
)
SELECT ST_AsText(ST_ClosestPoint(border, point)), -- POINT(40.7029702970297 12.029702970297)
    ST_DWithin(ST_ClosestPoint(border, point)::geography, border::geography, 700) -- false
FROM points;

It's actually over 700 metres from the line. I believe this happens because ST_ClosestPoint operates on a plane. How would I do the same thing on the spheroid (or at least on the sphere)?

I tried @MerseyViking's function (ported to plpgsql) and it mostly works very well, but in some cases returns a point that is not on the line. A simple example:

WITH points AS (
    SELECT ST_GeomFromEWKT('SRID=4326;POINT(10 3)') AS point,
    ST_GeomFromEWKT('SRID=4326;LINESTRING(10 5,10 15)') AS border
),
calc AS
(
    SELECT ST_ClosestPoint(border, point)::geography AS plane_closest,
        closest_point_to_line(border, point) AS sphere_closest
    FROM points
)
SELECT ST_AsText(plane_closest) AS plane_closest, -- POINT(10 5)
    ST_DWithin(plane_closest, border::geography, 1), -- true
    ST_AsText(sphere_closest) AS sphere_closest, -- POINT(10 3)
    ST_DWithin(sphere_closest, border::geography, 1), -- false
FROM points, calc;

Best Answer

Sadly there isn't a geographic version of ST_ClosestPoint, so you will have to write your own function. There are two ways of calculating the nearest point of a great circle: spherical trigonometry, or 3D vector algebra. Luckily for you I have just written such a function for the latter method; I've not attempted the former because my spherical trig is pretty poor.

I've written a PL/Python function that you can add to your database, although you will need to install the Python libraries and PL/Python wrappers, which if you're using Linux should be simply to install the postgresql-python-x.x package (where x.x is the version of PostgreSQL you're using).

This function isn't necessarily the most efficient, it includes no bounds or error checking, and it only uses a sphere, but it certainly works well enough with the data I've tried it with. Feel free to modify accordingly.

First I define my own simple long/lat type:

CREATE TYPE lon_lat AS (lon float, lat float);

The function looks like this:

CREATE OR REPLACE FUNCTION geog_nearest_point(lp1 lon_lat, lp2 lon_lat, p lon_lat) RETURNS lon_lat AS
$$
import math
def sph2cart(lon, lat, r):
    return (r * math.cos(lat) * math.cos(lon),
            r * math.cos(lat) * math.sin(lon),
            r * math.sin(lat))

def cart2sph(x, y, z):
    return (math.atan2(y,x),
            math.atan2(z, math.sqrt(x * x + y * y)),
            math.sqrt(x * x + y * y + z * z))

def cross_prod(v1, v2):
    return (v1[1] * v2[2] - v1[2] * v2[1],
            v1[2] * v2[0] - v1[0] * v2[2],
            v1[0] * v2[1] - v1[1] * v2[0])

lv1 = sph2cart(math.radians(lp1['lon']), math.radians(lp1['lat']), 1.0)
lv2 = sph2cart(math.radians(lp2['lon']), math.radians(lp2['lat']), 1.0)
pv = sph2cart(math.radians(p['lon']), math.radians(p['lat']), 1.0)

f = cross_prod(lv1, lv2)
g = cross_prod(pv, f)
h = cross_prod(f, g)

nearest_point = cart2sph(h[0], h[1], h[2])
return (math.degrees(nearest_point[0]), math.degrees(nearest_point[1]))
$$
LANGUAGE 'plpythonu' IMMUTABLE;

And giving it the data you provided gives me this:

geog_db=# SELECT geog_nearest_point((41, 15)::lon_lat, (40, 5)::lon_lat, (41, 12)::lon_lat);
      geog_nearest_point
-------------------------------
 (40.6960040707,12.0295700818)
(1 row)
  • It works by finding the point f on the sphere which is normal to the plane that describes the great circle. This point is unique for each great circle.
  • Next, it uses f and the input point p to define a plane perpendicular to the input great circle, and determines its unique normal point g. The vectors from the centre of the sphere to f and g define a plane that is perpendicular to the previous two planes.
  • Finding the normal of this third plane gives us the point at which the two great circles intersect.
  • Finally, this Cartesian coordinate is converted back to latitude and longitude.

If anyone has a spherical trig method, which I think will be slightly more efficient, then I'd love to see it.

Related Question