Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

Edit: 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;
share|improve this question
add comment

3 Answers

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.

share|improve this answer
    
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. –  Reality Exists Nov 24 '11 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)? –  Reality Exists Nov 24 '11 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. –  MerseyViking Nov 24 '11 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. –  Reality Exists Nov 25 '11 at 0:29
add comment

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;
share|improve this answer
add comment

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;
share|improve this answer
    
All that does is cast between geography and geometry, right? The closest point is still calculated on the plane, not on the sphere. –  Reality Exists Nov 29 '11 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/… –  LR1234567 Nov 29 '11 at 19:07
    
I see - thank you. I may try that out later. –  Reality Exists Dec 2 '11 at 0:19
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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