Given a table of points and a table of lines, snap all the points within 10 metres of the lines to the lines.My first thought was "PostGIS doesn't have that snapping function", but it actually does, hidden in the linear-referencing functions: ST_Line_Locate_Point(line, point).
OK, that returns a measure along the line, but I want a point! No problem, ST_Line_Interpolate_Point(line, measure) returns a point from a measure.
Great, so now all I need are, for each point within 10 metres of the lines, the nearest line. Yuck, finding the minimum. However, with the PostgreSQL
DISTINCT ON syntax and some ordering, it all pops out:SELECT
DISTINCT ON (pt_id)
pt_id,
ln_id,
ST_AsText(
ST_line_interpolate_point(
ln_geom,
ST_line_locate_point(ln_geom, vgeom)
)
)
FROM
(
SELECT
ln.the_geom AS ln_geom,
pt.the_geom AS pt_geom,
ln.id AS ln_id,
pt.id AS pt_id,
ST_Distance(ln.the_geom, pt.the_geom) AS d
FROM
point_table pt,
line_table ln
WHERE
ST_DWithin(pt.the_geom, ln.the_geom, 10.0)
ORDER BY
pt_id,d
) AS subquery;
The sub-query finds all the points/line combinations that meet the 10 meter tolerance rule, and returns them in sorted order, by point id and distance. The outer query then strips off the first entry for each distinct point id and runs the LRS functions on it to derive the new snapped point.
Snapperiffic!
9 comments:
Nifty...
Don't you have to ORDER the subquery by the distance in order to "ensure" that the DISTINCT ON picks the closest line?
Yes, that's right, so order by pt_id (which I'll be distincting on) and then by d, which is distance.
Ah right, I didn't see the "d".
This is a good one for the PostGIS wiki - or that "Expert PostGIS" that someone should write!
This doesn't quite look right. Don't you have your distinct on in the wrong place. I would think it should be - and I still like inner join over ,. Something about , just always smelled funny to me.
SELECT
pt_id,
ln_id,
ST_AsText(
ST_line_interpolate_point(
ln_geom,
ST_line_locate_point(ln_geom, vgeom)
)
)
FROM
(
SELECT DISTINCT ON (pt_id)
ln.the_geom AS ln_geom,
pt.the_geom AS pt_geom,
ln.id AS ln_id,
pt.id AS pt_id,
ST_Distance(ln.the_geom, pt.the_geom) AS d
FROM
point_table pt INNER JOIN
line_table ln
ON
ST_DWithin(pt.the_geom, ln.the_geom, 10.0)
ORDER BY
pt_id,d
) AS subquery;
Slight correction
DISTINCT ON(pt.id)
I couldn't get the ordering to happen before the DISTINCT ON, so I moved it out a level. If I could get ordering first, I could maybe ditch the whole outer query!
Well it is a neat idea never the less. I'll have to play with it sometime.
I don't think you want to get rid of the subquery even if you could because PostgreSQL is naive in the sense that it will process everything in the SELECT part even if it never pulls the record so you'd be running that expensive interpolate across all the records of a point. You want to do all the expensive stuff after you have the set of data you want. I made yet another typo.
The ORDER BY should have been pt.id, ST_Distance(ln.the_geom, pt.the_geom)
So should be
SELECT
pt_id,
ln_id,
ST_AsText(
ST_line_interpolate_point(
ln_geom,
ST_line_locate_point(ln_geom, pt_geom?)
)
)
FROM
(
SELECT DISTINCT ON (pt_id)
ln.the_geom AS ln_geom,
pt.the_geom AS pt_geom,
ln.id AS ln_id,
pt.id AS pt_id
FROM
point_table pt INNER JOIN
line_table ln
ON
ST_DWithin(pt.the_geom, ln.the_geom, 10.0)
ORDER BY
pt.id,ST_Distance(ln.the_geom, pt.the_geom)
) AS subquery;
(if you made the mistake I did when first trying it, maybe that was the problem)
Also where is vgeom defined? I just noticed that I don't know where that is coming from and assume it must be pt_geom.
Okay one more try and then I'm going to give up and admit I can't function without a debugger.
SELECT
pt_id,
ln_id,
ST_AsText(
ST_line_interpolate_point(
ln_geom,
ST_line_locate_point(ln_geom, pt_geom?)
)
)
FROM
(
SELECT DISTINCT ON (pt.id)
ln.the_geom AS ln_geom,
pt.the_geom AS pt_geom,
ln.id AS ln_id,
pt.id AS pt_id
FROM
point_table pt INNER JOIN
line_table ln
ON
ST_DWithin(pt.the_geom, ln.the_geom, 10.0)
ORDER BY
pt.id,ST_Distance(ln.the_geom, pt.the_geom)
) AS subquery;
Hi, great post.
But I've a problem: points do not come from a table but from a client, how can I use this feature?
thanks
Post a Comment