First we create a function to do our angle calculations, then use select to loop through all the records and do the calculations.
Within our function, first we find our first (in this case) five nearest streets using our distance operator , based on bounding boxes, then we further order by ST_Distance, with a LIMIT 1 to get just the closest street. By the way, this is about twice as fast as the last approach.
And our new code:
CREATE OR REPLACE FUNCTION angle_to_street (geometry) RETURNS double precision AS $$ WITH index_query as (SELECT ST_Distance($1,road.geom) as dist, degrees(ST_Azimuth($1, ST_ClosestPoint(road.geom, $1))) as azimuth FROM street_centerlines As road ORDER BY $1 <#> road.geom limit 5) SELECT azimuth FROM index_query ORDER BY dist LIMIT 1; $$ LANGUAGE SQL; DROP TABLE IF EXISTS address_points_rot; CREATE TABLE address_points_rot AS SELECT addr.*, angle_to_street(addr.geom) FROM address_points addr;
I have to credit Alexandre Neto with much of the code for this solution.