Skip to content

Instantly share code, notes, and snippets.

@jsanz
Last active January 9, 2022 23:57
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 2 You must be signed in to fork a gist
  • Save jsanz/bcd7d6ae7f8e2c2ee3d7 to your computer and use it in GitHub Desktop.
Save jsanz/bcd7d6ae7f8e2c2ee3d7 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