Given 5 (currently hardcoded) cities, this code uses Google Maps tosplit the world into 5 regions, each region being the points closer toa given city than to the other 4.It works, but is slow in two senses:The program takes a long time to run, especially when $minarea is small.The program generates many polygons, slowing down the Google Map above. With smaller '$minarea', I've even run out of JavaScript stack space.Thought: Would something w/ qhull be faster?The bvoronoi subroutine does most of the work:#!/bin/perl# Unusual approach to Voronoi diagram of Earth sphere: cut into 4# pieces and compare closest point for 4 vertices# TODO: this program is long and clumsy and can doubtless be improveduse POSIX;# defining constants here is probably bad$PI = 4.*atan(1);$EARTH_RADIUS = 6371/1.609344; # milesopen(A,>/home/barrycarter/BCINFO/sites/TEST/gvorbin.txt);# latitude and longitude of points%points = ( Albuquerque => 35.08 -106.66, Paris => 48.87 2.33, Barrow => 71.26826 -156.80627, Wellington => -41.2833 174.783333, Rio => -22.88  -43.28);# primartish colors%colors = ( Albuquerque => #ff0000, Paris => #00ff00, Barrow => #0000ff, Wellington => #ffff00, Rio => #ff00ff, BORDER => #000000);# stop at what gridsize$minarea = .5;# the four psuedo-corners of the globe$nw = bvoronoi(0,90,-180,0);$ne = bvoronoi(0,90,0,180);$sw = bvoronoi(-90,0,-180,0);$se = bvoronoi(-90,0,0,180);for $i (split(\n,$nw\n$ne\n$sw\n$se)) {  # create google filled box  my($latmin, $latmax, $lonmin, $lonmax, $closest) = split(/\s+/, $i);  # build up the coords  print A << MARK;var myCoords = [ new google.maps.LatLng($latmin, $lonmin), new google.maps.LatLng($latmin, $lonmax), new google.maps.LatLng($latmax, $lonmax), new google.maps.LatLng($latmax, $lonmin), new google.maps.LatLng($latmin, $lonmin)];myPoly = new google.maps.Polygon({ paths: myCoords, strokeColor: $colors{$closest}, strokeOpacity: 1, strokeWeight: 0, fillColor: $colors{$closest}, fillOpacity: 0.5});myPoly.setMap(map);MARK;}# workhorse function: given a square (on an equiangular map),# determine the closest point of 4 vertices; if same, return square# and point; otherwise, break square into 4 squares and recursesub bvoronoi {  # Using %points as global is ugly  my($latmin, $latmax, $lonmin, $lonmax) = @_;  my($mindist, $dist, %closest);  # compute distance to each %points for each corner  # TODO: this is wildly inefficient, since I just need relative  # distance, not exact!  for $lat ($latmin,$latmax) {    for $lon ($lonmin,$lonmax) {      # TODO: has to be easier way to do this?      $mindist = 0; $dist= 0;      for $point (keys %points) {        my($plat,$plon) = split(/\s+/, $points{$point});        $dist = gcdist($lat, $lon, $plat, $plon);        if ($dist < $mindist || !$mindist) {          $mindist = $dist;          $minpoint = $point;        }      }      # this point is closest to one vertex of the square      # TODO: should abort loop if we already have two different closest points      $closest{$minpoint} = 1;    }  }  # if there's just one point closest to all four corners, return it  my(@keys) = keys %closest;  # if @keys has length 1, return it  unless ($#keys) {    return $latmin $latmax $lonmin $lonmax $keys[0];  }  # if we've hit a border point, return it (area too small)  my($area) = ($latmax-$latmin)*($lonmax-$lonmin);  if ($area <= $minarea) {    return $latmin $latmax $lonmin $lonmax BORDER;  }  # split square and recurse  my($latmid) = ($latmax+$latmin)/2.;  my($lonmid) = ($lonmax+$lonmin)/2.;  my(@sub) = ();  push(@sub, bvoronoi($latmin, $latmid, $lonmin, $lonmid));  push(@sub, bvoronoi($latmid, $latmax, $lonmin, $lonmid));  push(@sub, bvoronoi($latmin, $latmid, $lonmid, $lonmax));  push(@sub, bvoronoi($latmid, $latmax, $lonmid, $lonmax));  return join(\n, @sub);}=item gcdist($x,$y,$u,$v)Great circle distance between latitude/longitude x,y andlatitude/longitude u,v in miles Source: http://williams.best.vwh.net/avform.htm=cutsub gcdist {    my(@x)=@_;    my($x,$y,$u,$v)=map {$_*=$PI/180} @x;    my($c1) = cos($x)*cos($y)*cos($u)*cos($v);    my($c2) = cos($x)*sin($y)*cos($u)*sin($v);    my($c3) = sin($x)*sin($u);    return ($EARTH_RADIUS*acos($c1+$c2+$c3));}Another failed approachEDIT: Thanks to whoever upvoted this. I've now fixed the http://test.barrycarter.info/gmap8.php link and am including a screenshot below:
Spherical Voronoi diagram, binary splitting approach
performance;perl;geospatial;google maps
In short, try cssgrid from NCAR Graphics.At length:There's a similar question on StackOverflow, but it has no accepted answer.So here's an answer, if you'll permit two apologies. 1. I'll leave the Google Maps work to you. 2. I didn't read your code thoroughly.The algorithm's running time should be no higher in order than sorting, O(n log n). At least, Voronoi diagrams on the sphere, openly available as a Utrecht University report, claims that that result is in Kevin Quintin Brown's dissertation. Furthermore On the construction of the Voronoi mesh on a sphere claims in its abstract to construct the mesh in O(n).It looks like qhull doesn't do all that you need, but I found something else. The cssgrid package is found in the GPL-licensed NCAR Graphics software and does what you need, with C, Fortran, and NCAR-specific bindings. You can download all of NCAR Graphics from its page at the (United States) National Center for Atmospheric Research. If you wish, you can use SWIG to make the relevant C binding available to Perl.The cssgrid package has its own documentation page online. The functions *csvoro* (Voronoi) and *css2c* (spherical to cartesian, that is, latitude-longitude to X-Y-Z coordinates) are the most relevant ones.The cssgrid package is based on STRIPACK by permission of STRIPACK's author. STRIPACK comes up higher in web searches and is cost-free for noncommercial use. However, STRIPACK itself as an ACM TOMS algorithm is neither public domain nor GPL-compatible. Also it is Fortran 90 code without a provided C binding, although that difficulty might be surmountable using FortWrap.