Rui Ying

2021-11-19

Introduction

This post is to simply introduce how to make a map in R environment, or more officially, geographic/spatial visualisation. However, I just found numerous methods to do this in R, depending on which data type you have. So this post will not introduce the basic of R and GIS. After all, I am also not an expert in this.

To make a map, at least two types of data are needed:

Geographic Information System (GIS) data for position
Research data per se for value

Sometimes these two are combined to a unified data (e.g., population data as an attribute in a shapefile). However, in most cases for those who are not studying geography (like me), we have to join the research data (in csv, netCDF) to GIS fields based on shared “key/id column”. The relevant GIS data probably can be found in third-party packages such as maps or mapdata, but for more demands like custom resolution, one have to download in other sources.

CSV file with lat and long columns

There are two methods to plot a csv data frame on a map:

plot layers using ggplot2
transform to sf

Because the latter one is same to sections below, so I will introduce ggplot2, which is very (but not limitedly) proficient at analysing data frames.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
library(ggplot2)
library(maps) # outline data for continents,countries
library(mapdata) # extra map datasets
library(mapproj) # map projection
#library(maptools)

# a quickest method to plot a global map
# map("world")

# Here I create a simple dataframe, which can be replaced to your data by read.csv()
toy_df <- data.frame(city = c("Bristol", "Shangrao"),
                     long = c(51.45, 28.45),
                     lat = c(-2.59, 117.94))

world <- map_data("world") # a ggplot2 function to get a world data.frame,see also borders

# notice x is lat and y is long, which differs normal impression
ggplot(world) +
  geom_polygon(aes(x = long, y = lat, group=group)) + 
  geom_label(data= toy_df, aes(lat, long, label=city)) +
  coord_map("mollweide", xlim=c(-180,180)) + # change the projection
  theme_minimal()

Actuallyggfortify provides the autoplot function, which is a quicker wrapper for ggplot2. So we can do the same by using less command.

1
2
3
4
5
library(ggfortify)
world <-  map('world', plot = FALSE, fill = TRUE)
autoplot(world, geom = 'polygon') +
  geom_point(data = toy_df, aes(lat, long),
             color="red", size=2)

There is another geom type: geom_map in ggplot. The difference between geom_map and geom_polygon is that the latter just plot the polygon (or position) and does not care about the value, but the geom_map contains both position and value information. Therefore, there must be a id columns shared by both polygon and value data frames, but no need to add lat/long in aes().

An example

1
2
3
4
5
6
7
8
fr_map <- map_data("france")
fr_data <- data.frame(region = unique(fr_map$region), 
                      value     = rnorm(length(unique(fr_map$region))))

ggplot(fr_data, aes(map_id = region)) +
  geom_map(aes(fill = value), map=fr_map) +
  expand_limits(x = fr_map$long, y = fr_map$lat) +
  coord_fixed()

sf object + geom_sf

shp file can be read in R through sf or rgal, and then plot using ggplot .

1
2
#from https://r-spatial.org/r/2018/10/25/ggplot2-sf.htm
library(sf)

1
## Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE

1
2
3
4
5
6
library(rnaturalearth) #another public map data package
#library(rnaturalearthdata)

it_data =st_as_sf(map("italy", plot = FALSE, fill = TRUE)) #polygons
it_sf <- ne_countries(scale = "medium", country="italy",returnclass = "sf") #an outline
ggplot(it_sf)+geom_sf()+geom_sf(data = it_data, fill = NA) 

1
2
#and you can add more geom_sf() based on data in your  hand
#plot(uk_sf$geometry) #uk_sf$geometry equals st_geomety(uk_sf)

sf object + tmap

tmap has a similar layer-based syntax to plot map, but should be more professional and specilised than ggplot2.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
#from https://geocompr.robinlovelace.net/adv-map.html
library(spData)
#library(spDataLarge)
library(tmap)
library(sf)

world <- spData::world
toy_sf <- st_as_sf(toy_df, coords = c("lat", "long"), crs=4326) #WGS84

tm_shape(world)+tm_borders()+tm_fill()+
  tm_shape(toy_sf) + tm_dots(size=.5)

ggmap

ggmap is another friendly package for people from other research background, it provides a function get_map() to download some map data from google, open street map and so on.

1
2
3
4
5
6
7
8
9
library(ggmap)

#downloading raster data by get_map(). REQUIRE GOOGLE API
gmap <- get_stamenmap(bbox = c(left = -10, bottom = 50, right = 5, top = 60),
                      zoom = 6,
                      maptype='watercolor',
                      messaging = FALSE)

ggmap(gmap) + theme_void()

raster + levelplot/geom_raster

netCDF data can be transformed to raster format.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
library(ncdf4)
library(terra)
library(rasterVis)

# an example file from UCAR
# https://www.unidata.ucar.edu/software/netcdf/examples/files.html
nc <- nc_open("tos_O1_2001-2002.nc")
#names(nc$var) # print all variables
sst_K <- ncvar_get(nc, varid="tos") #in Kelvin
sst <- sst_K - 273.15 #in Celsius
#dim(sst) # check dimension
#image(sst[,,1])

sst_raster <- flip(rast(t(sst[,,1])))
plot(sst_raster)

We can also use ggplot to do the plotting. But notice that the array from netCDF is not containing CRS data, that’s why we have x-y axis with 0-1 range.

1
2
3
4
5
6
library(ggplot2)

#to use geom_raster/geom_tile, we should transfer raster to a data.frame
sst_df <- sst_raster |> terra::as.data.frame(xy=TRUE) #pipe requires R > 4.1, otherwise use %>% in tidyverse
names(sst_df)[3] <- "sst"
ggplot(sst_df) + geom_raster(aes(x=x, y=y, fill=sst)) + theme_light() + scale_fill_viridis_c()

1
2
# what's the difference between |> and %>%?
#https://stackoverflow.com/questions/67633022/what-are-the-differences-between-rs-new-native-pipe-and-the-magrittr-pipe

As for the difference between geom_tile, geom_raster, let’s see the offical answer:

geom_rect() and geom_tile() do the same thing, but are parameterised differently: geom_rect() uses the locations of the four corners (xmin, xmax, ymin and ymax), while geom_tile() uses the center of the tile and its size (x, y, width, height). geom_raster() is a high performance special case for when all the tiles are the same size.

More resources

rworldmap for extra intersting data
cowplot for subplot added on ggplot objects
gganimate for animated plot
leaflet for interactive map

Cartography in R