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

chicagosupermarkets = 

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_}] := 
      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 @@@ 
      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?


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 = 

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]

enter image description here

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 = 

Out[] = 292.428

Compare with the geodesic distance (mileage):

In[] := dG = GeoDistance[

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]

enter image description here

  • $\begingroup$ This is very nifty! When I integrate within each cell, I weight points by the population density of their corresponding census tracts. But given that my census tracts are just polygons with coordinate vertices, I should be able to apply the same stereographic projection to their vertices and regenerate the polygons. The level of error you mention is not a concern at all. I'm going to leave the question open for a little bit to see if there are other interesting options, but thank you immensely for taking the time to offer such a thoughtful solution! $\endgroup$ – Shane Jul 22 '17 at 14:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.