3

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
;
9
  • Update 1: It took about 45 minutes to run on my Intel® Core™ i3-3217U CPU @ 1.80GHz × 4, 4GB Mem, Ubuntu 13.04. I think I'll need to clean up some of the self-joins by using a Postgres window function (stackoverflow.com/questions/13325583/…). I'm also a little scared that my ST_DWithin function is doing something weird cause the points are all on major highways: geocommons.com/maps/297082#
    – Matt Moehr
    Oct 5 '13 at 18:51
  • What's the diffrence between table2 and table3? Oct 6 '13 at 8:24
  • table2 picks the minimum distance for each point. table2 is based on table1. table1 and table3 are exactly the same, so basically table2 is just there to pick the shortest distance, which is what i would call a "self-join" and is what i'm hoping to optimize a bit with the window function in that link in Update 1. thoughts?
    – Matt Moehr
    Oct 6 '13 at 17:45
  • 1
    Read up on CTE, put indexes on ST_Transform functions if you haven't already. Consider a correlated query instead of a cross-join: gis.stackexchange.com/questions/71592/… Oct 7 '13 at 8:01
  • Yeah I reworked the spatial columns a bit and they don't need ST_Transform anymore. I thought I was doing a pretty good job of using common table expressions (CTE)? I tried following the advice in this good post: workshops.boundlessgeo.com/postgis-intro/…
    – Matt Moehr
    Oct 8 '13 at 18:12
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.

3
  • I assume by FME, you mean: safe.com/fme/fme-technology/fme-desktop/overview
    – Matt Moehr
    Oct 5 '13 at 18:32
  • That's not really going to work because of the software costs, but it looks useful in general.
    – Matt Moehr
    Oct 5 '13 at 18:34
  • Your question about FME, yes, the software by Safe... If it is a one off issue, you could request a two week trial but otherwise, sorry I couldn't be more help. Oct 6 '13 at 0:52
1

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.

5
  • -> Nested Loop (cost=0.00..1343261576.55 rows=1 width=706) Join Filter: ((st_closestpoint(l.geom, dade_random_points.wkb_geometry) && st_expand(dade_random_points.wkb_geometry, 200::double precision)) AND (dade_random_points.wkb_geometry && st_expand(st_closestpoint(l.geom, dade_random_points.wkb_geometry), 200::double precision)) AND _st_dwithin(st_closestpoint(l.geom, dade_random_points.wkb_geometry), dade_random_points.wkb_geometry, 200::double precision) AND st_intersects(st_buffer(st_closestpoint(l.geom, dade_random_points.wkb_geometry), 1e-05::double precision), l.geom))
    – Matt Moehr
    Oct 8 '13 at 18:08
  • This "works" but the join is still too slow to really implement this as a solution in my processing flow. Relevant lines from explain are in my above comment. I guess I'm going to have to write a function that does incremental table inserts...
    – Matt Moehr
    Oct 8 '13 at 18:10
  • I assume you have spatial indexes on your geometry columns; did you try adding a subquery to just get the row with rank 1?
    – djq
    Oct 8 '13 at 18:44
  • Yes on spatial indexes. Do you mean: select * from <your query> where ranked_distance = 1; ??? Cause that doesn't help performance at all. I think the loop still has to go over all roads for each point.
    – Matt Moehr
    Oct 8 '13 at 19:08
  • maybe ORDER BY rank and LIMIT 1 might help (not sure - have not tried)
    – djq
    Oct 8 '13 at 19:42
1

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
;
1
  • 1
    Performance is diffrent not better or worse. Correlated queries tend to be slower but they are more memory efficient. Oct 9 '13 at 12:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.