first commit:

import from MAGEO course
nested grid function updated

.gitignore

@@ -0,0 +1,4 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata

R/as.nb.sgbp.R

@@ -0,0 +1,9 @@
+# Conversion function suggested by R Bivand between sgbp and nb objects
+# https://cran.r-hub.io/web/packages/spdep/vignettes/nb_sf.html
+as.nb.sgbp <- function(x, ...) {
+  attrs <- attributes(x)
+  x <- lapply(x, function(i) { if(length(i) == 0L) 0L else i } )
+  attributes(x) <- attrs
+  class(x) <- "nb"
+  x
+}

R/get.lux.data.R

@@ -0,0 +1,87 @@
+# get.lux.data.R
+#
+# Get some open data per municipalities of Luxembourg----
+
+## Polygons----
+# Download and convert to sf communal polygon vector data (geojson) from geoportail of Luxembourg:
+# source: https://data.public.lu/fr/datasets/limites-administratives-du-grand-duche-de-luxembourg/
+url<-"https://data.public.lu/fr/datasets/r/16103fa4-7ff1-486a-88bc-5018353957ea"
+lux102sf<-sf::st_read(url) #102 because it is the latest vintage with 102 communes, crs is 4326
+#save as geopackage for future use
+sf::st_write(lux102sf,"data/Communes102_4326.gpkg", delete_dsn = TRUE)
+
+## Attributes----
+# Download population density since 1821! from STATEC
+apiURL_1821_2023 <- "https://lustat.statec.lu/rest/data/LU1,DF_X020,1.0/.A?startPeriod=1821&endPeriod=2023"
+data_1821_2023 <- data.frame(rsdmx::readSDMX(apiURL_1821_2023))
+saveRDS(data_1821_2023, "data/data_1821_2023.rds")
+
+## Cleaning attributes----
+# For some reason, the 'libelles' is not given although asked for on the website,
+# only a specification code is given to uniquely identify communes (and cantons)
+# Hopefully the "libelles" are present via csv download.
+# We use a csv with population 2015 to 2023, downloaded manually from same source
+lu1<-read.csv("data/LU1_DF_X021_1.0_A..csv") 
+
+lu<-unique(lu1[,"SPECIFICATION..Spécification"])
+#from which we see starting characters "CT" is for Cantons, which we remove
+lu102<-lu[stringr::str_sub(lu, 1, 2)=="CM"]
+#lu102 has now 102 communes.
+
+#A unique id is available with 7 first characters,
+# while commune names start from character 10.
+# We thus split lu102 in 2 columns``
+CMCODE<-stringr::str_sub(lu102, 1, 7)
+CMNAME<-stringr::str_sub(lu102, 10)
+CMLU102<-data.frame(cbind(CMCODE,CMNAME))
+
+# which we can merge to  columns of interest in data_1821_2023,
+#, which is in long format (years as rows)
+data1821sub<-data_1821_2023[,c("SPECIFICATION","obsTime","obsValue")]
+# and removing also the cantons and country aggregates (selecting CM only)
+data1821sub2<-data1821sub[stringr::str_sub(data1821sub[,"SPECIFICATION"], 1, 2)=="CM",]
+
+df1821m<-merge(data1821sub2,CMLU102,by.x="SPECIFICATION", by.y="CMCODE", all.x = TRUE,sort=FALSE)
+#reorganize and rename SPECIFICATION as CMCODE, obsValue as Density and obsTime as Year
+df1821mo<-df1821m[,c(1,4,2,3)]
+names(df1821mo)<-c("CMCODE","CMNAME","Year","Density")
+#Density as numeric
+df1821mo$Density<-as.numeric(df1821mo$Density)
+
+#Cast attribute table (df1821mo) in wide format and adapt column names
+df1821wide<-reshape2::dcast(df1821mo, CMCODE + CMNAME ~ Year, value.var="Density")
+names(df1821wide)[-c(1,2)]<-paste0("Density",names(df1821wide)[-c(1,2)])
+
+## Merge with sf----
+
+# Note that the CMCODE does not match the European LAU2 convention used by the 
+#  geoportail for the spatial data hence our lux102sf. Matching
+#  needs to be made on names before matching CMCODE with LAU2.
+#  But names are dangerous beasts, and we need to worked 3 differences
+#   (out of 102) manually:
+setdiff(unique(df1821wide$CMNAME), unique(lux102sf$COMMUNE))
+#[1] "Lac de la Haute-Sûre" "Rosport - Mompach"   
+#[3] "Redange-sur-Attert"  
+setdiff(unique(lux102sf$COMMUNE),unique(df1821wide$CMNAME))
+#[1] "Lac de la Haute Sûre" "Rosport-Mompach"     
+#[3] "Redange"
+
+# Manual replacement (in attributes table, i.e. using names as of geoportail)
+df1821<-df1821wide
+df1821[df1821[,"CMNAME"]=="Lac de la Haute-Sûre","CMNAME"]<-"Lac de la Haute Sûre"
+df1821[df1821[,"CMNAME"]=="Rosport - Mompach","CMNAME"]<-"Rosport-Mompach"
+df1821[df1821[,"CMNAME"]=="Redange-sur-Attert","CMNAME"]<-"Redange"
+
+#merge with corrected names
+lux102sf_density1821<-merge(lux102sf,df1821,by.x="COMMUNE",by.y="CMNAME")
+
+## Save ----
+#   sf with attributes for use in R
+saveRDS(lux102sf_density1821,"data/lux102sf_density1821.rds")
+#   sf as geopackage with attributes
+sf::st_write(lux102sf_density1821,"data/lux102_4326_density1821.gpkg",delete_dsn=TRUE)
+#   LAU2 to CMCODE correspondance table (with names) as a csv
+write.csv(sf::st_drop_geometry(lux102sf_density1821[,c("COMMUNE","LAU2","CMCODE")]),"data/LAU2_CMCODE_table.csv")
+
+
+

R/ggplot.themap.R

@@ -0,0 +1,66 @@
+#' ggplot.themap.R
+#' 
+#' A function to make choropleth maps for continuous data using
+#' an sf input where the attributes to be mapped are included in the dataframe.
+#' Uses ggplot2 (scale_fill_manual) and classInt for flexible discretization
+#' (and allowing for own (fixed) breaks and colours to be set)
+#' 
+#' @param sf A polygon sf including the variable to be mapped (choropleth)
+#' @param varname Name of variable to be mapped
+#' @param n Number of classes 
+#' @param style A discretization style as used in classInt::classIntervals
+#' @param fixedBreaks If style is "fixed", breaks are chosen and input here
+#'  (similar to classInt::classIntervals)
+#' @param low.colour Polygon fill colour for lower end values.
+#'  Other colours interpolated up to high.colour. See classInt::findColours
+#' @param high.colour Polygon fill colour for upper end values
+#' @param cl.colours To supply a full list of n colours rather than low and high colours
+#' @param outline.colour Polygon outline colour
+#' @param outline.width Polygon outline width
+#' @param n.digits Digits for legend display
+#' @param leg.title Character string for legend title
+#' @param main.title Character string for map title
+#' @param sub.title Character string for map sub title
+#' @param ggtheme ggplot theme
+#'
+#' @return A ggplot object using geom_sf for polygons and scale_fill_manual
+#'  for rendering the discretized values
+#'  
+#' @export
+#' @examples See ggplot.themap.examples.R
+#'
+ggplot.themap<-function(sf,varname,
+                        n=5, style='quantile', fixedBreaks=NULL,
+                        low.colour="lightyellow",
+                        high.colour="darkred",
+                        cl.colours=attr(classInt::findColours(cl.intvl, c(low.colour,high.colour)),"palette"),
+                        outline.colour="#ffce00",
+                        outline.width=0.2,
+                        n.digits=2,
+                        leg.title=varname,
+                        main.title=paste(substitute(sf),varname),
+                        sub.title=attr(cl.intvl, "style"),
+                        ggtheme=ggplot2::theme_bw()
+                        ){
+  svar<-sf::st_drop_geometry(sf)[,varname]
+  cl.intvl<-classInt::classIntervals(svar, n=n, style=style, fixedBreaks=fixedBreaks)
+  n.classes<-length(cl.intvl$brks)-1
+  cl.value<-factor(classInt::findCols(cl.intvl))
+  leg.labels<-paste(format(round(cl.intvl$brks[1:n.classes],digits=n.digits), nsmall=2),
+                    format(round(cl.intvl$brks[2:(n.classes+1)],digits=n.digits), nsmall=2),
+                    sep=" - ")
+  sc.f<-ggplot2::scale_fill_manual(name = leg.title,
+                          breaks = seq(1:n.classes),
+                          values=cl.colours,
+                          labels=leg.labels)
+  themap<-ggplot2::ggplot()
+  themap<-themap+ggplot2::geom_sf(data=sf,ggplot2::aes(fill=cl.value),
+            colour=outline.colour,
+            size=outline.width)
+  themap<-themap+sc.f
+  themap<-themap+ggtheme
+  themap<-themap+ggplot2::ggtitle(label=main.title, subtitle=sub.title)
+
+  return(themap)
+}
+

R/ggplot.themap.f.R

@@ -0,0 +1,48 @@
+#' ggplot.themap.f.R
+#'
+#' A function to make choropleth maps for categorical data using
+#' an sf input where the attributes to be mapped are included in the dataframe.
+#' Uses ggplot2 in a similar way as ggplot.themap (for discretized continuous
+#'  data) using ggplot2::scale_fill_manual
+#' 
+#' @param sf A polygon sf including the variable to be mapped (choropleth)
+#' @param varname Name of variable to be mapped
+#' @param cl.colours To supply a vector of colours. Named vector using
+#'  categories values as names to have a match of ex-ante defined colours 
+#' @param outline.colour Polygon outline colour
+#' @param outline.width Polygon outline width
+#' @param leg.title Character string for legend title
+#' @param main.title Character string for map title
+#' @param sub.title Character string for map sub title
+#' @param ggtheme ggplot theme
+#'
+#' @return
+#' @export
+#' @examples See ggplot.themap.f.examples.R
+#'
+ggplot.themap.f<-function(sf,varname,
+                          cl.colours=NULL,
+                          outline.colour="#ffce00",
+                          outline.width=0.2,
+                          leg.title=varname,
+                          main.title=paste(substitute(sf),varname),
+                          sub.title=NULL,
+                          ggtheme=ggplot2::theme_bw()
+){              
+  svar<-sf::st_drop_geometry(sf)[,varname]
+  sc.f<-if (!is.null(cl.colours)) { #check if there is a color palette provided
+    ggplot2::scale_fill_manual(values=cl.colours, name = leg.title)
+  } else {
+    ggplot2::scale_fill_hue(name = leg.title)
+  }
+  themap<-ggplot2::ggplot()
+  themap<-themap+ggplot2::geom_sf(data=sf,ggplot2::aes(fill=svar),
+            colour=outline.colour,
+            size=outline.width)
+  themap<-themap+sc.f
+  themap<-themap+ggtheme
+  themap<-themap+ggplot2::ggtitle(label=main.title,subtitle=sub.title)
+
+  return(themap)
+}
+

R/ggplot.themap.f.lux.R

@@ -0,0 +1,62 @@
+#' ggplot.themap.f.examples.R
+#' 
+#' Examples of using ggplot.themap.f with Luxembourg data
+
+#function
+source("R/ggplot.themap.f.R")
+
+#data
+lux102sf_density1821<-readRDS("data/lux102sf_density1821.rds")
+
+p<-ggplot.themap.f(lux102sf_density1821,"CANTON") #default colours and titles
+p
+
+# with Color Brewer or Lego bricks colours
+pb<-ggplot.themap.f(lux102sf_density1821,"CANTON",
+                cl.colours = RColorBrewer::brewer.pal(12, "Set3"),
+                main.title="Luxembourg cantons",leg.title=NULL)
+pb
+
+plego<-ggplot.themap.f(lux102sf_density1821,"CANTON",
+                cl.colours = legocolors::legocolors[2:13,]$hex,
+                outline.colour="white",
+                main.title="Luxembourg cantons",leg.title=NULL)
+plego
+
+#If a named vector is provided, enforces a match of values and colours
+pb2<-ggplot.themap.f(lux102sf_density1821, "DISTRICT", cl.colours = RColorBrewer::brewer.pal(3, "Set3"))
+pb2
+
+district.colours=c("Grevenmacher" = "darkolivegreen1", "Luxembourg" = "darkolivegreen3", "Diekirch"="darkolivegreen4", "Wallonie"="orange") 
+pdistr<-ggplot.themap.f(lux102sf_density1821, "DISTRICT",
+                        cl.colours = district.colours,
+                        main.title="Luxembourg former districts", leg.title=NULL)
+#As we know Wallonie is not a valid district of Luxembourg and is thus ignored.
+# This is going to be useful for mapping when you don't know before if a category is present or not (e.g. LISA maps)
+pdistr
+
+#Since the retuned outpu is a ggplot object, we can still modify ex-post:
+#For example removing the legend:
+pdistr+ theme(legend.position = "none")
+# This case is useful if one maps many categories,
+#  such as all communes
+pcom<-ggplot.themap.f(lux102sf_density1821,"LAU2",
+                        main.title="Luxembourg 102 communes")
+#The legend would take up the whole map space, so:
+pcom<-pcom+ theme(legend.position = "none")
+
+#Note, theoretically the map could work with only 4 colours.
+#The minimum number of colours is implemented in tmap.
+# tm_shape(lux102sf_density1821) +
+#  tm_polygons(col = "MAP_COLORS", minimize = TRUE)
+
+
+pdf(file="output/Lux_atlas_categorical.pdf")
+print(p)
+print(pb)
+print(plego)
+print(pb2)
+print(pdistr)
+print(pcom)
+dev.off()
+

R/ggplot.themap.lux.R

@@ -0,0 +1,83 @@
+#' ggplot.themap.examples.R
+#' 
+#' Examples of using ggplot.themap with Luxembourg data
+
+#
+library(ggplot2)
+
+#function
+source("R/ggplot.themap.R")
+
+#data
+lux102sf_density1821<-readRDS("data/lux102sf_density1821.rds")
+
+# Compare ggplot2 basic mapping with this function
+g<-ggplot()+
+  geom_sf(data=lux102sf_density1821,
+          aes(fill=Density2023))
+g
+
+p<-ggplot.themap(lux102sf_density1821,"Density1821")
+p
+
+# Improved map titles
+myvar<-"Density1821"
+p<-ggplot.themap(lux102sf_density1821,myvar,
+                 leg.title= "inh. /sq. km",
+                 main.title=paste("Population density",
+                                  stringr::str_sub(myvar, 8)))
+p
+
+# all into pdf - quantile maps
+pdf(file="output/Lux_atlas_density_qt.pdf")
+for (i in 6:62) {
+  myvar<-names(lux102sf_density1821)[i]
+  p<-ggplot.themap(lux102sf_density1821,myvar,
+                   leg.title= "inh. /sq. km",
+                   main.title=paste("Population density",
+                                    stringr::str_sub(myvar, 8)))
+  print(p)
+}
+dev.off()
+
+# all into pdf - jenks maps
+pdf(file="output/Lux_atlas_density_jenks.pdf")
+for (i in 6:62) {
+  myvar<-names(lux102sf_density1821)[i]
+  p<-ggplot.themap(lux102sf_density1821,myvar, n=6,
+                   leg.title= "inh. /sq. km",
+                   main.title=paste("Population density",
+                                    stringr::str_sub(myvar, 8)))
+  print(p)
+}
+dev.off()
+
+# all into pdf - fixed breaks to visualize increase across time with
+#  single legend. ALso example of colorbrewer for palette
+
+#make fixed breaks from log sequence
+mybrks<-exp(seq(log(15), log(2600), length.out = 11)) #reproduced from lseq in emdbook pkg
+
+pdf(file="output/Lux_atlas_density_fixed.pdf")
+for (i in 6:62) {
+  myvar<-names(lux102sf_density1821)[i]
+  p<-ggplot.themap(lux102sf_density1821,myvar,
+                   style="fixed",fixedBreaks=mybrks,
+                   cl.colours=rev(RColorBrewer::brewer.pal(11, "Spectral")),
+                   leg.title= "inh. /sq. km",
+                   main.title=paste("Population density",
+                                    stringr::str_sub(myvar, 8)))
+  print(p)
+}
+dev.off()
+
+# Animated gif from pdf (written next to file)
+pdf2animgif<-function(pdfpath,fps=1){
+  mypdf<-magick::image_read_pdf(paste0(pdfpath,".pdf"))
+  mygif<-magick::image_animate(mypdf,fps = fps)
+  magick::image_write(mygif, path = paste0(pdfpath,".gif"))
+}
+
+pdf2animgif("output/Lux_atlas_density_fixed")
+
+