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;

3 Answers 3


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]))

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);
(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.

  • Excellent! Your calculated point appears to lie on the great circle, because even ST_DWithin(0.00000000001) returns true. Just one question: is there any reason Python is needed? It doesn't seem to use any math functions that PostgreSQL doesn't have natively.
    – EM0
    Commented Nov 24, 2011 at 0:26
  • 1
    OK, I ported this to plpgsql and tried it out. It generally work very well, but there are some cases it doesn't handle - see edit. Also, it only works for a line with 2 points - is there a way to extend it for a LINESTRING with multiple points (other than just trying the line between every 2 points)?
    – EM0
    Commented Nov 24, 2011 at 5:34
  • The only reason I used Python was because my knowledge of PL/pgSQL is all but non-existant, whereas I use Python quite a lot. Kudos to you for converting it. Any chance you could post the code for future reference? I'm not sure why it doesn't work if the point is already on the line, you may need to do a distance check first. And yes, the way the code stands, you'll need to check every segment one at a time. You may be able to tap into the spatial index to eliminate the majority of segments, but that's beyond my ken. Commented Nov 24, 2011 at 13:28
  • Just posted the PL/pgSQL code. OK, I think I could extend it to handle multiple points (even if naively), but I'm not sure how to fix that incorrect result. The point is not already on the line in that case, but the closest point is the start or end point of the line. Of course, I could just add a hack to check for that specific case, but I'm reluctant to do that without really understanding the problem - who knows if it fixes all cases.
    – EM0
    Commented Nov 25, 2011 at 0:29

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)
    lon float;
    lat float;
    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);

CREATE OR REPLACE FUNCTION _cartesian_to_point(x float, y float, z float)
RETURNS geometry(Point) AS $BODY$
    lon float;
    lat float;
    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);

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)
    x := y1 * z2 - z1 * y2;
    y := z1 * x2 - x1 * z2;
    z := x1 * y2 - y1 * x2;

CREATE OR REPLACE FUNCTION closest_point_to_line(line geography(LineString), point geography(Point))
RETURNS geography(Point) AS $BODY$
    num_points integer;
    lv1 RECORD;
    lv2 RECORD;
    pv RECORD;
    f RECORD;
    g RECORD;
    h RECORD;
    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;

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)) $$
        COST 500;
  • All that does is cast between geography and geometry, right? The closest point is still calculated on the plane, not on the sphere.
    – EM0
    Commented Nov 29, 2011 at 9:18
  • The closest point is calculated on the plane, but it transforms the geography to the in theory best planar coordinate for the region and then transforms back after computing the closest. So its not a straight cast. The secret sauce as Paul would say is the _ST_BestSRID call that finds the best spatial ref compromise given the two input geographies. That's how the ST_Intersection for geography is implemented. Function on this page with a grin use this approach postgis.org/documentation/manual-svn/…
    – Regina Obe
    Commented Nov 29, 2011 at 19:07
  • you can use the _ST_BestSRID?
    – ziggy
    Commented Dec 8, 2017 at 13:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.