Population Weighted Density - Lebanon exmple.R

require(raster, quietly=T)
require(dplyr, quietly=T, warn.conflicts=F)

download.file('https://data.worldpop.org/GIS/Population/Global_2000_2020/2020/LBN/lbn_ppp_2020_UNadj.tif',
              destfile = 'lebanon-pops.tif')
download.file('https://data.worldpop.org/GIS/Pixel_area/Global_2000_2020/LBN/lbn_px_area_100m.tif',
              destfile = 'lebanon-areas.tif')

pop_grid <- raster('lebanon-pops.tif')
area_grid <- raster('lebanon-areas.tif') / 1000000  # square meters to kms
density_grid <- pop_grid / area_grid

p <- values(pop_grid)
a <- values(area_grid)
d <- values(density_grid)
f <- !is.na(p)
d <- d[f]; a = a[f]; p = p[f]

# Conventional pop density
( pd <- sum(p) / sum(a) )

# PWD Arithmetic mean
( pwd_a <- sum(d * p) / sum(p) )

# PWD Geometric mean
p2 <- p[p > 0]  # all values must be positive for geometric mean
d2 <- d[d > 0]
( pwd_g <- exp(sum(p2 * log(d2)) / sum(p2)) )

# PWD percentiles
percentile = sum(p) / 100

df = tibble(pop = p, density = d) %>% 
  arrange(pop) %>% 
  mutate(cum_pop = cumsum(pop),
         pop_percentile = floor(cum_pop / percentile),
         percentile_diff = c(0, diff(pop_percentile))) %>% 
  filter(percentile_diff > 0)
print(df)

# PWD median
( pwd_m <- df$density[df$pop_percentile == 50] )

require(ggplot2)

df %>% ggplot(aes(pop_percentile, density)) + 
  geom_hline(yintercept = pd, col = 'orangered', size = 0.7) +
  annotate(x = 90, y = pd*1.3, label = paste('Pop density:', scales::comma(round(pd))), 
           geom = 'text', col = 'orangered', size = 4.5) +
  geom_hline(yintercept = pwd_g, col = 'royalblue1', size = 0.7) +
  annotate(x = 90, y = pwd_g*0.8, label = paste('PWD-G:', scales::comma(round(pwd_g))), 
           geom = 'text', col = 'royalblue1', size = 4.5) +
  geom_hline(yintercept = pwd_m, col = 'springgreen4', size = 0.7) +
  annotate(x = 90, y = pwd_m*1.3, label = paste('PWD-M:', scales::comma(round(pwd_m))), 
           geom = 'text', col = 'springgreen4', size = 4.5) +
  geom_hline(yintercept = pwd_a, col = 'violetred3', size = 0.7) +
  annotate(x = 90, y = pwd_a*1.3, label = paste('PWD-A:', scales::comma(round(pwd_a))), 
           geom = 'text', col = 'violetred3', size = 4.5) +
  geom_line() + geom_point() +
  scale_y_log10(labels = scales::comma) +
  labs(title = 'Lebanon population density by population percentile, 2020',
       subtitle = 'Data: worldpop.org, 3 arcsec grid',
       x = 'Percentile', y = 'Density (pop/km, log scale)')