I can plot the population of zip codes in NY with something like GeoRegionValuePlot[GeoEntities[New York, United States (administrative division), "ZIPCode"]->"Population"]

I am looking to cluster the zip codes into groups of equal population, and then display those groups on a map.

I am confident that I can do the last part given lists of entities, but I cannot seem to get the clustering to work.

Additionally, I want to cluster based on location, but I want to restrict the size of the clusters by population. Is that possible, and if so, how?

Possible illustrated result: enter image description here

This is a map of NY's current congressional districts. They have the same population within some amount, and they are contiguous. I am not looking to generate this map per se, but rather other maps that meet those criteria (would also be nice if I can seed clusters if possible). I am aware that there is error with using ZIPCode (they cross state lines) but I'm looking for the process, not the result.

Update: I now have a list of lists of size two, where the first contains the polygon information, and the second contains the population. I did that by using:

EntityValue[GeoEntities[New York, United States (administrative division), "ZIPCode"], {"Polygon", "Population"}]

The problem is that I have no way to define a distance function over those pairs (I can't even seem to get the regionCentriod of the polygon)

  • $\begingroup$ This is an interesting problem but the question reads a bit dry. Perhaps some kind of illustration (graphical or otherwise) of what you want would make it more approachable and increase attention. $\endgroup$
    – Mr.Wizard
    Commented Dec 22, 2014 at 18:24
  • $\begingroup$ @Mr.Wizard is this better? $\endgroup$
    – soandos
    Commented Dec 22, 2014 at 18:32
  • $\begingroup$ I think so; we'll see I guess. :-) $\endgroup$
    – Mr.Wizard
    Commented Dec 22, 2014 at 18:34

1 Answer 1


To address the Centroid issue.

 zip = EntityValue[Entity["AdministrativeDivision", {"NewYork", "UnitedStates"}], 
       EntityProperty["AdministrativeDivision", "ZIPCodes"]];
list = {#, EntityValue[#, {"Polygon", "Population"}]} & /@ zip;

Eliminate Zip codes without Polygons and remove the GeoPosition head to allow us to calculate the centroids.

list2 = Complement[list, 
    Cases[list, {_, {Missing["NotApplicable"], _}}]] /. 
   GeoPosition -> Sequence;
list2[[All, 2, 1]] = First[#] & /@ list2[[All, 2, 1]];
centroids = 
 RegionCentroid[DiscretizeGraphics[#]] & /@ list2[[All, 2, 1]]
(*{{41.2703, -71.9886}, {40.7505, -73.9974}, {40.716, -73.9865}, \
{40.7317, -73.9891}, {40.6944, -74.0158}, {40.7061, -74.0086}, \
{40.7096, -74.013}, {40.7137, -74.0083}}*)

You may get some error messages regarding issues with the graphics primitives as follows:

DiscretizeGraphics::invprim: The graphics primitive FilledCurve[{{Line[{{40.6129,-74.0935},{40.6155,-74.0922},{40.6167,-74.0943},{40.6183,-74.0932},<<36>>,{40.61,-74.0968},{40.6124,-74.0994},{40.6157,-74.0962},{40.614,-74.0963}}]},{Line[{<<1>>}]}}] is not valid. >>

We'll replace Indeterminate coordinates by the position of the first coordinate of the Graphics that represents the Zip Code. For that we use the following code.

y = Flatten@Position[centroids, {Indeterminate, Indeterminate}];
x = list2[[#, 2, 1]] & /@ y;
x = x[[#, 1, 1, 1, 1, 1]] & /@ Range[78];
(centroids[[y[[#]]]] = x[[#]]) & /@ Range@Length@x

Now that you have the centroids of each Zip Code you can proceed with the clustering.

data = {centroids[[#, 1]], centroids[[#, 2]], 
     QuantityMagnitude[list2[[#, 2, 2]]]} & /@ Range@Length@centroids;
clusters = FindClusters[%];
ListPointPlot3D[clusters, PlotStyle -> PointSize[Large]]

enter image description here

Please read the documentation on FindClusters, you may want to build your own custom DistanceFunction to give the weight you want to the distance between centroids vs the population in the Zip Code. You can also establish the number of clusters you want. In this example I left the DistanceFunction as the default, which gives the weight to the population (as you can see on the chart above, the clusters are layered along the Z axis).

Now we can project the clusters into the map

 GeoRegionValuePlot[{list2[[#, 1]], list2[[#, 2, 2]]} & /@ 
   Flatten[Position[centroids, #] & /@ clusters[[level, All, ;; 2]]], 
  GeoBackground -> "StreetMapNoLabels"], {{level, 1}, 1, 
  Length@clusters, 1}, SynchronousInitialization -> False, 
 SynchronousUpdating -> False]

enter image description here

Regarding the bounty questions:

1) is there a way to limit cluster size by population (i.e. each cluster can only contain at most N people)?

Yes, you can build your own cluster analysis and can incorporate whatever conditions you may want to use. Please check out the wikipedia article on Cluster Analysis to guide you in the possible algorithms that you may want to use as the base.

2) is there a way to force the clusters to be contigious?

Yes, Please check the code below. I'm using Atlanta due to network performance issues. A custom function coords[zipA,zipB] will return 1 if the zip codes are contiguous, otherwise 0.

polygons = EntityValue[#, "Polygon"] & /@ zips;
zips = Delete[zips, Position[polygons, _Missing]];
polygons = Delete[polygons, Position[polygons, _Missing]];
coordinates = 
     1] & /@ (Cases[#, {a_, b_} -> {a, b}, Infinity] & /@ 
       Cases[#, GeoPosition[a_] -> a, Infinity] & /@ polygons);
db = {zips[[#, 2]], coordinates[[#]]} & /@ Range@Length@coordinates;
coords[zip_] := Cases[db, {zip, {b__}} -> b]
contiguous[listA_, listB_] := 
 If[Length[Intersection[listA, listB]] > 0, 1, 0]

Test enter image description here

{contiguous[coords["30309"], coords["30327"]], 
 contiguous[coords["30309"], coords["30308"]]}
  • $\begingroup$ Any ideas on how to make the clusters the same size in terms of population? (still looking through the rest) $\endgroup$
    – soandos
    Commented Dec 22, 2014 at 21:22
  • $\begingroup$ And is there a way to project this onto a map? $\endgroup$
    – soandos
    Commented Dec 23, 2014 at 6:55
  • $\begingroup$ Firstly, thank you for this, it is most helpful, and I will bounty when I can (~7 hours). Second, I know I get a distance function in terms of the distance between the centriods, but is there a way to make sure that when I include the population that all my regions are contigious. Is there some function that can test that on the regions for me, and is there a way to integrate that into the distance function (maybe make it return infinity if its not connected?) $\endgroup$
    – soandos
    Commented Dec 23, 2014 at 21:00
  • $\begingroup$ And is there a way to restrict the size of the clusters (say the sum of populations in the cluster is not outside of n +- 1000? $\endgroup$
    – soandos
    Commented Dec 23, 2014 at 21:11
  • $\begingroup$ See bounty for more info $\endgroup$
    – soandos
    Commented Dec 26, 2014 at 2:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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