compute the voronoi diagram of a set of points given in geojson format

gj-voronoi.py

#!/usr/bin/python3
# coding=utf-8
# -*- coding: utf-8 -*-
# example output at umap: https://umap.openstreetmap.fr/zh-tw/map/python3-gj-voronoipy-tkecgeojson_774525

import argparse, scipy, warnings, json, scipy.spatial, math
import numpy as np

def voronoi(data):
    if type(data) is dict:
        if data['type'] == 'FeatureCollection':
            points = data['features']
        else:
            warnings.warn('voronoi can only process a geojson FeatureCollection or a geojson array of points')
            return {}
    elif type(data) is list:
        points = data
    else:
        warnings.warn('voronoi can only process a geojson FeatureCollection or a geojson array of points')
        return {}
    scaled_coords = np.array([ pt['geometry']['coordinates'] for pt in points ])
    (max_x, min_x) = (max(scaled_coords[:,0]), min(scaled_coords[:,0]))
    (max_y, min_y) = (max(scaled_coords[:,1]), min(scaled_coords[:,1]))
    width = max_x - min_x
    height = max_y - min_y
    mid_x = (max_x + min_x) / 2
    mid_y = (max_y + min_y) / 2
    scale = math.cos(mid_y/180*math.pi)
    scaled_coords[:,0] = (scaled_coords[:,0] - mid_x) * scale
    vor = scipy.spatial.Voronoi(np.array(scaled_coords))
    ret = []
    for point_index in range(len(vor.points)):
        # assumption: point indices/ordering in vor.points remain the same as input.
        assert(
            abs(scaled_coords[point_index][0] - vor.points[point_index][0])/width +
            abs(scaled_coords[point_index][1] - vor.points[point_index][1])/height < 1e-7
        )
        ret.append(points[point_index])
        reg_index = vor.point_region[point_index]
        any_infinity = False
        reg_coords = []
        for i in vor.regions[reg_index]:
            if i < 0:
                any_infinity = True
            else:
                scale_back = (vor.vertices[i][0] / scale + mid_x, vor.vertices[i][1])
                reg_coords.append(scale_back)
        if any_infinity:
            reg_feature = {
                'type':'Feature',
                'properties': points[point_index]['properties'].copy(),
                'geometry': {
                    'type': 'LineString',
                    'coordinates': reg_coords
                }
            }
        else:
            reg_coords.append(reg_coords[0])
            reg_feature = {
                'type':'Feature',
                'properties': points[point_index]['properties'].copy(),
                'geometry': {
                    'type': 'Polygon',
                    'coordinates': [reg_coords]
                }
            }
        if 'name' in points[point_index]['properties']:
            reg_feature['properties']['name'] += ' voronoi region'
        ret.append(reg_feature)
    return ret if type(data) is list else { 'type':'FeatureCollection', 'features':ret }

parser = argparse.ArgumentParser(
    description='compute the voronoi diagram of a geojson file of points',
    formatter_class=argparse.ArgumentDefaultsHelpFormatter)
parser.add_argument('-a', '--area', type=str,
    default='area', help='field name of area')
parser.add_argument('-c', '--centroid', type=str,
    default='centroid', help='field name of centroid')
parser.add_argument('FILE', help='the geojson file')
args = parser.parse_args()

with open(args.FILE, 'r') as gjfile:
    data = json.load(gjfile)

# https://stackoverflow.com/questions/18337407/saving-utf-8-texts-in-json-dumps-as-utf8-not-as-u-escape-sequence
r = voronoi(data);
print(json.dumps(r, ensure_ascii=False, indent=2)) #.encode('utf8'))

tkec.geojson

numpy.asscalar(a) is deprecated since NumPy v1.16 and fails in newer versions of Python.

I changed the following line:

vert_sites = [set([np.asscalar(s) for e in elist for s in vor.ridge_points[e]]) for elist in vert_ridges]

to

vert_sites = [set([s.item() for e in elist for s in vor.ridge_points[e]]) for elist in vert_ridges]

now the script works again

Thanks for the fix, @Robert-112 . Now I see more problems with my program. For one thing, the library certainly would assume that distances between points is sqrt(x*x+y*y). But this formula is incorrect for points given as longitude and latitude. One should convert each longitude x to x*cos(th) where th is the rough latitude of the entire data set (assuming the points are close enough to use the same th for all points as an approximation), call the voronoi library, and then convert the longitude values back. That part is easy to fix. Secondly, and more problematically, my program does not correctly process vertices at infinity. I was trying my program on a small data set such as A=(0,0) B=(1.732,0) C=(0,1) D=(-1,1) and it produces an empty list. I looked at the new version of the scipy.spatial library and have not found an easy way to fix that yet.

OK, I uploaded a test data file tkec.geojson . It contains the locations of all branches of tk3c, a major electronic chain store in Taiwan. I also updated the code to implement the cos(th) fix. The output now contains both the voronoi regions and the input points. The result is displayed at this umap Closed voronoi regions are output in the geojson file as "Polygon" features and should now be correctly displayed. Open voronoi regions are output in the geojson file as LineString features and may not be quite correct.