Zooming in on maps with sf and ggplot2

zoom.R

# Source for blog post "Zooming in on maps with sf and ggplot2"
# URL: https://datascience.blog.wzb.eu/2019/04/30/zooming-in-on-maps-with-sf-and-ggplot2/
#
# Markus Konrad <markus.konrad@wzb.eu>
# Wissenschaftszentrum Berlin für Sozialforschung
# April 30, 2019
# 

#### world map ####

library(ggplot2)
library(sf)
library(rnaturalearth)


worldmap <- ne_countries(scale = 'medium', type = 'map_units', returnclass = 'sf')
head(worldmap[c('name', 'continent')])

ggplot() + geom_sf(data = worldmap) + theme_bw()

ggsave('zoom1_world.png')

#### select france ####

france <- worldmap[worldmap$name == 'France',]
ggplot() + geom_sf(data = france) + theme_bw()

ggsave('zoom2_france.png')


#### select europe ####

europe <- worldmap[worldmap$continent == 'Europe',]
ggplot() + geom_sf(data = europe) + theme_bw()

ggsave('zoom3_europe.png')

#### crop europe ####

europe_cropped <- st_crop(worldmap, xmin = -20, xmax = 45, ymin = 30, ymax = 73)
ggplot() + geom_sf(data = europe_cropped) + theme_bw()

ggsave('zoom4_europe.png')

#### crop europe -- variant: remove "gaps" / padding ####

ggplot() + geom_sf(data = europe_cropped) + coord_sf(expand = FALSE) + theme_bw()


#### crop europe -- variant: apply mollweide projection ####

ggplot() + geom_sf(data = europe_cropped) + coord_sf(crs = st_crs('+proj=moll')) + theme_bw()
ggsave('zoom4_europe_moll.png')


#### restrict display window on europe ####

ggplot() + geom_sf(data = worldmap) +
    coord_sf(xlim = c(-20, 45), ylim = c(30, 73), expand = FALSE) +
    theme_bw()
ggsave('zoom5_europe.png')

#### restrict display window on europe: fail ####

# this fails because xlim and ylim must be specified in the CRS that is used for projection (here: Mollweide)
ggplot() + geom_sf(data = worldmap) +
    coord_sf(crs = st_crs('+proj=moll'), xlim = c(-20, 45), ylim = c(30, 73), expand = FALSE) +
    theme_bw()
ggsave('zoom5_europe_fail.png')


#### restrict display window on europe using a different projection -- overview ####

target_crs <- '+proj=moll'
worldmap_trans <- st_transform(worldmap, crs = target_crs)

# display window corner points in the following order:
# top left, top right, bottom right, bottom left, top left (again, to close the shape)
# "crs = 4326" means that coordinates are supplied as WGS84 longitude/latitude coordinates
disp_win <- st_sfc(st_point(c(-20, 73)), st_point(c(45, 73)),
                   st_point(c(45, 30)), st_point(c(-20, 30)),
                   st_point(c(-20, 73)), crs = 4326)
disp_win

disp_win_trans <- st_transform(disp_win, crs = target_crs)
disp_win_lines <- st_sfc(st_linestring(st_coordinates(disp_win_trans)), crs = target_crs)

# corner points at bottom left and top right
# supply as WGS84 and directly transform to target CRS (Mollweide)
cornerpts <- data.frame(label = c('A', 'B'))
cornerpts$geometry <- st_transform(st_sfc(st_point(c(-20, 30)), st_point(c(45, 73)),
                                          crs = 4326),
                                   crs = target_crs)
cornerpts <- st_as_sf(cornerpts)

ggplot() + geom_sf(data = worldmap_trans) +
    geom_sf(data = disp_win_lines, color = 'red') +
    geom_sf_label(aes(label = label), data = cornerpts, color = 'red') +
    coord_sf(datum = target_crs) + theme_bw() + theme(axis.title = element_blank())

ggsave('zoom6_europe_moll_dispwin.png')

#### restrict display window on europe using a different projection ####

target_crs <- '+proj=moll'
worldmap_trans <- st_transform(worldmap, crs = target_crs)

disp_win_wgs84 <- st_sfc(st_point(c(-20, 30)), st_point(c(45, 73)), crs = 4326)
disp_win_wgs84

disp_win_trans <- st_transform(disp_win_wgs84, crs = target_crs)
disp_win_trans

disp_win_coord <- st_coordinates(disp_win_trans)

ggplot() + geom_sf(data = worldmap_trans) +
    coord_sf(xlim = disp_win_coord[,'X'], ylim = disp_win_coord[,'Y'], datum = target_crs, expand = FALSE) +
    theme_bw()

ggsave('zoom6_europe_moll.png')

#### zoom to berlin ####

zoom_to <- c(13.38, 52.52)  # Berlin

zoom_level <- 3

lon_span <- 360 / 2^zoom_level
lat_span <- 180 / 2^zoom_level

lon_bounds <- c(zoom_to[1] - lon_span / 2, zoom_to[1] + lon_span / 2)
lat_bounds <- c(zoom_to[2] - lat_span / 2, zoom_to[2] + lat_span / 2)


ggplot() + geom_sf(data = worldmap) +
    geom_sf(data = st_sfc(st_point(zoom_to), crs = 4326), color = 'red') +
    coord_sf(xlim = lon_bounds, ylim = lat_bounds) + theme_bw()

ggsave('zoom7_berlin.png')


#### zoom to kamchatka, different projection ####

zoom_to <- c(159.21, 56.41)  # ~ center of Kamchatka
zoom_level <- 2.5

# Lambert azimuthal equal-area projection around center of interest
target_crs <- sprintf('+proj=laea +lon_0=%f +lat_0=%f', zoom_to[1], zoom_to[2])


C <- 40075016.686   # ~ circumference of Earth in meters
x_span <- C / 2^zoom_level
y_span <- C / 2^(zoom_level+1)   # also sets aspect ratio

zoom_to_xy <- st_transform(st_sfc(st_point(zoom_to), crs = 4326), crs = target_crs)
zoom_to_xy

disp_window <- st_sfc(st_point(st_coordinates(zoom_to_xy - c(x_span / 2, y_span / 2))),
                      st_point(st_coordinates(zoom_to_xy + c(x_span / 2, y_span / 2))),
                      crs = target_crs)

ggplot() + geom_sf(data = worldmap) +
    geom_sf(data = zoom_to_xy, color = 'red') +
    coord_sf(xlim = st_coordinates(disp_window)[,'X'],
             ylim = st_coordinates(disp_window)[,'Y'],
             crs = target_crs,
             datum = target_crs) +
    theme_bw()

ggsave('zoom8_kamchatka.png')

# or, using the default graticule:

ggplot() + geom_sf(data = worldmap) +
    geom_sf(data = zoom_to_xy, color = 'red') +
    coord_sf(xlim = st_coordinates(disp_window)[,'X'],
             ylim = st_coordinates(disp_window)[,'Y'],
             crs = target_crs) +
    theme_bw()