points by mkl 6 years ago

No I don't. All you need, though, is the ability to draw straight lines (well, great circle lines) between points. You could do the region generation once offline, in Python:

    import numpy as np
    from scipy.spatial import SphericalVoronoi

    import json

    with open('volcanos.json', 'rt') as f:
        volcanos = json.load(f)
        #There are some nulls in the data:
        volcanos = [v for v in volcanos if v['longitude'] != None and v['latitude'] != None]

    def longlat_to_xyz(longlat):
        longlat = longlat.reshape((-1, 2))*(np.pi/180.)
        r = np.cos(longlat[:, 1])
        return np.column_stack((r*np.cos(longlat[:, 0]), r*np.sin(longlat[:, 0]), np.sin(longlat[:, 1])))
    
    #go through set to remove duplicates:
    volcano_longlat = np.array(list(set((v['longitude'], v['latitude']) for v in volcanos)))
    volcano_xyz = longlat_to_xyz(volcano_longlat)

    sv = SphericalVoronoi(volcano_xyz)
    sv.sort_vertices_of_regions() #make the vertices go around each shape

    #find unique edges (they're shared between regions):
    edges = set()
    for region in sv.regions:
        region.append(region[0])
        for i in range(len(region)-1):
            edges.add((region[i], region[i+1]))

    long_min = -30.
    long_max = 330.
    def xyz_to_longlat(xyz, long_min=long_min, long_max=long_max):
        long = np.arctan2(xyz[:, 1], xyz[:, 0])*(180./np.pi)
        long[long<long_min] += 360.
        long[long>long_max] -= 360.
        return np.column_stack((long, np.arcsin(xyz[:, 2])*(180./np.pi)))

    vertices = xyz_to_longlat(sv.vertices)

    #To draw the regions on Google Maps etc., save edges and vertices as JSON and load into JS.

    #To draw the linked map with Matplotlib (very slow for hi-res, and it was fiddly setting up the right conda environment):
    import matplotlib.pyplot as plt
    from mpl_toolkits.basemap import Basemap

    fig, ax = plt.subplots(figsize=(16, 10))
    m = Basemap(projection='robin', lon_0=-200, resolution='h', ax=ax)
    m.drawcoastlines(linewidth=.1)
    m.drawcountries(linewidth=.1)
    volcano_longlat2 = xyz_to_longlat(volcano_xyz)
    x, y = m(volcano_longlat2[:, 0], volcano_longlat2[:, 1])
    m.plot(x, y, 'ro', markersize=.5)

    for edge in edges:
        m.drawgreatcircle(vertices[edge[0], 0], vertices[edge[0], 1],
                          vertices[edge[1], 0], vertices[edge[1], 1],
                          color='b', linewidth=.3)
    fig.tight_layout()
    fig.savefig('volcano_voronoi.png', dpi=600)