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 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
x.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 * v2 - v1 * v2, v1 * v2 - v1 * v2, v1 * v2 - v1 * v2) 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, h, h) return (math.degrees(nearest_point), math.degrees(nearest_point)) $$ 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 geography
LineStringtypes 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;