Creating an Animated Cartogram (QGIS3)

Cartogram is a type of map visualization where the shape of each feature is distorted in proportion to a variable. Cartograms make it easy to see large variations in the data. The simplest method to create a cartogram is by scaling the size of each region according to a variable. This method retains the original shape of the polygon and only changes the size. These are known as Non-contiguous isomorphic cartograms. In this tutorial, we will learn how to use QGIS expressions to create a cartogram and use the Temporal Controller to create an animation that gradually transforms the features to the target size.

This tutorial is builds on Hans van der Kwast's excellent cartogram tutorial. The formula for scale factor used here is from the original paper Noncontiguous Area Cartograms by Judy M Olson.

Overview of the task

We will take a layer of states in the US and create an animated cartogram by scaling each state by population. The resulting map will have each state's area proportional its population.

Get the data

United States Census Bureau provides cartographic boundary files along with demographic datasets. We will download the data and process them to create a data layer suitable for our task.

1. Visit the Cartographic Boundary Files - Shapefile data page on the US Census Bureau website. Scroll down and download the States shapefile cb_2018_us_state_20m.zip. This has the state polygons we need for the visualization.

2. The state polygons do not have any demographic data. This data needs to be downloaded separately and joined with the shapefile for use in a GIS. Visit the State Population Totals and Components of Change: 2020-2023 page and download the Annual Population Estimates, Estimated Components of Resident Population Change, and Rates of the Components of Resident Population Change for the United States, States, District of Columbia, and Puerto Rico: April 1, 2020 to July 1, 2023 dataset which will be downloaded as NST-EST2023-ALLDATA.csv.

3. Open QGIS. Locate the cb_2018_us_state_20m.zip file in the QGIS Browser and expand it. Drag and drop the cb_2018_us_state_20m.shp layer to the canvas.

4. A new layer cb_2018_us_state_20m will be added to the Layers panel. Next, we will load the CSV file. Click the Open Data Source Manager button.

5. Switch to the Delimited Text tab. Click ... next to File name and browse to the downloaded NST-EST2023-ALLDATA.csv file. Expand the Geometry Definition section and select No geometry (attribute only table). Click Add.

6. A new layer NST-EST2023-ALLDATA will be added to the Layers panel. Right-click and select Open Attribute Table. This table has the population counts in the POPESTIMATE2023 column. Each state has a unique id in the STATE column that we will use to join this table with the polygon layer.

7. Right-click the cb_2018_us_state_20m layer and select Open Attribute Table. The state ids are contained in the GEOID column. While the numbers are the same, they are formatted as a 2-digit 0-padded number. To match these with the STATE column from our population table, we need similarly formatted numbers. Select the NST-EST2023-ALLDATA. Go to Processing ‣ Toolbox. Search and locate the Vector table ‣ Field calculator algorithm. Double-click to open it.

8. In the Field calculator dialog, select NST-EST2023-ALLDATA as the Input layer. Enter GEOID as the Field name and set the Result field type to Text (string). We will now take the numbers from the STATE field and use the lpad() function to create a 2-digit 0-padded string. Enter the following expression and click Run.

lpad("STATE", 2, '0')

9. A new layer Calculated will be added to the Layers panel. Right-click and select Open Attribute Table. Note that the newly created GEOID column has correctly formatted identifiers. We can now use this field to join this table with the states layer. Search and locate the Vector general ‣ Join attributes by field value algorithm. Double-click to open it.

10. In the Join attributes by field value dialog, select cb_2018_us_state_20m as the Input layer. Select GEOID as the Table field. For Input layer 2, select our table Calculated and Table field 2 select GEOID. The table has many columns but we only need the population for the latest year. Click the ... button for Layer 2 fields to copy and select only the POPESTIMATE2023 field. Leave other options to their default value and click Run.

11. A new layer Joined layer will be added to the Layers panel. Before using this layer for our cartogram, let's reproject it to a projected CRS. Search and locate the Vector general ‣ Reproject layer algorithm. Double-click to open it.

12. In the Reproject layer dialog, select Joined layer as the Input layer. For the Target CRS, click the Select CRS button. Search for the North_America_Albers_Equal_Area_Conic CRS and select it. This is our final layer so we will save it to disk. Click the ... button next to Reprojected and select Save to File....

13. Enter the name of the layer as us_states_with_population.gpkg and select Save. Click Run to create the GeoPackage file with the reprojected data.

We will use this layer in the next section. For convenience, you may directly download a copy of the above layer from below:

us_states_with_population.gpkg

Data Source [USCENSUS]

Procedure

1. Start a new QGIS project. Locate the us_states_with_population.gpkg file in the Browser and expand it. Drag and drop the us_states_with_population layer to the empty canvas.

2. A new layer us_states_with_population will be added to the Layers panel. Right-click and select Open Attribute Table. For our cartogram, we need to use the population density as the variable. We will use the population values in the POPESTIMATE2023 field. Close the attribute table.

3. Search and locate the Vector table ‣ Field calculator algorithm from the Processing Toolbox. Double-click to open it.

4. In the Field calculator dialog, select us_states_with_population as the Input layer. Enter density as the Field name. Enter the following expression to compute the density. As the function area(@geometry) calculates the area in the unit of the CRS (which is in meters), we apply the conversion factor to convert it to square kilometers. Click Run.

1000*1000* "POPESTIMATE2023" /
area(@geometry)

5. Click the ... button next to Calculated and select Save to File.... Enter the name of the layer as us_states_population_density.gpkg and select Save. Click Run.

6. A new layer us_states_population_density will be added to the Layers panel. Right-click and select Open Attribute Table. We need to pick an Anchor feature against which all other features will be scaled. Ideally you would pick the feature with the highest value of the variable you want to use for the cartogram. This will ensure there are no overlapping areas. Double-click the density column header to sort the table by density. You will notice that the highest density value in our dataset is quite large compared to other values and it belongs to a fairly small state. This will result all the features being scaled to very small sizes. We can pick the feature the second highest density which has a relatively large area and its density is comparable to other features.

7. In Area cartograms, the scale factor determines how much the feature's area is reduced. We must reduce each feature's area so that the population density of the feature is the same as the population density of the anchor feature. The formula for scale factor is the ratio of the square-root of the feature's value against the square-root of the value of the anchor feature. Open the Vector table ‣ Field calculator algorithm from the Processing Toolbox. In the Field calculator dialog, select us_states_population_density as the Input layer. Enter scale_factor as the Field name. Enter the following expression to compute the scale factor. The expression calculates the ratio of the square root of the feature's density against the square root of the density of the second largest density value. Click the ... button next to Calculated and select Save to File.... Enter the name of the layer as us_states_scale_factor.gpkg and select Save. Click Run.

sqrt("density")/array_get(array_agg(
   expression:=sqrt("density"),
   order_by:=sqrt("density")), -2)

8. A new layer us_states_scale_factor will be added to the Layers panel. Right-click and select Open Attribute Table. The scale_factor field now contains the ratio by which each feature must be scaled to have the same population density as the anchor feature.

9. We only need the us_states_scale_factor layer for the final visualization. Select the remaining layers, right-click and select Remove Layer.

10. Select the us_states_scale_factor layer and click the Open the layer styling panel button in the Layers panel. Select Simple Fill and open the drop-down selector for Symbol layer type. Set the Symbol layer type to Outline: Simple Line and select a Color of your choice. This symbol layer will be a reference for our map when we resize the polygons.

11. Click the Add Symbol layer (+) button. A new symbol layer Simple Fill will be added. Set the Fill color to the same color as the lines and the Stroke color to be a slightly darker color.

12. Next, open the drop-down selector for Symbol layer type and select Geometry Generator as the Symbol layer type. Geometry generator allows us to modify the geometry for rendering using expressions. Click the Expression Builder button.

13. We will use the scale() function which resizes the given geometry by X- and Y-scaling factors. For our cartogram, we want to resize each polygon by the ratio of its population to the highest population. Enter the following expression to apply this scaling and click OK.

scale(
   @geometry,
   "scale_factor",
   "scale_factor",
   centroid(@geometry)
)

14. You will see the state polygons are now sized by the proportion of each state's population to the highest population. Many large states with low population densities are now much smaller than their original size. You will notice that polygons with irregular shapes are off-center after being scaled. This is because the anchor point of the scaling is the centroid of the geometry which is often the representative point for the polygon. Let's update our expression to fix this. Click the Expression Builder button.

15. The scale() function takes an optional parameter to specify the scaling center point. We will use the pole_of_inaccessibility() function to find a representative anchor point for each polygon. This is similar to a centroid, but it is guaranteed to be inside of the polygon whereas a centroid can fall outside for certain shapes. Update the expression as shown below which calculates the pole of inaccessibility of the geometry with a small tolerance value and click OK.

scale(
    @geometry,
    "scale_factor",
    "scale_factor",
    pole_of_inaccessibility(@geometry, 100)
)

16. Now the scaled polygons will have much better placement. We see another problem. Many features in the layer are Multipolygons, i.e. they have more than one part. Such features have 2 or more polygons that are part of the same geometry. With our current expression, both are scaled with the same anchor point computed from the combined geometry. This is not ideal. For example, a large feature with multiple islands should be scaled such that each island is scaled with its own center point. To fix this, we update our expression to iterate over each part of the geometry and scale it with its own center. Click the Expression Builder button.

17. Here we use the array_foreach() function to iterate over each part of the geometry and generate the scaled versions of them. Finally collect_geometries() function combines each scaled part into a single multi-polygon geometry. Update the expression as shown below and click OK.

collect_geometries(
    array_foreach(generate_series(1, @geometry_part_count),
            scale(geometry_n(@geometry,@element),
                    "scale_factor",
                    "scale_factor",
pole_of_inaccessibility(geometry_n(@geometry,@element), 100)
            )
    )
)

18. The result is much better scaling for multi-part features.

19. Our cartogram is ready. This map shows concentration of the population on the eastern half of the US and a striking absence of population in states west of the Mississippi river.

20. We can make an improved visualization by creating an animation that slowly transforms the original rendering to the final size. Click the Temporal Control Panel button with the Clock icon in the Map Navigation Toolbar. Select the Animated temporal navigation button.

21. The default Animation range will be populated with a 24-hour window in the increment of 1-hour. This is fine for our use case as we will get 24-frames of animation. You can adjust this if you want a slower/faster animation. Right-click the us_states_with_population layer and select Properties.

22. Select Temporal tab and enable the Dynamic Temporal Control. This layer will be updated using expression so we don't need to configure it here. Just select Redraw Layer Only so that the layer is refreshed after each time step and rendered with the updated values from the expression.

23. Let's update our Geometry Generator expression to use the animation time steps and gradually scale the geometry. Click the Open the layer styling panel button in the Layers panel. Select Geometry Generator followed by the Expression Builder button.

24. Here we want to start with a scale factor of 1 and end up with the final scale factor value for the feature in the field scale_factor. We use the scale_linear() function which takes the time of the current time-step and calculates the scale factor using the start and end times. Update the expression as shown below and click OK.

collect_geometries(
    array_foreach(generate_series(1, @geometry_part_count),
            scale(geometry_n(@geometry,@element),
                    scale_linear(
                            epoch(@map_start_time),
                            epoch(@animation_start_time),
                            epoch(@animation_end_time),
                            1,
                            "scale_factor"),
                    scale_linear(
                            epoch(@map_start_time),
                            epoch(@animation_start_time),
                            epoch(@animation_end_time),
                            1,
                            "scale_factor"),
pole_of_inaccessibility(geometry_n(@geometry,@element), 100)
            )
    )
)

25. Back in the Temporal Controller panel, click the Play button to see the animation. You should see the shape of each polygon gradually scaled after each frame.

26. Once you are happy with the configuration, we can export the animation. Click the Export Animation button.

27. In the Export Map Animation dialog, click ... next to Output directory and browse to any folder on your computer. Keep all the other options to their default value and click Save.

28. The individual frames of the animation will be exported as images. We can create a video or animated GIF from these frames. I recommend using the website ezgif.com which allows you to create GIFs from individual images easily. Visit Ezgif Animated GIF Maker. Browse to the exported animation frames and click Upload files!.

29. Configure the GIF options by setting the Delay time to 5. Check crossfade frames effect and set the Fader delay and Fader count to 2. Click Make a GIF!/

30. Click the save button to download the animation as a GIF file.

---

If you want to give feedback or share your experience with this tutorial, please comment below. (requires GitHub account)