Getting ST_ClosestPoint for Geography type?

Question

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;

Answer 1

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.

Answer 2

This is the PL/pgSQL version of MerseyViking's code. It also uses PostGIS geography Point and LineString types rather than a custom type to represent coordinates.

CREATE OR REPLACE FUNCTION _point_to_cartesian(point geometry(Point), radius float, OUT x float, OUT y float, OUT z float)
RETURNS RECORD AS $BODY$
DECLARE
    lon float;
    lat float;
BEGIN
    lon := radians(ST_X(point));
    lat := radians(ST_Y(point));

    x := radius * cos(lat) * cos(lon);
    y := radius * cos(lat) * sin(lon);
    z := radius * sin(lat);
END
$BODY$ LANGUAGE plpgsql IMMUTABLE;


CREATE OR REPLACE FUNCTION _cartesian_to_point(x float, y float, z float)
RETURNS geometry(Point) AS $BODY$
DECLARE
    lon float;
    lat float;
BEGIN
    lon := degrees(atan2(y, x));
    lat := degrees(atan2(z, sqrt(x * x + y * y)));
    -- Z coordinate ignored: degrees(sqrt(x * x + y * y + z * z));

    RETURN ST_MakePoint(lon, lat);
END
$BODY$ LANGUAGE plpgsql IMMUTABLE;


CREATE OR REPLACE FUNCTION _cross_product(x1 float, y1 float, z1 float, x2 float, y2 float, z2 float,
    OUT x float, OUT y float, OUT z float)
RETURNS RECORD AS $BODY$
BEGIN
    x := y1 * z2 - z1 * y2;
    y := z1 * x2 - x1 * z2;
    z := x1 * y2 - y1 * x2;
END
$BODY$ LANGUAGE plpgsql IMMUTABLE;


CREATE OR REPLACE FUNCTION closest_point_to_line(line geography(LineString), point geography(Point))
RETURNS geography(Point) AS $BODY$
DECLARE
    num_points integer;
    lv1 RECORD;
    lv2 RECORD;
    pv RECORD;
    f RECORD;
    g RECORD;
    h RECORD;
BEGIN
    num_points := ST_NumPoints(line::geometry);
    IF num_points != 2 THEN
        RAISE EXCEPTION 'Only two points are currently supported, but the line has %.', num_points;
    END IF;

    lv1 := _point_to_cartesian(ST_PointN(line::geometry, 1), 1.0);
    lv2 := _point_to_cartesian(ST_PointN(line::geometry, 2), 1.0);
    pv := _point_to_cartesian(point::geometry, 1.0);

    f := _cross_product(lv1.x, lv1.y, lv1.z, lv2.x, lv2.y, lv2.z);
    g := _cross_product(pv.x, pv.y, pv.z, f.x, f.y, f.z);
    h := _cross_product(f.x, f.y, f.z, g.x, g.y, g.z);

    RETURN _cartesian_to_point(h.x, h.y, h.z)::geography;
END
$BODY$ LANGUAGE plpgsql IMMUTABLE;

Answer 3

You could use a similar hack to the ST_Intersection for geography. Would look like this:

CREATE OR REPLACE FUNCTION st_closestpoint(geography, geography)
     RETURNS geography AS
       $$SELECT  geography(ST_Transform(ST_ClosestPoint(ST_Transform(geometry($1),           
    _ST_BestSRID($1,$2)),ST_Transform(geometry($2), _ST_BestSRID($1,$2)) ),4326)) $$
        LANGUAGE sql IMMUTABLE STRICT
        COST 500;