Best way to assign points to nearest roads in postgis?

Question

I'm trying to implement a simplified type of dasymetric mapping for census data. I've downloaded the Census's American Community Survey individual-level data and randomly assigned every household to a point within a block group. Now I want to move these random points to the nearest road. I've downloaded the OSM data and loaded it via osm2pgsql. So I have two tables:

dade_random_points: a point for each household. 'wkb_geometry' is the geo field.

dade_planet_osm_roads: a linestring for every road. 'way' is the geo field.

(By the way, I'm doing a proof-of-concept with just Miami-Dade County for right now.)

After reading How to get the nearest point on a linestring to a given point? and http://www.postgis.org/docs/ST_ClosestPoint.html

I wrote this enormous, multi-nested query. It's still churning on my laptop, but it seems to work when I put in 'limit' statements to do smaller pieces.

Surely there has to be a better way? (Or maybe there's a postgis dasymetric library out there somewhere that I haven't been able to find?)

drop table if exists dade_points_on_roads cascade;
create table dade_points_on_roads as
select
  table4.id,
  table4.osm_id,
  table4.geoid10,
  table4.name,
  ST_ClosestPoint( 
    ST_Transform(table4.way, 4326), 
    ST_Transform(table4.wkb_geometry, 4326)
    ) as point_on_road
from
  (
  select
    --table3.*,
    points.*,
    roads.*
  from
    (
    select
      table1.id,
      min(table1.distance) as distance
    from
      (
      select
        p.id,
        r.osm_id,
        ST_Distance(
          ST_Transform(r.way, 4326),
          ST_Transform(p.wkb_geometry, 4326)
          ) as distance
      from 
        dade_planet_osm_roads as r,
        dade_random_points as p
      where 
        ST_DWithin( 
          ST_Transform(r.way, 4326), 
          ST_Transform(p.wkb_geometry, 4326), 
          500
          )  
      --limit 1000
      ) as table1
    --order by distance asc
    group by table1.id
    --limit 10
    ) as table2,
    (
    select
      p.id,
      r.osm_id,
      ST_Distance(
        ST_Transform(r.way, 4326),
        ST_Transform(p.wkb_geometry, 4326)
        ) as distance
    from 
      dade_planet_osm_roads as r,
      dade_random_points as p
    where 
      ST_DWithin( 
        ST_Transform(r.way, 4326), 
        ST_Transform(p.wkb_geometry, 4326), 
        500
        )  
    --limit 1000
    ) as table3,
    dade_planet_osm_roads as roads,
    dade_random_points as points
  where
    table2.id = table3.id
    and table2.distance = table3.distance
    and table2.id = points.id
    and table3.osm_id = roads.osm_id
  ) as table4
;

Answer 1

I acknowledge it isn't PostGIS, but a very easy solution would be to use FME for this process. FME plays nice (as it was specifically intended) with hundreds of formats, including PostGIS/PostGre.

My two cents, if you do use or have access to it, I would simply connect up a NeighborFinder transformer with the Points (ie the houses) to the Candidate port and the lines (ie the roads) to the Base. Enter a max search distance (ie 500 ft) and it will "snap" the points to the line. You will need to add a 2D Point Replacer and re-generate the Points. Then you can dump them into whatever database/file you would like.

I routinely work with datasets in the hundreds of thousands and the smaller the search radius the faster it works.

Answer 2

I don't full understand your query but one thing that you can do is ORDER BY distance and then LIMIT 1 to get the nearest road segment to a point. ORDER BY gives you the results in ascending order by default.

Considering a different approach, you could try the following:

• find a point on the closest lines within a certain distance
• buffer that point a tiny amount and intersect it with the line to select it
• rank the distances for each line

Here is a suggested implementation:

SELECT
 p.id as point_id,
 l.id as line_id,
 l.geom,
 ROW_NUMBER() OVER(PARTITION BY l.id ORDER BY ST_Distance(ST_ClosestPoint(l.geom, p.geom), p.geom)) as ranked_distance
FROM
 point as p,
 lines as l
WHERE
 ST_Intersects( ST_Buffer(ST_ClosestPoint(l.geom, p.geom), 0.00001), l.geom)
AND 
 ST_DWithin(ST_ClosestPoint(l.geom, p.geom), p.geom, 1000);

Then, to get all of the nearest lines you would select everything that has a rank equal to 1. This could be a subquery in the above query if you wanted to do it in one step.

Answer 3

Based on @Jakub recommendations I re-wrote this with a correlated query instead of self-join. The performance is maybe a little better (maybe not?), but it's cleaner and easier to get additional fields in the outermost select statement. I think I will work on a function that wraps this basic query and does a separate "insert into" for each point_on_road instead of the "create table" approach. Yes that will be slow and klunky, but I can build in better logging and turn the thing lose on the whole country. And if it fails, I'll be able to restart it at the point of failure.

drop table if exists dade_points_on_roads cascade;
create table dade_points_on_roads as
select 
  p.id,
  p.geoid10, 
  l.osm_id,
  l.name, 
  l.geom as road,
  p.wkb_geometry as household,
  ST_ClosestPoint(l.geom, p.wkb_geometry) as point_on_road
from
  ( 
  select 
    p.id as g1,
    (
    select l.osm_id
    from dade_planet_osm_roads as l
    order by ST_Distance(l.geom, p.wkb_geometry) ASC 
    limit 1
    ) as g2
  FROM 
    dade_random_points AS p
  --limit 10
  OFFSET 0
  ) AS table1
  inner join dade_random_points as p
  ON table1.g1=p.id
  inner join dade_planet_osm_roads as l
  ON table1.g2=l.osm_id
;