Paris.R

library(sf)
library(spatstat)
library(sp)
library(maptools)
library(raster)
library(cartography)
library(SpatialPosition)

## import dataset
feat <- sf::st_read("https://gist.githubusercontent.com/rCarto/747164575e3f216a123c3092d0ce9162/raw/f12390464f255b5f9760c577ab6bf5456cf61a40/iris75.geojson")
# Use french projection
feat <- sf::st_transform(feat, 2154)

## Compute density raster
# from sf to sp
feat <- as(feat,'Spatial')
# get coodinates
coords <- coordinates(feat)
# from sp to spatstat
pts <- ppp(x = coords[,1], y = coords[,2], window = as.owin(feat, 10), 
           marks = feat$P14_POP)
# compute density
ds <- density.ppp(x = pts, sigma = 300, eps = 100, weights = pts$marks)
# spatstat to raster, in inhab per sq. km
ras <- raster(ds, crs = proj4string(feat)) * 1000 * 1000
ras[is.na(ras)] <- 0

## Compute contour polygons
# raster break values 
bks <- c(seq(0,50000, 5000), 54100)
iso_pop <- rasterToContourPoly(r = ras, breaks = bks, mask = feat)

## prepare for tanaka contours
# from sp to sf
iso_pop <- st_as_sf(iso_pop)
# order iso surfaces
iso_pop <- iso_pop[order(iso_pop$center),]
# color palette à la viridis
pal <-c("#000004FF", "#170C3AFF", "#420A68FF", "#6B186EFF", "#932667FF", 
        "#BB3754FF", "#DD513AFF", "#F3771AFF", "#FCA50AFF", "#F6D645FF", 
        "#FCFFA4FF")
########### Figure 3
png(filename = "paris3.png", width = 800, height = 500, res = 100, bg = NA)
# custom margin
par(mar = c(0,0,1.2,0))
# plot area
plot(st_geometry(iso_pop), col = NA, border = NA, bg = "ivory1")
## tanaka contours + colors
# plot each contour layer successively
for(i in 1:nrow(iso_pop)){
  p <- st_geometry(iso_pop[i,])
  # plot light contour polygons with a Nort-West shift
  plot(p + c(-30, 30), add=T, border = "#ffffff90",col = "#ffffff90")
  # plot dark contour polygons with a South-East shift
  plot(p + c(50, -50),  col = "#00000099", add=T, border = "#00000099")
  # plot colored contour polygons in place
  plot(p, col = pal[i], border = "NA", add=T)
}
# legend
legendChoro(pos = "topright", breaks = bks, col = pal, nodata = F,
            title.txt = "Inhabitants\nper sq. km *", cex = 0.8)
# layout
layoutLayer(title = "Smoothed Population Density, Paris 2014",
            col = "ivory1", tabtitle = T, coltitle = "black",
            frame = T, scale = 1,
            sources = "T. Giraud - 2018",
            author = "Contours...Iris® - IGN 2017, Recensements de la population - Insee 2017")
north(pos = c(661000,6857900 ))
# annotations
text(x = 654500, y = 6856900, adj = 0, cex = 0.6, font = 3,
     labels='(*) Kernel Density Estimation with\n    a gaussian kernel (sigma = 300 m)')
dev.off()