Tuesday, April 08, 2008

Snapping Points in PostGIS

Fun question on the #postgis IRC channel today, just hard enough to be interesting and just easy enough to not be overwhelming:
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:

Dr JTS said...

Nifty...

Don't you have to ORDER the subquery by the distance in order to "ensure" that the DISTINCT ON picks the closest line?

Paul Ramsey said...

Yes, that's right, so order by pt_id (which I'll be distincting on) and then by d, which is distance.

Dr JTS said...

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!

Regina Obe said...

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;

Regina Obe said...

Slight correction

DISTINCT ON(pt.id)

Paul Ramsey said...

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!

Regina Obe said...

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.

Regina Obe said...

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;

markux said...

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

About Me

My Photo
Victoria, British Columbia, Canada

Followers

Blog Archive

Labels

bc (32) it (26) postgis (17) icm (11) enterprise IT (9) sprint (9) open source (8) osgeo (8) video (8) management (6) cio (5) enterprise (5) foippa (5) gis (5) spatial it (5) foi (4) mapserver (4) outsourcing (4) bcesis (3) foss4g (3) oracle (3) politics (3) COTS (2) architecture (2) boundless (2) esri (2) idm (2) natural resources (2) ogc (2) open data (2) opengeo (2) openstudent (2) postgresql (2) rant (2) technology (2) vendor (2) web (2) 1.4.0 (1) HR (1) access to information (1) accounting (1) agile (1) aspen (1) benchmark (1) buffer (1) build vs buy (1) business (1) business process (1) cathedral (1) cloud (1) code (1) common sense (1) consulting (1) contracting (1) core review (1) crm (1) custom (1) data warehouse (1) deloitte (1) design (1) digital (1) email (1) essentials (1) evil (1) exadata (1) fcuk (1) fgdb (1) fme (1) foocamp (1) foss4g2007 (1) ftp (1) gds (1) geocortex (1) geometry (1) geoserver (1) google (1) google earth (1) government (1) grass (1) hp (1) iaas (1) icio (1) industry (1) innovation (1) integrated case management (1) introversion (1) iso (1) isss (1) isvalid (1) javascript (1) jts (1) lawyers (1) mapping (1) mcfd (1) microsoft (1) mysql (1) new it (1) nosql (1) opengis (1) openlayers (1) oss (1) paas (1) pirates (1) policy (1) portal (1) proprietary software (1) qgis (1) rdbms (1) recursion (1) regression (1) rfc (1) right to information (1) saas (1) salesforce (1) sardonic (1) seibel (1) sermon (1) siebel (1) snark (1) spatial (1) standards (1) svr (1) tempest (1) texas (1) tired (1) transit (1) twitter (1) udig (1) uk (1) uk gds (1) verbal culture (1) victoria (1) waterfall (1) wfs (1) where (1) with recursive (1) wkb (1)