finding ley lines in postgres/postgis.sql

query.sql

-- STEP 1 make all alignments at least 20km
-- get every pub-pub connection ((N^2)/2)-N, where N is number of pubs
-- use increasing order of OSM id to limit to a single line between each pair of pubs
-- as we don't care about the direction
-- we also drop connections from each pub to itself (hence the -N bit)

create table aligns as (
    with    
        pubfrom as (select osm_id, amenity, name, way from planet_osm_point where amenity='pub'),
        pubto as (select osm_id as osm_id2, amenity, name, way from planet_osm_point where amenity='pub')
    SELECT 
        pubfrom.osm_id,
        pubto.osm_id2,
        ST_MakeLine(pubfrom.way, pubto.way) as geom
    FROM 
        pubfrom, 
        pubto
    WHERE 
        pubfrom.osm_id <> pubto.osm_id2
        AND pubfrom.osm_id < pubto.osm_id2
    AND ST_Length(ST_MakeLine(pubfrom.way, pubto.way)) > 20000
)

create index ix_aligns on aligns using gist(geom)

-- get start, end, middle pubs within 20 meters, using a .01% sample of all possible alignments

with 
    linez as (select * from aligns where random() < 0.0001) ,
    pubs as (select * from planet_osm_point where amenity='pub')
select
    linez.osm_id as pub_start_id,
    linez.osm_id2 as pub_end_id,
    pubs.osm_id as pub_id
from 
    linez, pubs
where 
    ST_Distance(pubs.way,linez.geom) < 20 and
    linez.osm_id <> pubs.osm_id and
    linez.osm_id2 <> pubs.osm_id
order by
    linez.osm_id, pubs.osm_id

-- STEP TWO find all pub to pub alignments with at least 4 (where start and end are included)
-- using a 5% random sample. can tweak this for larger/smaller datasets

create table leys as (
    with distances as (
        with 
            linez as (select * from aligns where random() < 0.05) ,
            pubs as (select * from planet_osm_point where amenity='pub')
        select
            linez.osm_id as pub_start_id,
            linez.osm_id2 as pub_end_id,
            pubs.osm_id as pub_id,
            ST_Transform(linez.geom,27700) as geom
        from 
            linez, pubs
        where 
            ST_Distance(pubs.way,linez.geom) < 5 and
            linez.osm_id <> pubs.osm_id and
            linez.osm_id2 <> pubs.osm_id
        order by
            linez.osm_id, pubs.osm_id
    ) select
        distances.pub_start_id, distances.pub_end_id, count(*) as ct, geom
      from distances
      group by distances.pub_start_id, distances.pub_end_id, geom
      having count(*) > 1
      order by ct desc
)
    
 -- STEP THREE export as shp
 pgsql2shp -f ~/infoviz/publeylines/scottishleys.shp osmscotland public.leys