7
\$\begingroup\$

EDIT: My code DOES not work, see https://gis.stackexchange.com/questions/294380/if-rectangle-corner-points-have-same-nearest-neighbor-does-whole-interior

Given 5 (currently hardcoded) cities, this code uses Google Maps to split the world into 5 regions, each region being the points closer to a 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 improved

use POSIX;

# defining constants here is probably bad
$PI = 4.*atan(1);
$EARTH_RADIUS = 6371/1.609344; # miles

open(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 recurse

sub 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 and
latitude/longitude u,v in miles Source: http://williams.best.vwh.net/avform.htm

=cut

sub 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 approach

EDIT: Thanks to whoever upvoted this. I've now fixed the http://test.barrycarter.info/gmap8.php link and am including a screenshot below:

enter image description here

\$\endgroup\$
0

1 Answer 1

Reset to default
4
\$\begingroup\$

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.

\$\endgroup\$