Getting ST_ClosestPoint for Geography type?

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;

Sadly there isn't a geographic version ofST_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 thepostgresql-python-x.xpackage (wherex.xis 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.

This is the PL/pgSQL version of MerseyViking's code. It also uses PostGIS geographyPointandLineStringtypes 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;

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;