lakes.R

library(rnaturalearth)
library(cartography)
library(sf)

lakes <- ne_download(scale = 10, type = "lakes", category = c("physical"), 
                     destdir = tempdir(), load = TRUE, returnclass = c("sf"))
countries <- ne_download(scale = 10, type = "countries", 
                         category = c("cultural"), destdir = tempdir(), 
                         load = TRUE, returnclass = c("sf"))
countries <- st_transform(countries, "ESRI:54017")
lakes <-  st_transform(lakes, "ESRI:54017")


lake_i <- st_intersection(lakes, countries)
lake_i$surf <- as.numeric(st_area(lake_i)) / 1000000
x <- aggregate(lake_i$surf, by = list(lake_i$ADMIN), FUN = length)
x$sh <- 100 * x$x / sum(x$x)
x <- x[order(x$x, decreasing = T),][1:5,]
x[3,1] <- "USA"
x$lab <- paste0(x$Group.1, " ", round(x$sh, 0), "%")

y <- aggregate(lake_i$surf, by = list(lake_i$ADMIN), FUN = sum)
y$sh <- 100 * y$x / sum(y$x)
y <- y[order(y$sh, decreasing = T),][1:5,]
y[2,1] <-'USA'
y[4,1] <- 'Tanzania'
y$lab <- paste0(y$Group.1, " ", round(y$sh, 0), "%")


logo2 <- png::readPNG("out.png")
pp <- dim(logo2)[2:1] * 20000
xref <- -16466598
yref <- -6769754


par(mar = c(0,0,0,0))
ghostLayer(lakes, bg = "ivory2")
plot(countries$geometry, col = "ivory4", border = "white",
     lwd = .5, lty = 3, add = TRUE)
plot(lakes$geometry, add = T, col = "lightblue", border = "blue", lwd = .5)
rasterImage(image = logo2, xleft = xref, ybottom = yref, 
            xright = xref + pp[1], ytop = yref + pp[2]) 
text(x = xref, y = 7520242,
     labels = '"Statistically, 60% of all lakes are located in Canada."', 
     font = 2, col = "black", cex = 1.1, adj = 0)
text(x = -10659541, y = 3715418,
     labels = 'Nope!', srt = 45, 
     font = 2, col = "red", cex = 2, adj = 0)
mtext("Source: Natural Earth, 1:10,000,000 countries and lakes - T. Giraud, 2020",
      side = 1, line = -1, 
      cex = .8, adj = c(.9))
text(-2985930,  -216075.5, 
     labels = paste0("# of lakes: \n",paste0(x$lab[1:5], collapse = "\n")), 
     adj = c(0,1), cex = .8)
text(5641697, -216075.5, 
     labels = paste0("surface: \n",paste0(y$lab[1:5], collapse = "\n")), 
     adj = c(0,1), cex = .8)