Skip to content

Instantly share code, notes, and snippets.

@stevefaeembra
Last active August 29, 2015 14:07
Show Gist options
  • Save stevefaeembra/bec7f845f3f33bd7959e to your computer and use it in GitHub Desktop.
Save stevefaeembra/bec7f845f3f33bd7959e to your computer and use it in GitHub Desktop.
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