tanakapop.R

library(raster)
library(sf)
library(viridis)
library(cartography)
library(tanaka)

temp <- tempfile()
data_url <- "http://cidportal.jrc.ec.europa.eu/ftp/jrc-opendata/GHSL/GHS_POP_GPW4_GLOBE_R2015A/GHS_POP_GPW42015_GLOBE_R2015A_54009_1k/V1-0/GHS_POP_GPW42015_GLOBE_R2015A_54009_1k_v1_0.zip"
download.file(data_url, temp)
unzip(temp, exdir = "pop")
pop2015 <- raster("pop/GHS_POP_GPW42015_GLOBE_R2015A_54009_1k_v1_0/GHS_POP_GPW42015_GLOBE_R2015A_54009_1k_v1_0.tif")
center <- st_as_sf(data.frame(x=425483.8, y=5608290), 
                   coords=(c("x","y")), crs = st_crs(pop2015))
center <- st_buffer(center, dist = 800000)
ras <- crop(pop2015, st_bbox(center)[c(1,3,2,4)])
mat <- focalWeight(x = ras, d = c(10000), type = "Gauss")
rassmooth <- focal(x = ras, w = mat, fun = sum, pad = TRUE, padValue = 30)
bks <- c(0,25,50,100,250,500,750,1000,1750,2500,5000, 7500,10000)
png(filename = "circle.png", width = 800, height = 700, res = 100)
par(mar = c(0,0,1.2,0))
tanaka(x = rassmooth, 
       breaks = bks, 
       mask = center, 
       col = inferno(12),shift = 2500,
       legend.pos = "topleft",
       legend.title = "Inhabitants\nper km2")
plot(st_geometry(center), add = T, border = "white", lwd = 6)
layoutLayer(title = "Smoothed Population Density", 
            author = 'Data : European Commission, Joint Research Centre (JRC); Columbia University, CIESIN (2015): GHS population grid, derived from GPW4.', 
            sources = 'T. Giraud, 2019', scale = F, frame = F, tabtitle = TRUE)
text(-374516.2 ,6408290.0, "Gaussian smoothing, sigma = 10km", adj = 0, font = 3, cex = .8 )
dev.off()