world_maps.Rmd

---
title: "World maps"
output:
  html_document:
    df_print: paged
---

```{r echo = FALSE, message = FALSE}
library(tidyverse)
library(sf)
```

## Longitude--latitude

```{r world-longlat, fig.width = 8.5}
world_sf <- sf::st_as_sf(rworldmap::getMap(resolution = "low"))

crs_longlat <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

ggplot(world_sf) + 
  geom_sf(fill = "#E69F00B0", color = "black", size = 0.5/.pt) + 
  coord_sf(expand = FALSE, crs = crs_longlat) + 
  scale_x_continuous(name = "longitude", breaks = seq(-120, 120, by = 60)) +
  scale_y_continuous(name = "latitude", breaks = seq(-60, 60, by = 30)) +
  ggtitle("Cartesian longitude and latitude") +
  theme_minimal() +
  theme(
    panel.background = element_rect(fill = "#56B4E950"),
    panel.grid.major = element_line(color = "gray30", size = 0.25),
    axis.ticks = element_line(color = "gray30", size = 0.5/.pt)
  )
```

## Mercator

```{r world-mercator, fig.width = 8.5}
crs_mercator <- "+proj=merc"

# calculate bounding box in transformed coordinates
mercator_bbox <- 
  rbind(c(-180, -85), c(180, 85)) %>%
  st_multipoint() %>%
  st_sfc(
    crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
  ) %>% 
  st_transform(crs = crs_mercator)

ggplot(world_sf) + 
  geom_sf(fill = "#E69F00B0", color = "black", size = 0.5/.pt) + 
  scale_x_continuous(name = NULL, breaks = seq(-120, 120, by = 60)) +
  scale_y_continuous(name = NULL, breaks = seq(-60, 60, by = 30)) +
  coord_sf(
    xlim = mercator_bbox[[1]][, 1],
    ylim = mercator_bbox[[1]][, 2],
    expand = FALSE,
    crs = crs_mercator
  ) + 
  ggtitle("Mercator") +
  theme_minimal() +
  theme(
    panel.background = element_rect(fill = "#56B4E950", color = "white", size = 1),
    panel.grid.major = element_line(color = "gray30", size = 0.25)
  )
```


## Robinson

```{r world-robin, fig.asp = 0.6}
crs_robin <- "+proj=robin +lat_0=0 +lon_0=0 +x0=0 +y0=0"

# projection outline in long-lat coordinates
lats <- c(90:-90, -90:90, 90)
longs <- c(rep(c(180, -180), each = 181), 180)

robin_outline <- 
  list(cbind(longs, lats)) %>%
  st_polygon() %>%
  st_sfc(
    crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
  ) %>% 
  st_transform(crs = crs_robin)

# bounding box in transformed coordinates
xlim <- c(-18494733, 18613795)
ylim <- c(-9473396, 9188587)
robin_bbox <- 
  list(
    cbind(
      c(xlim[1], xlim[2], xlim[2], xlim[1], xlim[1]), 
      c(ylim[1], ylim[1], ylim[2], ylim[2], ylim[1])
    )
  ) %>%
  st_polygon() %>%
  st_sfc(crs = crs_robin)

# area outside the earth outline
robin_without <- st_difference(robin_bbox, robin_outline)

ggplot(world_sf) + 
  geom_sf(fill = "#E69F00B0", color = "black", size = 0.5/.pt) + 
  geom_sf(data = robin_without, fill = "white", color = NA) +
  geom_sf(data = robin_outline, fill = NA, color = "grey30", size = 0.5/.pt) +
  scale_x_continuous(name = NULL, breaks = seq(-120, 120, by = 60)) +
  scale_y_continuous(name = NULL, breaks = seq(-60, 60, by = 30)) +
  coord_sf(xlim = 0.95*xlim, ylim = 0.95*ylim, expand = FALSE, crs = crs_robin, ndiscr = 1000) + 
  ggtitle("Robinson") +
  theme_minimal() +
  theme(
    panel.background = element_rect(fill = "#56B4E950", color = "white", size = 1),
    panel.grid.major = element_line(color = "gray30", size = 0.25)
  )
```


## Interrupted Goode homolosine

```{r world-goode, fig.width = 8.5}
crs_goode <- "+proj=igh"

# projection outline in long-lat coordinates
lats <- c(
  90:-90, # right side down
  -90:0, 0:-90, # third cut bottom
  -90:0, 0:-90, # second cut bottom
  -90:0, 0:-90, # first cut bottom
  -90:90, # left side up
  90:0, 0:90, # cut top
  90 # close
)
longs <- c(
  rep(180, 181), # right side down
  rep(c(80.01, 79.99), each = 91), # third cut bottom
  rep(c(-19.99, -20.01), each = 91), # second cut bottom
  rep(c(-99.99, -100.01), each = 91), # first cut bottom
  rep(-180, 181), # left side up
  rep(c(-40.01, -39.99), each = 91), # cut top
  180 # close
)

goode_outline <- 
  list(cbind(longs, lats)) %>%
  st_polygon() %>%
  st_sfc(
    crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
  ) %>% 
  st_transform(crs = crs_goode)

# bounding box in transformed coordinates
xlim <- c(-21945470, 21963330)
ylim <- c(-9538022, 9266738)
goode_bbox <- 
  list(
    cbind(
      c(xlim[1], xlim[2], xlim[2], xlim[1], xlim[1]), 
      c(ylim[1], ylim[1], ylim[2], ylim[2], ylim[1])
    )
  ) %>%
  st_polygon() %>%
  st_sfc(crs = crs_goode)

# area outside the earth outline
goode_without <- st_difference(goode_bbox, goode_outline)

ggplot(world_sf) + 
  geom_sf(fill = "#E69F00B0", color = "black", size = 0.5/.pt) +
  geom_sf(data = goode_without, fill = "white", color = NA) +
  geom_sf(data = goode_outline, fill = NA, color = "grey30", size = 0.5/.pt) +
  scale_x_continuous(name = NULL, breaks = seq(-120, 120, by = 60)) +
  scale_y_continuous(name = NULL, breaks = seq(-60, 60, by = 30)) +
  coord_sf(xlim = 0.95*xlim, ylim = 0.95*ylim, expand = FALSE, crs = crs_goode, ndiscr = 1000) + 
  ggtitle("Interrupted Goode homolosine") +
  theme_minimal() +
  theme(
    panel.background = element_rect(fill = "#56B4E950", color = "white", size = 1),
    panel.grid.major = element_line(color = "gray30", size = 0.25)
  )
```


## Winkel tripel

```{r world-winkel-tri, fig.width = 8.5}
# Winkel tripel projection needs to be done manually, is not
# supported by sf.
crs_wintri <- "+proj=wintri +datum=WGS84 +no_defs +over"

# world
world_wintri <- lwgeom::st_transform_proj(world_sf, crs = crs_wintri)

# graticule
grat_wintri <- sf::st_graticule(lat = c(-89.9, seq(-80, 80, 20), 89.9))
grat_wintri <- lwgeom::st_transform_proj(grat_wintri, crs = crs_wintri)

# earth outline
lats <- c(90:-90, -90:90, 90)
longs <- c(rep(c(180, -180), each = 181), 180)
wintri_outline <- 
  list(cbind(longs, lats)) %>%
  st_polygon() %>%
  st_sfc(
    crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
  ) %>% 
  lwgeom::st_transform_proj(crs = crs_wintri)

ggplot() + 
  geom_sf(data = wintri_outline, fill = "#56B4E950", color = NA) +
  geom_sf(data = grat_wintri, color = "gray30", size = 0.25/.pt) + 
  geom_sf(data = world_wintri, fill = "#E69F00B0", color = "black", size = 0.5/.pt) + 
  geom_sf(data = wintri_outline, fill = NA, color = "grey30", size = 0.5/.pt) +
  ggtitle("Winkel tripel") +
  coord_sf(datum = NA) +
  theme_minimal()
```