Honest maps: classification, missingness and distortion

mapdf <- attach_geometry(snap, geometry = "polygon")

A world choropleth makes four claims before it says anything about your data: that the classes it drew are the natural ones, that grey means nothing rather than zero, that a rate over eleven thousand people is worth as much attention as one over a billion, and that the shapes on screen are the shapes on the ground. All four are usually false. This vignette is the tour of what countryatlas gives you for each.

1. The classification is doing the talking

Brewer & Pickle (2002) ran 56 subjects over nine series of mortality maps and found quantiles among the most accurately read classifications, with natural breaks (Jenks) below 70% as accurate. That is the reverse of the common GIS default, and it matters because the choice is not cosmetic: it decides what the reader concludes.

classify_compare() draws the same data under several methods at once.

cmp <- classify_compare(mapdf, gdp_per_capita,
                        methods = c("quantile", "jenks", "equal", "pretty"),
                        ncol = 2)
cmp

The picture is only half of it. The counts are attached to the plot:

attr(cmp, "countryatlas_classification")
#> # A tibble: 20 × 7
#>    method   class              n   share   gvf   tai max_class_share
#>    <chr>    <chr>          <int>   <dbl> <dbl> <dbl>           <dbl>
#>  1 quantile 269 to 1.68K      40 0.201   0.626 0.672           0.201
#>  2 quantile 1.68K to 4.65K    40 0.201   0.626 0.672           0.201
#>  3 quantile 4.65K to 10.3K    39 0.196   0.626 0.672           0.201
#>  4 quantile 10.3K to 30.1K    40 0.201   0.626 0.672           0.201
#>  5 quantile 30.1K to 247K     40 0.201   0.626 0.672           0.201
#>  6 jenks    269 to 13.1K     129 0.648   0.958 0.756           0.648
#>  7 jenks    13.1K to 34.8K    38 0.191   0.958 0.756           0.648
#>  8 jenks    34.8K to 68.1K    25 0.126   0.958 0.756           0.648
#>  9 jenks    68.1K to 117K      6 0.0302  0.958 0.756           0.648
#> 10 jenks    117K to 247K       1 0.00503 0.958 0.756           0.648
#> 11 equal    269 to 49.6K     181 0.910   0.787 0.433           0.910
#> 12 equal    49.6K to 99K      15 0.0754  0.787 0.433           0.910
#> 13 equal    99K to 148K        2 0.0101  0.787 0.433           0.910
#> 14 equal    148K to 198K       0 0       0.787 0.433           0.910
#> 15 equal    198K to 247K       1 0.00503 0.787 0.433           0.910
#> 16 pretty   0 to 50K         181 0.910   0.787 0.433           0.910
#> 17 pretty   50K to 100K       15 0.0754  0.787 0.433           0.910
#> 18 pretty   100K to 150K       2 0.0101  0.787 0.433           0.910
#> 19 pretty   150K to 200K       0 0       0.787 0.433           0.910
#> 20 pretty   200K to 250K       1 0.00503 0.787 0.433           0.910

Equal-interval and pretty breaks put over 90% of countries into a single class, because GDP per capita is strongly right-skewed and the top of the range is one country. A map like that is technically correct and communicates nothing. Quantiles put roughly 40 countries in each class.

You do not need the comparison to get the report: any classified world_map() will produce it. (A continuous colourbar has no classes, so asking there returns nothing and says why.)

p <- world_map(mapdf, gdp_per_capita, style = "quantile",
               classification_report = TRUE)
attr(p, "countryatlas_classification")
#> # A tibble: 5 × 6
#>   method   class              n share   gvf   tai
#>   <chr>    <chr>          <int> <dbl> <dbl> <dbl>
#> 1 quantile 269 to 1.68K      40 0.201 0.626 0.672
#> 2 quantile 1.68K to 4.65K    40 0.201 0.626 0.672
#> 3 quantile 4.65K to 10.3K    39 0.196 0.626 0.672
#> 4 quantile 10.3K to 30.1K    40 0.201 0.626 0.672
#> 5 quantile 30.1K to 247K     40 0.201 0.626 0.672

Jenks still earns its place: on a strongly clustered distribution, quantiles will split a natural group across two colours where Jenks keeps it together. The point is to look, not to take the default.

2. Grey is not a value

The default no-data grey reads as “low” to a lot of readers, which is precisely the wrong inference. na_style gives you three alternatives, and "hatched" (via the optional ggpattern) is the one that survives both colour-blindness and a black-and-white printer.

world_map(mapdf, co2_per_capita, style = "quantile",
          na_style = "hatched", footnote = "auto")

footnote = "auto" writes the coverage line into the caption, so the map cannot quietly overstate what it covers.

When the missingness is the story, map it directly:

coverage_map(mapdf, co2_per_capita)

audit_coverage() is the same question as a table, and is the better tool when you want to act on the answer rather than look at it.

3. Small denominators shout

Any per-capita or per-100k figure computed over a tiny population is mostly noise, and on a choropleth it gets exactly as much ink as a figure computed over a billion people. There are two answers. The cartogram distorts geometry until area matches the denominator; value-by-alpha (Roth, Woodruff & Johnson 2010) leaves the geometry alone and spends opacity instead.

value_by_alpha_map(mapdf, gdp_per_capita, population)

Countries fade toward the background in proportion to how little population stands behind their number. Compare with the cartogram answer to the same problem, cartogram_map() / dorling_map(), in Beyond the choropleth: the trade-off is that a cartogram makes the weighting unmissable but costs you the recognisable world.

4. The projection is doing the talking too

Every flat world map distorts something. projection_info() says what each of the thirteen preserves:

projection_info()[, c("projection", "property", "equal_area", "conformal")]
#> # A tibble: 13 × 4
#>    projection           property    equal_area conformal
#>    <chr>                <chr>       <lgl>      <lgl>    
#>  1 equal_earth          equal-area  TRUE       FALSE    
#>  2 robinson             compromise  FALSE      FALSE    
#>  3 mollweide            equal-area  TRUE       FALSE    
#>  4 natural_earth        compromise  FALSE      FALSE    
#>  5 plate_carree         equidistant FALSE      FALSE    
#>  6 mercator             conformal   FALSE      TRUE     
#>  7 winkel_tripel        compromise  FALSE      FALSE    
#>  8 eckert4              equal-area  TRUE       FALSE    
#>  9 gall_peters          equal-area  TRUE       FALSE    
#> 10 orthographic         perspective FALSE      FALSE    
#> 11 azimuthal_equal_area equal-area  TRUE       FALSE    
#> 12 north_polar          equal-area  TRUE       FALSE    
#> 13 south_polar          equal-area  TRUE       FALSE

For a choropleth the honest choice is equal-area, because the eye reads coloured area as quantity, and a projection that inflates Greenland makes Greenland’s value look more important than it is. Equal Earth is the package default and the recommendation (Šavrič, Patterson & Jenny 2019).

subset(projection_info(), equal_area)$projection
#> [1] "equal_earth"          "mollweide"            "eckert4"             
#> [4] "gall_peters"          "azimuthal_equal_area" "north_polar"         
#> [7] "south_polar"

Tissot’s indicatrix makes the cost visible. Each circle has the same radius on the ground; whatever the projection does to them, it is doing to your data.

tissot_map("mercator")

tissot_map("equal_earth")

Mercator keeps every circle round (it is conformal, so local shapes are right) and grows them without limit toward the poles. Equal Earth keeps every circle’s area and shears the shapes instead. Neither is wrong; they are answers to different questions, and only one of them belongs under a choropleth.

To see it on your own data, vary the CRS and hold everything else fixed:

attach_geometry(snap, geometry = "sf") |>
  projection_compare(gdp_per_capita, style = "quantile", labeller = "property")

5. Say what the map is

Everything above is a decision, and a published map should carry its decisions. map_provenance() reads them back off the plot.

world_map(mapdf, gdp_per_capita, style = "quantile", n_bins = 5,
          na_style = "hatched", footnote = "auto") |>
  map_provenance()
#> 
#> ── countryatlas map provenance
#> package: countryatlas 4.0.0 (snapshot 2024)
#> fill: gdp_per_capita
#> geometry: polygon backend, equal_earth
#> classification: quantile, 5 bins
#> missing data: hatched
#> coverage: 199 countries shown, 39 missing
#> breaks: 268.7 | 1684 | 4655 | 10250 | 30130 | 247200
#> data: World Bank WDI NY.GDP.PCAP.KD (constant 2015 US$), release 2026-07,
#> fetched 2026-10-02

Every field there was already known when the plot was built; the only new thing is that you can read it. Paired with footnote = "auto" on the plot itself and the classification report, that is most of a methods note.

Finally, citation("countryatlas") produces the package citation and the sources it reconciles: countrycode, the World Bank, Natural Earth, and the papers behind the methods used here. Citing the join layer without the data would be the last dishonest thing a map could do.

6. Where the data comes from, and when

Two more ways a country map goes quietly wrong, both added in 3.0.0.

The borders are not the borders. A 1950 map drawn on 2024 boundaries is simply a different world. attach_geometry(year = ) and [historical_geometry()] draw the real ones, from CShapes – including the colonies, without which most of Africa and Asia is absent:

attach_geometry(snap[, c("iso3c", "gdp_per_capita")], year = 1950) |>
  world_map(gdp_per_capita, style = "quantile",
            title = "1950 borders, 1950 world")

Note what this costs: ISO 3166 was published in 1974 and never covered colonies, so historical geometry is keyed on Gleditsch-Ward codes and iso3c is NA for every entity that never had one. country_join(key = "gwn") is the join that works before 1970.

Membership is a function of time. A snapshot silently misstates any panel that spans an accession:

c(`2016` = in_group("United Kingdom", "EU", as_of = 2016),
  `2021` = in_group("United Kingdom", "EU", as_of = 2021))
#>  2016  2021 
#>  TRUE FALSE

7. Islands are not missing at random

Land-border contiguity, morans_i()’s default weights before 4.0.0, gives an island no neighbours at all. On the bundled snapshot that removes a quarter of the countries with data: Japan, Australia, Madagascar, New Zealand, the Philippines, Cuba, Sri Lanka, Iceland and every small island state. They are not a random quarter, and the default is now the five nearest neighbours, so every country with data and a centroid takes part.

Not every island goes, either: the United Kingdom keeps its land border with Ireland, and Indonesia keeps its borders with Malaysia, Papua New Guinea and Timor-Leste. Which is rather the point – you cannot tell from the finished map who dropped out of the statistic.

rbind(
  contiguity = morans_i(snap, gdp_per_capita, n_perm = 0,
                        weights = country_weights("contiguity"))[c("i", "n", "n_excluded")],
  knn = morans_i(snap, gdp_per_capita, n_perm = 0)[c("i", "n", "n_excluded")]
)
#> # A tibble: 2 × 3
#>       i     n n_excluded
#> * <dbl> <int>      <int>
#> 1 0.600   146         53
#> 2 0.452   199          0

Both numbers are defensible; only one of them is global. country_weights() also takes "distance", and "custom" – which is how an adjacency that is not geographic at all (trade volume, migration, shared language) goes through the same API.

References

Brewer, C. A. & Pickle, L. (2002). Evaluation of methods for classifying epidemiological data on choropleth maps in series. Annals of the Association of American Geographers 92(4), 662–681.

Roth, R. E., Woodruff, A. W. & Johnson, Z. F. (2010). Value-by-alpha maps: an alternative technique to the cartogram. The Cartographic Journal 47(2), 130–140.

Šavrič, B., Patterson, T. & Jenny, B. (2019). The Equal Earth map projection. International Journal of Geographical Information Science 33(3), 454–465.