A quick go at producing equal area cartograms based on English public health data from Public Health England’s fingertips website. The hexmapr package uses the Hungarian algorithm to assign hexagons (or rectangles) to original shpaefile polygons and can be found over on github: https://github.com/sassalley/hexmapr/.

Some people will absolutely hate this style of map (if you can truly call it a map), personally I think they’re very appealing but obviously need to be used carefully to ensure incorrect judgments and assumptions aren’t made.

0.1 Load Libraries

library("tidyverse"); theme_set(theme_bw())
library("broom") # to put shapefile into a dataframe format
library("devtools") # to enable installing from github
library("fingertipsR") # to access PHE fingertips data
#install_github("sassalley/hexmapr") # install hexmapr if not installed already
library("hexmapr")
library("viridis") # fancy colour palettes
library("gridExtra") # to arrange multiple plots

0.2 Load Data

I only look a the most recent data for a single region, the South East (Thames Valley represent). For now I will filter to males only as allowing for a merge which doubles the polygon count of a shapefile causes a big crash in hexmapr. We download the data at local authority district level and only look at one indicator, life expectancy at birth.

full_data <- fingertips_data(IndicatorID=90366,
                        AreaTypeID = 101,
                        stringsAsFactors=F)


# Filter!
data <- filter(full_data,
                  TimeperiodSortable==20130000 &
                  ParentName=="South East region" &
                  Sex=="Male")

0.3 Load Shapefile

Shapefile is obtained from the ONS geoportal. Contains OS data © Crown copyright and database right 2017

# Download shapefile from here
# http://geoportal.statistics.gov.uk/datasets/local-authority-districts-december-2016-generalised-clipped-boundaries-in-great-britain
# and load using hexmapr's read_polygons function

file_path <- "shp\\LAD_2016_Generalised_Clipped\\Local_Authority_Districts_December_2016_Generalised_Clipped_Boundaries_in_Great_Britain.shp"
LAD_shp <- read_polygons(file_path)
## Reading layer `Local_Authority_Districts_December_2016_Generalised_Clipped_Boundaries_in_Great_Britain' from data source `F:\Statistics\Hexmapr_1\shp\LAD_2016_Generalised_Clipped\Local_Authority_Districts_December_2016_Generalised_Clipped_Boundaries_in_Great_Britain.shp' using driver `ESRI Shapefile'
## Simple feature collection with 380 features and 10 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 5513.183 ymin: 5342.376 xmax: 655646.4 ymax: 1220304
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs
plot(LAD_shp)

Now we can merge the fingertips data to the shapefile.

LAD_shp_data <- sp::merge(LAD_shp,
                      data,
                      by.x="lad16cd",
                      by.y="AreaCode",
                      duplicateGeoms=T,
                      all.x=F)

0.4 Tidy Shapefile and Plot Map

Before we can plot the first map, we need to convert the shapefile into a dataframe object using a custom function which calls the tidy function (this process used to be called fortifying).

# define a function to handle tidying (which used to be called fortifying) and then joinng the data items back in
clean <- function(shape){
                          shape@data$id = rownames(shape@data)
                          shape.points = tidy(shape, region="id")
                          shape.df = inner_join(shape.points, shape@data, by="id")
                        }


# Apply function to shapefile
LAD_shp_data_tidy <- clean(LAD_shp_data)

# plot the map
(plot_1 <- ggplot(LAD_shp_data_tidy, aes(long.x,
                              lat.x,
                              fill=Value,
                              group=group)) +
  geom_polygon(col="white") +
  scale_fill_viridis(option="magma") +
  coord_equal() + theme_void())

0.5 Creating the Grid

Using different seed numbers, we can try different layouts of the hexagons which will be used in our final map.

LAD_shp_details <- get_shape_details(LAD_shp_data)

par(mfrow=c(3,3), mar = c(0,0,2,0))

for (i in 1:9){
new_cells <-  calculate_cell_size(LAD_shp_data, LAD_shp_details,0.03, 'hexagonal', i)
plot(new_cells[[2]], main = paste("Seed",i, sep=" "))
}

Based on this, three seems to be the best as it preserves the Isle of Wight’s location.

0.6 Create Final Grid

new_cells_hex <-  calculate_cell_size(LAD_shp_data, LAD_shp_details,0.03, 'hexagonal', 3)

# assigns data from the old shapefile to the new hexagonal polygons
resulthex <- assign_polygons(LAD_shp_data,new_cells_hex)

resulthex <- clean(resulthex)

0.7 Create Comparison Plot

# generate plot
plot_2 <- ggplot(resulthex, aes(long.x,
                              lat.x,
                              fill=Value,
                              group=group)) +
  geom_polygon(col="white") +
  scale_fill_viridis(option="magma") +
  coord_equal() + theme_void()

# arrange plot with comparison to the original
grid.arrange(plot_1, plot_2, nrow=1, ncol=2)

0.8 Create Large Final Plot

ggplot(resulthex, aes(long.x,
                              lat.x,
                              fill=Value,
                              group=group)) +
  geom_polygon(col="white") +
  geom_text(aes(V1,V2, label = substr(lad16nm,1,4)), size=3,color = "white") +
  scale_fill_viridis(option="magma",
                     begin = 0, end = 0.9,
                     name = "Male Life Expectancy at Birth") +
  coord_equal() + theme_void()+theme(legend.position="bottom") +
  guides(fill = guide_legend(title.position = "top")) +
  labs(title="Male Life Expectancy in the South East Region of England")