Message Boards Message Boards

[WSS19] Cartogram Function

<em>A visualisation of 500 iterations on a map of contiguous US states

A visualisation of 500 iterations on a map of contiguous US states


Introduction

The goal of this project is to create a function, which generates cartograms. Cartograms are geometrically distorted maps, where the distortion conveys a piece of information obtained from the due dataset.

Currently, more than 25 different algorithms for creating distortion can be found in the literature. As I really enjoy working with matrices, I decided to use the Cellular Automata method1. This algorithm was developed in 1990 by Daniel Dorling, an English social geographer and a professor of Geography at the University of Oxford. The limitations of this approach are that it does not work properly for non-contiguous areas and that it preserves the global borders. However, it is relatively easy to implement in the Wolfram Language and the output is reliable.



Algorithm

enter image description here

As checking the topology condition would be computationally too difficult, I decided to omit it. I will talk about the implications of this decision later.


Code

Cartogram[geoPolysi:{__Rule},numberOfIterations_:30,res_:100, plotAll_?BooleanQ]:=
    Module [
        {regions,masks,map,mapPairs,population,numberOfCells,cellDensity,gPairs,geoPolys},

        geoPolys = Keys[geoPolysi];

        regions = Range[Length[EntityValue[geoPolys,"Polygon"]]];

        map = (
            masks = Function [
                poly, ImageData [
                    Binarize @ Rasterize [
                        Style [
                            Graphics [
                                {White, MeshPrimitives[DiscretizeGraphics @ poly, 2]},
                                Background -> Black,
                                PlotRange -> Reverse @ GeoBounds [
                                    EntityValue[geoPolys,"Polygon"]
                                    ]
                                ],
                            Antialiasing -> False
                        ],
                        RasterSize -> res
                    ],
                    "Bit"
                ]
            ] /@ EntityValue[geoPolys, "Polygon"];

            Fold [
                #1 + Clip [
                        #2-(#2*#1),
                        {0,Infinity}
                    ]&,
                    MapIndexed [
                        #1 * First[#2]&,
                        masks
                    ]
            ]
        );

        mapPairs = BlockMap [
            {
                #[[2,2]],
                3*Count [
                    Position [
                        Flatten[#],
                        #[[2,2]]
                    ],
                    _?EvenQ,
                    2
                ] + Count [
                    Position [
                        Flatten[#],
                        #[[2,2]]
                    ],
                    _?OddQ,
                    2
                ] -1 
            }&,
            ArrayPad[map, 1, 0], 
            {3, 3}, 
            {1, 1}
        ];

        population = QuantityMagnitude[Values[geoPolysi]];

        numberOfCells = (Count[map, #, 2]) & /@ Range[Length[EntityList[geoPolys]]];

        cellDensity = N[Prepend[population / numberOfCells, 0]];

        gPairs = mapPairs /. {a_Integer, b_Integer} -> g[a, b];

        (
            f [{
                    {g[a1_, a2_], g[b1_, b2_], g[c1_, c2_]},
                    {g[d1_, d2_], g[e1_, e2_], g[f1_, f2_]},
                    {g[g1_, g2_], g[h1_, h2_], g[i1_, i2_]}
                }] :=

                Module [
                    {highest, nsecurity},
                    If [
                        numberOfCells[[e1]] > 1 && e1 =!= 0,

                        highest = First [
                            Sort [
                                {a1, b1, c1, d1, f1, g1, h1, i1},
                                cellDensity[[#1+1]] > cellDensity[[#2+1]]&
                            ]
                        ];
                        nsecurity = RandomChoice [
                            Cases [
                                {{a1, a2}, {b1, b2}, {c1, c2}, {f1, f2}, {i1, i2}, {h1, h2}, {g1, g2}, {d1, d2}},
                                {highest, s_} -> s,
                                {1},
                                1
                            ]
                        ];
                        If [
                            (cellDensity[[highest+1]] > cellDensity[[e1+1]] && e2 < 12) || e2 < 8,
                            g[highest, e2],
                            g[e1, e2]
                        ],

                        g[e1,e2]
                     ]
                ]
        );
        (
            regionStep[gPairs_] := 
                BlockMap [
                    f,
                    ArrayPad[gPairs, 1, g[0, 0]],
                    {3,3},
                    {1,1}
                ];

            recalcSecurities[gPairs_] := 
                BlockMap [
                    g [
                        #[[2,2]],
                        3*Count [
                            Position [
                                Flatten[#],
                                #[[2,2]]
                            ],
                            _?EvenQ,
                            2
                        ] 
                        + 
                        Count [
                            Position [
                                Flatten[#],
                                #[[2,2]]
                            ],
                            _?OddQ,
                            2
                        ] - 1
                    ]&,
                    ArrayPad[gPairs[[All, All, 1]], 1, "-"],
                    {3,3},
                    {1,1}
                ];

            recalcDensities[gPairs_] := Module [
                {numberOfCells},
                numberOfCells = (Count[gPairs[[All, All, 1]], #, 2])& /@ regions;
                cellDensity = N[Prepend[population / (numberOfCells + 1), 0]]
            ];
        );

        steps1 = NestList [
            (
                recalcDensities[#];
                recalcSecurities[regionStep[#]]
            )&,
            gPairs,
            numberOfIterations
        ];
        If[
           TrueQ[plotAll],
           ArrayPlot[#[[All, All, 1]]]& /@ steps1,
           ArrayPlot[Last[steps1][[All, All, 1]]]
           ]
    ];

Cartogram[geoPolysi:{__Rule},numberOfIterations_:30,res_:100]:=Cartogram[geoPolysi,numberOfIterations,res, False]

Cartogram[geoPolys_, opts___] := 
    Cartogram [
        Normal [
            DeleteMissing[
               EntityValue [
                   geoPolys,
                   "Population",
                   "EntityAssociation"
            ]
            ]
        ],
        opts
    ];

The function has got 4 arguments, 3 of which are optional and have default values. The first argument can takes an entity of the type Country. Initially, I was planning to use the Wolfram language built-in function CellularAutomata, but the function is not suitable due to delayed updating of values during iterations.


Output

Generating a plot of a region. When not given any specific values for countries, the function automatically calculates population densities and distorts the map based on that. enter image description here Running the function for a list of countries. enter image description here Running the function for a list of countries using the user-input values. enter image description here Running the function for just one iteration. enter image description here Running the function for 5 iterations with the width of 400 pixels. enter image description here Returns all steps of first 5 iterations. enter image description here


As promised in the Algorithm section, now I will explain the consequences of omitting the topology condition. When different regions (countries) reach equal densities, the areas adjacent to borders will keep oscillating. This is due to the fact that there are no restrictions on the number of times each cell in a neighbourhood can change its region, and thus bordering cells will do so during each iteration (as they have low security factors mentioned in the algorithm).

enter image description here Oscillations described above.


References

Dorling, D. (1996). Area cartograms. Norwich: Geo Abstracts Univ. of East Anglia, p.30.


Acknowledgments

Mentor: Christopher Wolfram

POSTED BY: Richard Lapin
4 Replies

@Richard, this is very impressive work and thanks for making the Wolfram Function Repository entry. But which data are plotted on the very top animation for USA? If these are population densities it is a bit strange Florida is so low - should be very high.

POSTED BY: Kapio Letto

enter image description here - Congratulations! This post is now featured in our Staff Pick column as distinguished by a badge on your profile of a Featured Contributor! Thank you, keep it coming, and consider contributing your work to the The Notebook Archive!

POSTED BY: Moderation Team

Very impressive and very useful. Thanks for posting.

Cheers,

Marco

POSTED BY: Marco Thiel

Thank you Marco, glad you find the function useful.

My submission to the Wolfram Function Repository has been approved, you can find it here: https://resources.wolframcloud.com/FunctionRepository/resources/Cartogram

Cheers,

Richard

POSTED BY: Richard Lapin
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract