In a previous post (I feel like I say that a lot), I wrote about rotating address points to match nearby roads in replicate the effect of USGS quads that represented small buildings with little squares that followed the nearby road alignment.

The function was effective:

ALTER TABLE jstein.address_points_all ADD COLUMN azimuth integer;
UPDATE jstein.address_points_all addr
SET azimuth = degrees(ST_Azimuth(addr.the_geom, ST_ClosestPoint(road.the_geom, addr.the_geom)))
FROM jstein.street_centerlines road;

but deadly slow when applied to all 500,000 address points. And so we iterate. First, I’ll show you our reasonably fast, but elegantly written solution. Later I’ll have a post on our very fast, but somewhat hackerish approach.

A reminder of our target aesthetic:

And our new code:

CREATE TABLE address_points_rot AS
WITH index_query AS
(
  SELECT
    addr.*, ST_Distance(addr.the_geom, road.the_geom) AS distance,
    degrees(ST_Azimuth(addr.the_geom, ST_ClosestPoint(road.the_geom, addr.the_geom))) as azimuth
  FROM address_pointssub as addr, street_centerlines as road
  WHERE ST_DWithin(addr.the_geom, road.the_geom, 2000) = true
  ORDER BY addr.the_geom <-> road.the_geom
)

SELECT DISTINCT ON(cartodb_id) * FROM index_query ORDER BY gid, distance

In the interest of full disclosure, I stole the structure of this from Paul Ramsey’s post on K Nearest Neighbor searches in PostGIS. AFAIK, you need Postgresql 9.1 and Postgis 2.0 or later to do this… . For 500,000 points this took 16 minutes. Narrowing that search distance would help a lot. I suspect there is a much better way to structure this using true KNN for each point, but my SQL-fu failed me here… . Suggestions welcome.