Skip to content

Instantly share code, notes, and snippets.

@mgottholsen
Forked from jsanz/README.md
Created August 11, 2017 12:26
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 mgottholsen/95c9b9bd69db58c7d48a79e0ee260afa to your computer and use it in GitHub Desktop.
Save mgottholsen/95c9b9bd69db58c7d48a79e0ee260afa to your computer and use it in GitHub Desktop.
Analysis example: Voronoi on CartoDB

This example loads a CartoDB layer using a crazy SQL from a small variation of this awesome Stack Overflow response. Note the use of CSS style tags for the CartoCSS to allow easy editing.

The uncompressed (and not perfect!) version of the SQL to draw the Voronoi diagram is:

WITH 
    Sample AS (
      SELECT st_setsrid(st_union(the_geom),0) as geom 
      FROM registro_centros_nz WHERE spanish = 1
    ),
    Edges AS (
    SELECT id,
        UNNEST(ARRAY['e1','e2','e3']) EdgeName,
        UNNEST(ARRAY[
            ST_MakeLine(p1,p2) ,
            ST_MakeLine(p2,p3) ,
            ST_MakeLine(p3,p1)]) Edge,
        ST_Centroid(ST_ConvexHull(ST_Union(
            ST_CurveToLine(REPLACE(
                ST_AsText(ST_LineMerge(ST_Union(ST_MakeLine(p1,p2),ST_MakeLine(p2,p3)))),
                'LINE','CIRCULAR'),15),
            ST_CurveToLine(REPLACE(
                ST_AsText(ST_LineMerge(ST_Union(ST_MakeLine(p2,p3),ST_MakeLine(p3,p1)))),
                'LINE','CIRCULAR'),15)
        ))) ct      
    FROM (
        SELECT id,
            ST_PointN(g,1) p1,
            ST_PointN(g,2) p2,
            ST_PointN(g,3) p3
        FROM    (
            SELECT (gd).Path id, ST_ExteriorRing((gd).Geom) g 
            FROM (SELECT (ST_Dump(ST_DelaunayTriangles(geom))) gd FROM Sample) a 
            )b
        ) c
    ), 
    Voronoi AS(
    SELECT *
    FROM (
        SELECT  
            ST_LineMerge(ST_Union(ST_MakeLine(
            x.ct,
            CASE 
            WHEN y.id IS NULL THEN
                CASE WHEN ST_Within(
                    x.ct, (SELECT ST_ConvexHull(geom) FROM sample)) 
                THEN 
                    ST_MakePoint(
                        ST_X(x.ct) + ((ST_X(ST_Centroid(x.edge)) - ST_X(x.ct)) * 2),
                        ST_Y(x.ct) + ((ST_Y(ST_Centroid(x.edge)) - ST_Y(x.ct)) * 2))
                END
            ELSE 
                y.ct
            END
            ))) v
        FROM Edges x 
            LEFT OUTER JOIN 
            Edges y ON x.id <> y.id AND ST_Equals(x.edge,y.edge)
        ) z
    ), 
    Convex as (
      SELECT 
        ST_Buffer(
            ST_ConvexHull(
              ST_Transform(
                ST_SetSrid(sample.geom,4326)
               ,3857)
            )
       ,20) as geom
      FROM sample
    ) 
SELECT 
    int4(row_number() OVER (ORDER BY v)) AS cartodb_id,
    ST_Intersection(
        convex.geom,
        ST_Transform(ST_SetSRID(v,4326),3857)
    ) AS the_geom_webmercator
FROM voronoi,convex
<!doctype html>
<html>
<head>
<meta charset="utf-8">
<meta http-equiv="X-UA-Compatible" content="IE=edge,chrome=1">
<title>NZ Spanish schools Voronoi</title>
<link rel="shortcut icon" href="http://cartodb.com/assets/favicon.ico" />
<meta name="viewport" content="width=device-width, initial-scale=1">
<link rel="stylesheet" href="http://libs.cartocdn.com/cartodb.js/v3/3.14/themes/css/cartodb.css" />
<style type="cartocss/text" id="cartocss-voronoi">
#registro_centros_nz::back{
line-color: #FF9900;
line-width: 6;
line-opacity: 0.8;
}
#registro_centros_nz{
line-color: #FF5C00;
line-width: 3;
line-opacity: 0.8;
}
</style>
<style type="cartocss/text" id="cartocss-schools">
#registro_centros_nz{
marker-fill: #D6301D;
marker-width: 10;
marker-line-color: #FFFFFF;
marker-line-width: 1.5;
marker-line-opacity: 1;
marker-fill-opacity: 0.9;
marker-type: ellipse;
marker-placement: point;
marker-allow-overlap: true;
[zoom < 14] {
marker-comp-op: multiply;
}
}
</style>
<style>
html, body, #map {
height: 100%;
padding: 0;
margin: 0;
}
</style>
</head>
<body>
<div id="map"></div>
<script src="https://cartodb-libs.global.ssl.fastly.net/cartodb.js/v3/3.14/cartodb.js"></script>
<script>
$( document ).ready(function() {
//Create the leaflet map
var map = L.map('map', {
zoomControl: true,
center: [-36.8697,174.8069],
zoom: 11
});
//Add basemap
var basemap = L.tileLayer('http://{s}.basemaps.cartocdn.com/light_nolabels/{z}/{x}/{y}.png', {
attribution: '&copy; <a href="http://www.openstreetmap.org/copyright">OpenStreetMap</a> contributors, &copy; <a href="http://cartodb.com/attributions">CartoDB</a>'
}).addTo(map);
//Add the layers
cartodb.createLayer(map, {
user_name: 'vehrka',
type: 'cartodb',
sublayers: [
{
// here be dragons, be brave!!
sql: "WITH Sample AS (SELECT st_setsrid(st_union(the_geom),0) as geom from registro_centros_nz where spanish = 1 ), Edges AS (SELECT id, UNNEST(ARRAY['e1','e2','e3']) EdgeName, UNNEST(ARRAY[ST_MakeLine(p1,p2) , ST_MakeLine(p2,p3) , ST_MakeLine(p3,p1)]) Edge, ST_Centroid(ST_ConvexHull(ST_Union(ST_CurveToLine(REPLACE(ST_AsText(ST_LineMerge(ST_Union(ST_MakeLine(p1,p2),ST_MakeLine(p2,p3)))),'LINE','CIRCULAR'),15), ST_CurveToLine(REPLACE(ST_AsText(ST_LineMerge(ST_Union(ST_MakeLine(p2,p3),ST_MakeLine(p3,p1)))),'LINE','CIRCULAR'),15) ))) ct FROM (SELECT id, ST_PointN(g,1) p1, ST_PointN(g,2) p2, ST_PointN(g,3) p3 FROM (SELECT (gd).Path id, ST_ExteriorRing((gd).Geom) g FROM (SELECT (ST_Dump(ST_DelaunayTriangles(geom))) gd FROM Sample) a )b ) c ), voronoi AS(SELECT * FROM (SELECT ST_LineMerge(ST_Union(ST_MakeLine(x.ct, CASE WHEN y.id IS NULL THEN CASE WHEN ST_Within(x.ct, (SELECT ST_ConvexHull(geom) FROM sample)) THEN ST_MakePoint(ST_X(x.ct) + ((ST_X(ST_Centroid(x.edge)) - ST_X(x.ct)) * 2),ST_Y(x.ct) + ((ST_Y(ST_Centroid(x.edge)) - ST_Y(x.ct)) * 2)) END ELSE y.ct END ))) v FROM Edges x LEFT OUTER JOIN Edges y ON x.id <> y.id AND ST_Equals(x.edge,y.edge) ) z ), convex as (SELECT ST_Buffer(ST_ConvexHull(ST_Transform(ST_SetSrid(sample.geom,4326) ,3857) ) ,20) as geom from sample ) select int4(row_number() OVER (ORDER BY v)) AS cartodb_id, ST_Intersection(convex.geom, ST_Transform(ST_SetSRID(v,4326),3857) )AS the_geom_webmercator from voronoi,convex",
cartocss: $("#cartocss-voronoi").html()
},
{
sql: 'SELECT * FROM registro_centros_nz where spanish = 1',
cartocss: $("#cartocss-schools").html()
}]
}).addTo(map);
});
</script>
</body>
</html>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment