finding ley lines in postgres/postgis.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