How can I efficiently generate regions from long lists of conditions?

Question

I have OSM coordinates on the location of 260 Chicagoland supermarkets:

chicagosupermarkets = 
  Join[
    Uncompress[BarcodeRecognize[Import["https://i.stack.imgur.com/kSqkM.png"]]], 
    Uncompress[BarcodeRecognize[Import["https://i.stack.imgur.com/kWFxR.png"]]]]

Note that the coordinates are listed with longitudes first, which is nonstandard but useful for simple 2D plots.

What I essentially want to do is to calculate the Voronoi cell of each of the supermarkets, but I cannot use the builtin VoronoiMesh or VoronoiDiagram because I need to use actual mileages instead of just Euclidean coordinate distances (which are not accurate measures), and the built-in functions don't take a distance measure as an optional argument. I calculate the mileage between two coordinates as follows:

spheremiles[{θ1_, ϕ1_}, {θ2_, ϕ2_}] := 
  3959 
    InverseHaversine[
      Haversine[ϕ1 Degree - ϕ2 Degree] + 
        Cos[ϕ1 Degree] Cos[ϕ2 Degree] Haversine[θ1 Degree - θ2 Degree]]

For reference, this is the same as the standard implementation listed here but for the fact that it takes the longitude as the first input and returns miles instead of kilometers.

Then I can, in principle, generate the cells as follows:

cells = 
  And @@@ 
    Table[
      spheremiles[chicagosupermarkets[[i]], {x, y}] <= 
        spheremiles[chicagosupermarkets[[j]], {x, y}], 
      {i, Length[chicagosupermarkets]}, 
      {j, Complement[Range[Length[chicagosupermarkets]], {i}]}];

And this does work to some extent: if I'm patient and careful, I can look at a particular Voronoi cell:

RegionPlot[cells[[1]], {x, -88, -87.5}, {y, 41.8, 42.2}]

But cells[[1]] is still a list of 259 conditions, not a well-defined region. I want to make it a region so that I can integrate over it (to find its centroid). Because I'm doing this to eventually implement Lloyd's algorithm, I need to be able to generate these regions and work with them as efficiently as possible. The regions don't have to be all that precise -- discretizations would be fine.

I've also tried Region[ImplicitRegion[cells[[1]], {x, y}]], but this tanked my computer -- I had to restart.

How do I generate regions for these Voronoi cells as efficiently as possible?

Answer 1

A possible way of dealing with your problem is by choosing a cartographic projection such that the geodesic distance between two points of the Earth can be well approximated by the Euclidean distance computed with the projected coordinates within our region of interest. Let us see how this works:

Start by importing the data:

chicagosupermarkets = 
Join[
 Uncompress[BarcodeRecognize[Import["https://i.stack.imgur.com/kSqkM.png"]]],
 Uncompress[BarcodeRecognize[Import["https://i.stack.imgur.com/kWFxR.png"]]]
];

Create a map showing the data:

map=GeoGraphics[GeoPosition /@ Reverse /@ chicagosupermarkets];

We can now extract the cartographic projection used in the map:

In[] := Options[map, GeoProjection]

Out[] = {GeoProjection -> {"Mercator", "Centering" -> GeoPosition[{41.9326, -87.9079}], "GridOrigin" -> {-87.9079, 0}}}

Change the default projection used:

In[]:= projdata = Append[Last@Last@%, "ReferenceModel" ->  
GeodesyData["ITRF00", "SemimajorAxis"]] /. "Mercator" -> "Stereographic"

Out[]={"Stereographic", "Centering" -> GeoPosition[{41.9326,-87.9079}], 
"GridOrigin" -> {-87.9079, 0}, "ReferenceModel" ->   
Quantity[6.37814*10^6, "Meters"]}

New map with the new projection:

GeoMarker /@ GeoPosition /@ Reverse /@ chicagosupermarkets;
GeoGraphics[%, GeoProjection -> projdata, GridLines -> Automatic, 
Frame -> True]

We can now compute the list of the positions of the supermarkets in the projected coordinates:

chicagosupermarketsStereographic = GeoGridPosition[#, projdata] & /@ GeoPosition /@ Reverse /@ chicagosupermarkets;

The geodesic distance (mileage) between supermarkets can be reasonably approximated by the Euclidean distance if we take as the Euclidean coordinates those computed with the stereographic projection. Example:

In[]:= dE = 
EuclideanDistance[
First@chicagosupermarketsStereographic[[1]],
First@chicagosupermarketsStereographic[[2]]
]

Out[] = 292.428

Compare with the geodesic distance (mileage):

In[] := dG = GeoDistance[
chicagosupermarketsStereographic[[1]], 
chicagosupermarketsStereographic[[2]]
]

Out[]= Quantity[292.067, "Meters"]

In[]:= Quantity[dE, "meters"] - dG
Out[]= Quantity[0.360356, "Meters"]

As we see the error is less than half a meter (note that the Euclidean distance computed from the data in the stereographic projection is always given in meters). If we take supermarkets longer apart the error will increase. It is possible to give an estimation of the largest value for the error in the region we are working if we compute the greatest distance within that region:

In[]:= bounds = GeoBounds[GeoPosition /@ chicagosupermarkets]
Out[]= {{-88.7755, -87.0403}, {41.2891, 42.5761}}

Geodesic distance:

In[]:= gD = GeoDistance[
{bounds[[2, 1]], bounds[[1, 1]]}, {bounds[[2, 2]], 
bounds[[1, 2]]}
]
Out[]= Quantity[202.834, "Kilometers"]

Euclidean distance in stereographic projected coordinates:

In[]:= GeoPosition@
{{bounds[[2, 1]], bounds[[1, 1]]}, {bounds[[2, 2]], 
bounds[[1, 2]]}};
In[]:= GeoGridPosition[%, projdata];
In[]:= eD = EuclideanDistance @@ %[[1]]
Out[]= 202911.

Comparison with the geodesic distance:

In[]:= Quantity[eD, "Meters"] - gD
Out[]= Quantity[76.3978, "Meters"]

Therefore the biggest error we should expect is 76 meters which will happen in the calculation of the distance between supermarkets separated 200 km apart. You will have to decide whether this is acceptable for your analysis.

With all this information in mind, we proceed now to the computation of the Voronoi mesh:

In[]:= First /@ chicagosupermarketsStereographic;
In[]:= VoronoiMesh[%, Frame -> True, Axes -> True]