1 Introduction

There are now many ways of obtaining maps of the world in R in vector format. These can be useful if you want to analyse some statistics at a country, or administrative level.

library(tmap)
library(raster)
library(dismo)
library(rayshader)
library(elevatr)
library(rnaturalearth)
library(mapedit)
library(mapview)
library(dplyr)
library(tmaptools)
library(cartogram)
library(sf)
library(rmapshaper)
library(geojsonio)
library(sp)

2 The Natural Earth layers

Natural Earth has some useful layers that can be loaded directly from the web into R with a line of code. Countries and adminstrative districts (i.e counties and states) together with the standard letter codes.

countries<-ne_countries(scale = 110, type = "countries", continent = NULL,
  country = NULL, geounit = NULL, sovereignty = NULL,
  returnclass = c("sf"))

admin<-ne_states(country = NULL, geounit = NULL, iso_a2 = NULL, spdf = NULL,
  returnclass = c("sf"))

mapview(countries)

3 Subsetting to a local area

admin %>% filter(name=="Dorset") -> dor
countries %>% filter(name=="United Kingdom") -> UK
mapview(UK) + mapview(dor)

4 Digitising

It is clearly not possible to digitise in an interactive makdown document. However running the following code will allow you to add points, polygons or lines in an interactive session.

admin %>% filter(name=="Dorset") -> dor

This code, which is not run in the markdown, will pop open a mapview in the console or in RStudio. Capture a polygon and save.

mapview(dor) %>% editMap() ->d     ## Run this in real time in a session to capture digitised features.
write_sf(d$drawn,"poly",driver= "Esri Shapefile") ## Write the digitised features to a file for later use.

5 Corfe

corfe<-read_sf("poly")
mapview(corfe)

6 Getting some elevation data

The getData function in the raster package used to return SRTM raster layers, but seems to be broken at the time of writing. However the get_elev_raster function in the package elevatr works. So if you go through the steps of finding an adminstrative district, zooming in and selecting a local region then running the code you will have an elevation model for your area.

elev1 <- get_elev_raster(corfe, z = 14,clip="bbox")
contours<-rasterToContour(elev1,nlevels=20)
mapview(contours,legend=TRUE)

7 All code together

The code chunk below is not run here, but just cutting and pasting the code in this chunk into a new document will do all the above operations, without saving the digitised layer locally.

library(tmap)
library(raster)
library(elevatr)
library(rnaturalearth)
library(mapedit)
library(mapview)
library(dplyr)


countries<-ne_countries(scale = 110, type = "countries", continent = NULL,
  country = NULL, geounit = NULL, sovereignty = NULL,
  returnclass = c("sf"))

admin<-ne_states(country = NULL, geounit = NULL, iso_a2 = NULL, spdf = NULL,
  returnclass = c("sf"))

admin %>% filter(name=="Dorset") -> dor
countries %>% filter(name=="United Kingdom") -> UK

mapview(dor) %>% editMap() ->d  

poly<-d$drawn
elev1 <- get_elev_raster(poly, z = 14,clip="bbox")

8 Ray shader visualisation

The Rayshader package produces dynamic 3gl plots, but can also produce static hillshade plots.

library(rayshader)

elevation_matrix<-t(as.matrix(elev1))
elevation_matrix %>%
  sphere_shade(texture = "imhof2") %>%
  
  #add_shadow(ambient_shade(elevation_matrix)) %>%
  plot_map()
## 
Calculating Surface Normal [-----------------------------------] ETA:  0s
Calculating Surface Normal [===================================] ETA:  0s

9 San Cristóbal

admin %>% filter(name=="Chiapas") -> chis
mapview(chis) %>% editMap() ->d     ## Uncomment these lines and run in real time
write_sf(d$drawn,"sclc",driver= "Esri Shapefile") ## Write the digitised polygon.
sclc<-read_sf("sclc")
mapview(sclc)
library(elevatr)
elev2 <- get_elev_raster(sclc, z = 14,clip="bbox")
contours<-rasterToContour(elev2,nlevels=100)
mapview(contours)
elevation_matrix<-t(as.matrix(elev2))
elevation_matrix %>%
  sphere_shade(texture = "imhof2") %>%
  
  #add_shadow(ambient_shade(elevation_matrix)) %>%
  plot_map()

10 Using Tmap

The package Tmap makes chroropleths relatively easily.

library(tmap)
data(World, metro)
World$color <- map_coloring(World, palette="Pastel2")

qtm(World, fill="HPI", fill.n = 9, fill.palette = "div",
    fill.title = "Happy Planet Index", fill.id = "name", 
    style = "classic", format = "World")

11 Conclusion

Once the elevation data is loaded in R there are many other things that can be done with it. This just provides a simple way of getting data for a specific region by choosing it using mapedit.