MÉMOVISUELDAGRONOMIE

Comparison of herds and human density

This map compares livestock and human density using bivariate choropleths.

A bivariate choropleth map is a type of statistical thematic map in which two quantitative variables are represented, blending two color gradients

1. Load data

The data we’re going to use here is available in the {frex} package I’ve created.

The goal of this package is to provide several layers of information for metropolitan France, particularly useful for analyzing agricultural systems. The basis of this package is a gridded map of France in hexagons of about 450 km2. For each hexagon, different information can be added using the package function, such as the surface area occupied by different types of crop, the nature of the soil or climatic data.

For the sake of simplicity, we won’t install the whole package, but only load the data we need, starting with the basemap of “gridded” France.

# Load packages
library(tidyverse)
library(sf)

# Load basemap
hex<-read_sf('https://github.com/BjnNowak/frex_db/raw/main/map/hex_grid.gpkg')

# Make plot
ggplot(hex)+
  geom_sf()

We will now load the data tables to fill the basemap (one table with the human population density, an other with the farm animal population density).

They may be loaded as shown below:

# Load data
# Human population (per km2)
human_pop<-read_csv('https://raw.githubusercontent.com/BjnNowak/frex_db/main/data/human/human_pop_density.csv')
# Farm animal density (per km2)
herd<-read_csv('https://raw.githubusercontent.com/BjnNowak/frex_db/main/data/herds/herds_density_municipality.csv')

If you look at the hex object, you will see that each hexagon has a unique identifiant (hex_id column), that is also present in the data files: hence we can join the data tables to the basemap.

# Join data to map
data<-hex%>%
  left_join(human_pop)%>%
  left_join(herd)

Before going further, we will make a simple choropleth showing milking cows density:

ggplot(data, aes(fill=milk_cow_km2))+
  geom_sf()

2. Prepare the data for bivariate choropleth

We will now see how we can reproduce the bivariate choropleth map comparing human and cow density

First, fe will create a single column for cow density (combining milk and meat cows).

# Prepare the data
data<-data%>%
  mutate(cow_km2=milk_cow_km2+meat_cow_km2)

One very important thing for the final rendering of the bivariate choropleth is to set the appropriate thresholds for the two variables:

In our case, it seems appropriate to set the same thresholds for both variables (but it’s not always like that). To obtain a roughly balanced number of hexagons with each color in the palette, the most conventional way is to set the thresholds using tiertiles

quantile(data$human_km2, na.rm=TRUE, probs=c(0.33,0.66))
##      33%      66% 
## 31.81796 80.92791
quantile(data$cow_km2, na.rm=TRUE, probs=c(0.33,0.66))
##       33%       66% 
##  4.374322 20.839668

Here, by taking into account both human and cow tiertiles, thresholds can be set at c(25,50)

We will now assign classes to the hexagons according to these thresholds.

To do so, we will consider that the human population density is the variable A (variable B for cows).

We will start with the human population density (variable A):

# Set thresholds
t1<-25
t2<-50

# Set classes for human pop
data<-data%>% 
  replace(is.na(.), 0)%>%
  mutate(cl_A=case_when(
    human_km2<t1~'A1',
    human_km2<t2~'A2',
    TRUE~'A3'
  )
)

Now the same for cow density (variable B):

# Set classes for cow density
data<-data%>% 
  mutate(cl_B=case_when(
    cow_km2<t1~'B1',
    cow_km2<t2~'B2',
    TRUE~'B3'
  ))

… and finally merging both variables in one cl column:

data<-data%>%
  mutate(cl=paste0(cl_A,cl_B))

Note: the preceding workflow may be but in a function for easier reuse (eventually renaming the columns to “VarA” and “VarB” and defining them as arguments, with the thresholds).

3. Plot bivariate choropleth

We are now ready to fill the hex based on the newly created column (cl):

p1<-ggplot(data, aes(fill=cl))+
  geom_sf()
p1

but we need a better color palette !

Bivariate color palettes require two colors that blend well together. So there are a relatively limited number of ‘classic’ bivariate color palettes. Some are available in the {pals} package:

Source Jakub Nowosad

With {pals}, a color palette may be created as follows:

library(pals)

pal<-stevens.pinkblue(9)

# Assigning color palette to our combinations:
pal_bivar<-c(
  'A1B1'=pal[1],
  'A1B2'=pal[2],
  'A1B3'=pal[3],
  'A2B1'=pal[4],
  'A2B2'=pal[5],
  'A2B3'=pal[6],
  'A3B1'=pal[7],
  'A3B2'=pal[8],
  'A3B3'=pal[9]
)

Adding color palette to ggplot():

p2<-p1+
  scale_fill_manual(values=pal_bivar)+
  # Hide legend
  guides(fill='none')+
  theme_void()


p2

Looks better right?

Note that you may also create custom color palette from this Observable app from BVG software, or use an online color picker to extract hex codes for another online palette.

4. Create legend

We will now add the color legend to the plot.

To do so, we will start by creating a tibble with various combinations of classes for both variables (human population and animal density):

tib<-tibble(
  var_A = rep(c("A1","A2","A3"),3),
  var_B = c(rep("B1",3),rep("B2",3),rep("B3",3)),
  value = paste0(var_A,var_B)
)

This tibble may be plotted on a grid as follows:

leg_tile<-ggplot(data=tib,aes(x=var_A,y=var_B,fill=value))+
  geom_tile()+
  scale_fill_manual(values=pal)+
  guides(fill='none')+
  coord_fixed()+
  labs(x="Human density",y="Animal density")+
  theme_minimal()+
  theme(axis.text=element_blank())

leg_tile

…also a slightly different version with points:

leg_pts<-ggplot(data=tib,aes(x=var_A,y=var_B,fill=value))+
  geom_point(pch=21,size=60,color="grey90")+
  scale_fill_manual(values=pal)+
  guides(fill='none')+
  coord_fixed()+
  labs(x="Human density",y="Animal density")+
  theme_minimal()+
  theme(axis.text=element_blank())

leg_pts

We may combine several ggplot() with {patchwork} or {cowplot}. An example with {patchwork}:

library(patchwork)

leg_tile+p2

Slightly different version with different sizes for plots:

layout <- c(
  area(t = 8, l = 1, b = 10, r = 2),
  area(t = 1, l = 1, b = 10, r = 10)
)

leg_tile+p2+ 
  plot_layout(design = layout)

Using the same workflow, or the {biscale} package, you may now try to create another bivariate choropleth map with the available data (comparing for example human and sheep density, sheep and goats density…) !