Skip to content

Instantly share code, notes, and snippets.

@plablo09

plablo09/natural_cities.sql

Last active August 14, 2019 19:35
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save plablo09/2a752bd79befc0e535a397c9a145b480 to your computer and use it in GitHub Desktop.
Save plablo09/2a752bd79befc0e535a397c9a145b480 to your computer and use it in GitHub Desktop.
Download ZIP
Pure SQL implementation of Bin Jiang Natural Cities Algorithm
Raw
natural_cities.sql
-- create delaunay lines
select st_delaunaytriangles(st_collect(the_geom),0.0000001, 1)from planet_osm_roads_vertices_pgr
-- create delaunay and dump the resulting multilinestrings
select (st_dump(st_delaunaytriangles(st_collect(the_geom),0.0000001, 1))).path[1] as id, (st_dump(st_delaunaytriangles(st_collect(the_geom),0.0000001, 1))).geom
from planet_osm_roads_vertices_pgr
-- add edge length in meters (over the spheroid)
select (st_dump(st_delaunaytriangles(st_collect(the_geom),0.0000001, 1))).path[1] as id,
(st_dump(st_delaunaytriangles(st_collect(the_geom),0.0000001, 1))).geom,
st_length((st_dump(st_delaunaytriangles(st_collect(the_geom),0.0000001, 1))).geom::geography) as length
from planet_osm_roads_vertices_pgr
-- store on a table
create table as (
select (st_dump(st_delaunaytriangles(st_collect(the_geom),0.0000001, 1))).path[1] as id,
(st_dump(st_delaunaytriangles(st_collect(the_geom),0.0000001, 1))).geom,
st_length((st_dump(st_delaunaytriangles(st_collect(the_geom),0.0000001, 1))).geom::geography) as length
from planet_osm_roads_vertices_pgr
)
-- select values smaller than average
select * from del_lines
where length < ALL(select avg(length) from del_lines)
-- Polygonize and simplify (dissolve). Sore the result in a table
with a as (
select st_exteriorring((st_dump(ST_Union(geom))).geom) as geom
FROM (select (st_dump(st_polygonize(geom))).path[1] as id, (st_dump(st_polygonize(geom))).geom
from(
select * from del_lines
where length < ALL(select avg(length) from del_lines)
) as foo) as bar
)
select (row_number() over())::integer as id, st_makepolygon(geom) as geom
from a
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment