Skip to content

Instantly share code, notes, and snippets.

@rCarto rCarto/Paris.R
Last active Jul 21, 2018

Embed
What would you like to do?
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.