For problems with the code, please open an issue on GitHub here.

Contact details: GitHub / seethedatablog / LinkedIn / Twitter

1 Summary

The main steps for this map consist in getting country polygons data, cropping and shifting the geographical coordinates for a Pacific view, labeling some entities (here some small landmasses in Pacific) and generating boundaries between them using the Voronoi tessellation method.

The Rmarkdown file with the R code used to generate this html report can be found here.

2 The R environment

2.1 Load packages

library(rgeos)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
library(dplyr)
library(ggplot2)
library(ggrepel)

2.2 R Session

sessionInfo()

## R version 4.0.0 (2020-04-24)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18362)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ggrepel_0.8.2           ggplot2_3.3.0           dplyr_0.8.5            
## [4] rnaturalearthdata_0.1.0 rnaturalearth_0.1.0     sf_0.9-3               
## [7] rgeos_0.5-2             sp_1.4-1               
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.4.6       compiler_4.0.0     pillar_1.4.4       class_7.3-16      
##  [5] tools_4.0.0        digest_0.6.25      gtable_0.3.0       evaluate_0.14     
##  [9] tibble_3.0.1       lifecycle_0.2.0    lattice_0.20-41    pkgconfig_2.0.3   
## [13] rlang_0.4.6        DBI_1.1.0          yaml_2.2.1         xfun_0.13         
## [17] e1071_1.7-3        withr_2.2.0        stringr_1.4.0      knitr_1.28        
## [21] vctrs_0.2.4        classInt_0.4-3     grid_4.0.0         tidyselect_1.0.0  
## [25] glue_1.4.0         R6_2.4.1           rmarkdown_2.1      purrr_0.3.4       
## [29] magrittr_1.5       scales_1.1.0       htmltools_0.4.0    ellipsis_0.3.1    
## [33] units_0.6-6        assertthat_0.2.1   colorspace_1.4-1   KernSmooth_2.23-16
## [37] stringi_1.4.6      munsell_0.5.0      crayon_1.3.4

2.3 Package versions

If you want to install a snapshot of the packages as they existed on CRAN at the date of the creation of this document, then can use the checkpoint package. The code below should guarantee the same versions of the packages as they existed on CRAN at the specified point in time.

library(checkpoint) # checkpoint, version 0.4.9
checkpoint(snapshotDate = "2020-05-16")
# Check that library path is set to ~/.checkpoint/2020-05-16/ ...
.libPaths()
grepl(pattern = "\\.checkpoint/2020-05-16/", x = .libPaths()[1]) # should be TRUE

The checkpoint functionality consists on scanning for package names in the scripts and text files of your project folder and its subfolder. It scans all R code (.R, .Rmd, and .Rpres files) for library() and require() statements. Then creates a local library into which it installs a copy of the packages required in the project as they existed on CRAN at the specified snapshot date. See details with ?checkpoint after you load library(checkpoint), or here

Warning - Installing older versions of the packages in the checkpoint local library, may take up some hundreds of MG in space.

3 Data preparation

3.1 Set ggplot2 theme

This is more convenient, because once set here, is valid for any generated ggplot in this report. Feel free to adapt and change your theme preferences.

theme_set(
  theme_bw() +
    theme(panel.background = element_rect(fill = "azure"),
          panel.grid.major = element_blank(),
          axis.title = element_blank(),
          axis.text = element_text(size = 8))
)

3.2 Natural Earth data

Get world spatial polygons from the rnaturalearth package and we’ll focus on the Pacific area.

world <- rnaturalearth::ne_countries(scale = 'medium', returnclass = "sp")
box_cut <- bbox2SP(n = 90, s = -90, w = -70, e = 120, proj4string = world@proj4string)
world_crop <- gDifference(world, box_cut)

The main steps consist in cutting out the red area and binding the two green ones with shifting the geographical coordinates for a Pacific view as illustrated below:

Next we will further crop the results from above to a focus area on the Pacific depicted in green below:

So, we take the world polygons obtained by cutting out the central red box, shift the geographical coordinates for a Pacific view and crop to focus on a more narrow area:

pacific_crop <- world_crop %>% 
  st_as_sf() %>% # change from sp to sf object/class
  st_shift_longitude() %>% 
  st_crop(c(xmin = st_bbox(.)[["xmin"]],
            xmax = st_bbox(.)[["xmax"]],
            ymin = -50,
            ymax = 30))

## although coordinates are longitude/latitude, st_intersection assumes that they are planar

ggplot() +
  geom_sf(data = pacific_crop)

For the labeling exercise we can use the rnaturalearthdata::tiny_countries50 dataset. This set of spatial points need the be shifted and cropped also.

tiny_countries <- rnaturalearthdata::tiny_countries50 %>% 
  st_as_sf() %>%
  st_shift_longitude() %>% 
  st_crop(c(xmin = 120, xmax = 250, ymin = -50, ymax = 30)) %>% 
  # Also adds the coordinates to be used for labeling with geom_text_repel
  bind_cols(st_coordinates(.) %>% as.data.frame())

## although coordinates are longitude/latitude, st_intersection assumes that they are planar

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

ggplot() +
  geom_sf(data = pacific_crop) +
  geom_sf(data = tiny_countries, size = 2, color = "red")

3.3 Voronoi tessellation

This will create boundaries between the land masses we decided to label. Cropping is needed here, because otherwise we get very long segments at the edges of the Voronoi grid.

Voronoi tessellation might not be your cup of tea, since the generated borders are not perfect, but they are a strong simplification of reality. You may prefer generating buffers around your locations, though there you might have problems with overlapping polygons.

voronoi_lines <- tiny_countries %>%
  st_as_sf() %>% 
  st_union() %>% # this seems to be needed, otherwise st_voronoi will fail
  st_voronoi(bOnlyEdges = TRUE) %>% 
  st_crop(st_bbox(pacific_crop))

## Warning in st_voronoi.sfc(., bOnlyEdges = TRUE): st_voronoi does not correctly
## triangulate longitude/latitude data

## Warning in st_is_longlat(x): bounding box has potentially an invalid value range
## for longlat data

## although coordinates are longitude/latitude, st_intersection assumes that they are planar

Note - ignore the warnings. Things do not have to be perfect here, we use the Voronoi grid for visual effects only.

4 The map

Put all layers together in a map

pacific_map_voronoi <- ggplot() +
  geom_sf(data = voronoi_lines, color = gray(level = 0.6), size = 0.5) +
  geom_sf(data = pacific_crop, size = 0.25) +
  geom_sf(data = tiny_countries, color = "black", size = 0.5) +
  # Here, coord_sf ensures that all layers are displayed within the same limits,
  # the extremes being dictated by the minimum and maximum of pacific_crop. Also
  # the expanding factor is turned off so that the content of the map can touch
  # the axes:
  coord_sf(xlim = c(120, 290), expand = FALSE) +
  geom_text_repel(data = tiny_countries,
                  aes(x = X, y = Y, label = name),
                  size = 2.5,
                  fontface = "bold",
                  segment.size = 0.25,
                  box.padding = 0.4,
                  min.segment.length = 0,
                  seed = 2020-5-16)
pacific_map_voronoi

If you need to further save the map as a png file, you can, of course, use ggsave():

ggsave(filename = "pacific_map_voronoi.png",
       plot = pacific_map_voronoi,
       width = 15,
       height = 10,
       units = "cm",
       dpi = 300)

Pacific centered map and Voronoi tessellation

by: Valentin Stefan

last update: 2020-05-23