D3 vs Scipy (Voronoi diagram implementation)

Background

I'm working with a set of 8000 geographical points contained in csv file. On one hand I create a visualisation of Voronoi diagrams built using these points - it's done using D3 library. On the other hand I calculate these Voronoi diagrams in Python using Scipy.

My work logic is simple - I mess with my data on Python's side, making heatmaps, analysis and so on and then I visualise effects using D3. But today I accidentally found that Voronoi diagrams made by Scipy and D3 are different. I noticed that after using geojson.io to plot GeoJsons of Voronois made in Python just to see if I can visualise everything there.

As I said, the Voronois were different - some of them had different angles and some even had additional vertices.

Question:

Why is that happening? Why Voronoi diagrams calculated by these 2 libraries (D3 and Scipy) differ?

Further description

How it is done on D3 side: Based on Chris Zetter example http://chriszetter.com/blog/2014/06/15/building-a-voronoi-map-with-d3-and-leaflet/ I translate latitude and longitude into custom projection to visualise it on the mapbox map.

var voronoi = d3.geom.voronoi()
.x(function(d) { return d.x; })
.y(function(d) { return d.y; })
.clipExtent([[N_W.x , N_W.y],[S_E.x, S_E.y]])

I create Voronoi based on points that are visible within map border + some padding (filteredPoints)

  filteredPoints = points.filter(function(d) {
  var latlng = new L.LatLng(d.latitude, d.longitude);
  if (!drawLimit.contains(latlng)) { return false };
  // this translates points from coordinates to pixels
  var point = map.latLngToLayerPoint(latlng);
  key = point.toString();
  if (existing.has(key)) { return false };
  existing.add(key);
  d.x = point.x;
  d.y = point.y;
  return true;
  });

voronoi(filteredPoints).forEach(function(d) { d.point.cell = d});

How it is done on Python side: I use scipy.spatial.Voronoi.

from scipy.spatial import Voronoi

def create_voronois():
    points = numpy.array(points_list)
    vor = Voronoi(points)

Where "points_list" is a list of my 8000 geographical points.

EDIT:

Screenshot from my visualisation - black borders are Voronois made with D3, white ones are made by scipy.spatial.Voronoi. As we can see scipy is wrong. Did anyone compare these 2 libraries before?

https://i.sstatic.net/Csehd.jpg

Code to run. It prints GeoJson with badly calculated Voronois.

import numpy
from scipy.spatial import Voronoi
from geojson import FeatureCollection, Feature, Polygon

points = [
[22.7433333333000,  53.4869444444000],
[23.2530555556000,  53.5683333333000],
[23.1066666667000,  53.7200000000000],
[22.8452777778000,  53.7758333333000],
[23.0952777778000,  53.4413888889000],
[23.4152777778000,  53.5233333333000],
[22.9175000000000,  53.5322222222000],
[22.7197222222000   ,53.7322222222000],
[22.9586111111000,  53.4594444444000],
[23.3425000000000,  53.6541666667000],
[23.0900000000000,  53.5777777778000],
[23.2283333333000,  53.4713888889000],
[23.3488888889000,  53.5072222222000],
[23.3647222222000   ,53.6447222222000]]

def create_voronois(points_list):
    points = numpy.array(points_list)
    vor = Voronoi(points)
    point_voronoi_list = []
    feature_list = []
    for region in range(len(vor.regions) - 1):
        vertice_list = []
        for x in vor.regions[region]:
            vertice = vor.vertices[x]
            vertice = (vertice[1], vertice[0])
            vertice_list.append(vertice)
            polygon = Polygon([vertice_list])
            feature = Feature(geometry=polygon, properties={})
            feature_list.append(feature)

    feature_collection = FeatureCollection(feature_list)
    print feature_collection

create_voronois(points)

Solution

Apparently your javascript code is applying a transformation to the data before computing the Voronoi diagram. This transformation does not preserve the relative distances of the points, so it does not generate the same result as your scipy code. Note that I'm not saying that your d3 version is incorrect. Given that the data are latitude and longitude, what you are doing in the javascript code might be correct. But to compare it to the scipy code, you have to do the same transformations if you expect to get the same Voronoi diagram.

The scripts below show that, if you preserve the relative distance of the input points, scipy's Voronoi function and d3.geom.voronoi generate the same diagram.

Here's a script that uses scipy's Voronoi code:

import numpy
from scipy.spatial import Voronoi, voronoi_plot_2d
import matplotlib.pyplot as plt


points = [
[22.7433333333000,  53.4869444444000],
[23.2530555556000,  53.5683333333000],
[23.1066666667000,  53.7200000000000],
[22.8452777778000,  53.7758333333000],
[23.0952777778000,  53.4413888889000],
[23.4152777778000,  53.5233333333000],
[22.9175000000000,  53.5322222222000],
[22.7197222222000,  53.7322222222000],
[22.9586111111000,  53.4594444444000],
[23.3425000000000,  53.6541666667000],
[23.0900000000000,  53.5777777778000],
[23.2283333333000,  53.4713888889000],
[23.3488888889000,  53.5072222222000],
[23.3647222222000,  53.6447222222000]]


vor = Voronoi(points)

voronoi_plot_2d(vor)
plt.axis('equal')
plt.xlim(22.65, 23.50)
plt.ylim(53.35, 53.85)
plt.show()

It generates this plot:

Now here's a javascript program that uses d3.geom.voronoi:

<html>
<head>
    <script type="text/javascript" src="http://mbostock.github.com/d3/d3.js"></script>
    <script type="text/javascript" src="http://mbostock.github.com/d3/d3.geom.js"></script>
</head>
<body>
  <div id="chart">
  </div>
  <script type="text/javascript">

  // This code is a hacked up version of http://bl.ocks.org/njvack/1405439

  var w = 800,
      h = 400;

  var data = [
    [22.7433333333000,  53.4869444444000],
    [23.2530555556000,  53.5683333333000],
    [23.1066666667000,  53.7200000000000],
    [22.8452777778000,  53.7758333333000],
    [23.0952777778000,  53.4413888889000],
    [23.4152777778000,  53.5233333333000],
    [22.9175000000000,  53.5322222222000],
    [22.7197222222000,  53.7322222222000],
    [22.9586111111000,  53.4594444444000],
    [23.3425000000000,  53.6541666667000],
    [23.0900000000000,  53.5777777778000],
    [23.2283333333000,  53.4713888889000],
    [23.3488888889000,  53.5072222222000],
    [23.3647222222000,  53.6447222222000]
  ];

  // Translate and scale the points.  The same scaling factor (2*h) must be used
  // on x and y to preserve the relative distances among the points.
  // The y coordinates are also flipped.
  var vertices = data.map(function(point) {return [2*h*(point[0]-22.5), h - 2*h*(point[1]-53.4)]})

  var svg = d3.select("#chart")
    .append("svg:svg")
      .attr("width", w)
      .attr("height", h);
  var paths, points;

  points = svg.append("svg:g").attr("id", "points");
  paths = svg.append("svg:g").attr("id", "point-paths");

  paths.selectAll("path")
      .data(d3.geom.voronoi(vertices))
    .enter().append("svg:path")
      .attr("d", function(d) { return "M" + d.join(",") + "Z"; })
      .attr("id", function(d,i) { 
        return "path-"+i; })
      .attr("clip-path", function(d,i) { return "url(#clip-"+i+")"; })
      .style("fill", d3.rgb(230, 230, 230))
      .style('fill-opacity', 0.4)
      .style("stroke", d3.rgb(50,50,50));

  points.selectAll("circle")
      .data(vertices)
    .enter().append("svg:circle")
      .attr("id", function(d, i) { 
        return "point-"+i; })
      .attr("transform", function(d) { return "translate(" + d + ")"; })
      .attr("r", 2)
      .attr('stroke', d3.rgb(0, 50, 200));

  </script>
</body>
</html>

It generates:

Based on a visual inspection of the results, I'd say they are generating the same Voronoi diagram.