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 | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment