Vrftsjtryk

Suppose we define the "area of control" (AOC) for a single point $P_1$ as a circle of radius $r=500$ meters $C_1$. This could, for example, define the "customer base" of a small bodega in a big city. The area of this circle is $pi r^2 = pi *500^2 $ or about $ 0.7854$km$^2$.

Now suppose we have multiple points (bodegas, if you will) ${P_1, P_2, ..., P_n} in mathbb{R}^2$. We define the "area of control" for each point $P_k$ as the set of points $p in C_k cap mathbb{R}^2$ s.t. $||p-P_k|| < ||p-P_j||, forall jneq k$. In other words, the intersection of the Thiessen polygon for the point $P_k$ with a circle of $500$m around that point. Here is an example of what that might look like with quite a few points (h/t Reddit user /u/TheGoodShipArbitrary):

Is there a quick way to calculate--or even approximate--the area of each AOC given the common radius $r$ of each $C_k$ and the coordinates of the points? It seems that each AOC is entirely determined by $r$ and the list of coordinates.

Background: I am currently working on a coding project that requires finding the geographical "clientele" for sets of establishments which may fall close enough to each other to "interfere" with one another. My current solution involves using GIS (ArcPy) scripts in Python to calculate Thiessen polygons and intersect them with circles, but this is fairly slow. Ultimately, a fast numerical approximation (even a rough one) would be "good enough" for my needs.

An approach using Python would be ideal, but any help is appreciated.

For $r=infty$ you get a Voronoi diagram. There is extensive literature on how to compute these, often with a detour through the corresponding Delaunay triangulation. Different algorithms have different strengths and weaknesses, in terms of code simplicity, time complexity, numeric stability, incremental updates and so on. I'd start with that Voronoi diagram, then make sure to intersect every Voronoi cell with a disk of radius $r$.

On a second read, I see this appears to be mostly what you are doing already. (“Thiessen polygon” had been unfamiliar to me so far.) So let's investigate room for improvement.

Firstly for accurate results. The Voronoi diagram appears hard to avoid: this is essentially a combinatorial problem, and small changes in point placement can lead to huge changes in resulting area. And for big radii, the problem is equivalent to Voronoi diagrams, so in general it can't become simpler if you throw smaller radii into the mix (even though for some specific problem instances it might in fact become simpler). So you could hope that ArcPy does things less efficient than possible, and thus switching implementation will help speed things up. I guess I'd try writing a proof-of-concept implementation e.g. in C++ using CGAL, and see whether that does noticably better than ArcPy. If so, you can dig down into ways to improve this. If not, I wouldn't waste too much time trying more alternatives.

In terms of approximation, here are things that come to mind:

• If ArcPy operates on the sphere, projecting to the map first might help.


• Approximating the disk as a regular $n$-gon might help.


• Compute Voronoi cell for a given point only using other points no more than $2r$ away. As this requires different computations for each cell, you end up with $O(n)$ Voronoi computations, each taking $O(nlog n)$ time, so a net loss. But for specific arrangements this might still help, particularly if most regions do not interfere. Alas, your example picture suggests this will likely not be the case.

You might also want to benchmark whether the Voronoi triangulation or the shape intersections are your bottleneck, and then focus on addressing that.