Compare the speed of ggplot creation: geom_poly() VS geom_map()

microbenchmark-poly-vs-map.R

################################################################################
#                                                    
# ikashnitsky.github.io 2017-07-17
# geom_map VS geom_polygon: the speed of maps creation
# Ilya Kashnitsky, ilya.kashnitsky@gmail.com
#                                                  
################################################################################

# load required packages
library(tidyverse)      # data manipulation and viz
library(ggthemes)       # themes for ggplot2
library(viridis)        # the best color palette
library(rgdal)          # deal with shapefiles
library(microbenchmark) # measure the speed of executing
library(extrafont)      # nice font
myfont <- "Roboto Condensed"
library(eurostat)       # grab data
library(rgdal)          # deal with shapefiles
library(rgeos)
library(maptools)



################################################################################
# Download geodata
# Eurostat official shapefiles for regions
# http://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/administrative-units-statistical-units

# geodata will be stored in a directory "geodata"
ifelse(!dir.exists('geodata'),
       dir.create('geodata'),
       paste("Directory already exists"))

f <- tempfile()
download.file("http://ec.europa.eu/eurostat/cache/GISCO/geodatafiles/NUTS_2013_60M_SH.zip", destfile = f)
unzip(f, exdir = "geodata/.")
NUTS_raw <- readOGR("geodata/NUTS_2013_60M_SH/data/.", "NUTS_RG_60M_2013")

# colnames to lower case
names(NUTS_raw@data) <- tolower(names(NUTS_raw@data))

NUTS_raw@data %>% View

# let's have a look
plot(NUTS_raw)

# change coordinate system to LAEA Europe (EPSG:3035)
# check out https://espg.io
epsg3035 <- "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

NUTS <- spTransform(NUTS_raw, CRS(epsg3035))

# subset NUTS-3 regions
NUTS3 <- NUTS[NUTS$stat_levl_==3,]

# make the geodata ready for ggplot
gd <- fortify(NUTS3, region = "nuts_id") 


################################################################################
# Get stat data

# Find the needed dataset code 
# http://ec.europa.eu/eurostat/web/regions/data/database

# download fertility rates for countries
df <- get_eurostat("demo_r_frate3")

# if the automated download does not work, the data can be grabbed manually at
# http://ec.europa.eu/eurostat/estat-navtree-portlet-prod/BulkDownloadListing

tfr15 <- df %>% filter(geo %>% paste() %>% nchar() == 5, 
                       time == "2015-01-01")





################################################################################
# test speed with microbenchmark


# pre-join data for geom_poly
gdf <- left_join(gd, tfr15, by = c("id" = "geo"))


test_map <- microbenchmark(
        
        "geom_poly" = {
                left_join(gd, tfr15, by = c("id" = "geo")) %>% 
                ggplot()+
                coord_equal(ylim = c(1350000,5450000), 
                            xlim = c(2500000, 6600000), 
                            expand = c(0,0))+
                geom_polygon(aes(x = long, y = lat, group = group, 
                                 fill = values))+
                scale_fill_viridis()+
                theme_map(base_family = myfont)+
                theme(legend.position = c(1, 1),
                      legend.justification = c(1, 1))
        }
        ,
        
        "poly - pre-join" = {
                gdf %>% 
                ggplot()+
                coord_equal(ylim = c(1350000,5450000), 
                            xlim = c(2500000, 6600000), 
                            expand = c(0,0))+
                geom_polygon(aes(x = long, y = lat, group = group, 
                                 fill = values))+
                scale_fill_viridis()+
                theme_map(base_family = myfont)+
                theme(legend.position = c(1, 1),
                      legend.justification = c(1, 1))
        }
        ,
        
        "geom_map" = {
                tfr15 %>% 
                ggplot()+
                coord_equal(ylim = c(1350000,5450000), 
                            xlim = c(2500000, 6600000), 
                            expand = c(0,0))+
                geom_map(map = gd, 
                         aes(map_id = geo, fill = values))+
                scale_fill_viridis()+
                theme_map(base_family = myfont)+
                theme(legend.position = c(1, 1),
                      legend.justification = c(1, 1))
        }
        ,
        
        times = 1000
)


test_map
      

autoplot(test_map)+
        aes(fill = expr)+
        scale_fill_viridis(discrete = T)+
        theme_bw(base_size = 15, base_family = myfont)+
        theme(legend.position = "none",
              axis.text = element_text(size = 15))+
        labs(title = "The speed of creating a map with ggplot2")



################################################################################
# smaller map - Baltic countries

gd_balt <- fortify(NUTS, region = "nuts_id") %>% 
        filter(id %in% c("EE", "LT", "LV"))

df_balt <- df %>% filter(geo %in% c("EE", "LT", "LV"),
                         time == "2015-01-01")

# pre-join data for geom_poly
gdf_balt <- left_join(gd_balt, df_balt, by = c("id" = "geo"))


test_balt <- microbenchmark(
        
        "geom_poly" = {
                left_join(gd_balt, df_balt, by = c("id" = "geo")) %>%  
                ggplot()+
                coord_equal(xlim = c(49e5, 54e5), 
                            ylim = c(35e5, 42e5), 
                            expand = c(0,0))+
                geom_polygon(aes(x = long, y = lat, group = group, 
                                 fill = values),
                             color = "grey50")+
                scale_fill_viridis()+
                theme_void(base_family = myfont)+
                theme(legend.position = "right")
        }
        ,
        
        "poly - pre-join" = {
                gdf_balt %>% 
                ggplot()+
                coord_equal(xlim = c(49e5, 54e5), 
                            ylim = c(35e5, 42e5), 
                            expand = c(0,0))+
                geom_polygon(aes(x = long, y = lat, group = group, 
                                 fill = values),
                             color = "grey50")+
                scale_fill_viridis()+
                theme_void(base_family = myfont)+
                theme(legend.position = "right")
        }
        ,
        
        "geom_map" = {
                df_balt %>% 
                ggplot()+
                coord_equal(xlim = c(49e5, 54e5), 
                            ylim = c(35e5, 42e5), 
                            expand = c(0,0))+
                geom_map(map = gd_balt, 
                         aes(map_id = geo, fill = values),
                         color = "grey50")+
                scale_fill_viridis()+
                theme_void(base_family = myfont)+
                theme(legend.position = "right")
        }
        ,
        
        times = 1000
)


test_balt


autoplot(test_balt)+
        aes(fill = expr)+
        scale_fill_viridis(discrete = T)+
        theme_bw(base_size = 15, base_family = myfont)+
        theme(legend.position = "none",
              axis.text = element_text(size = 15))+
        labs(title = "The speed of creating a map with ggplot2")