waffle.R

library(sf)
library(raster)
library(dplyr)
library(tidyr)


waffle <- function(f, the_stats, the_res = 10000) {
  
  ## Create Raster
  r <- raster(extent(f))
  res(r) <- the_res
  values(r) <- NA
  f.r <- rasterize(f, r)
  
  ## Polygonize
  f.pols <- rasterToPolygons(f.r)
  
  ## Get coordinates
  f.pols <- st_as_sf(f.pols)
  coords <- f.pols %>% st_centroid %>% st_coordinates
  
  f.pols$x <- coords[, 1]
  f.pols$y <- coords[, 2]
  
  ## Sort features horizontally then vertically
  ## (Vertically or horizontally could be chosen by the user)
  f.pols <- f.pols %>% arrange(x, y)
  
  ## Number of squares for each class
  pc <- the_stats / sum(the_stats) ## Get percentages
  n <- floor(pc * nrow(f.pols))
  n <- cumsum(n)
  n <- c(1, n)
  to_add <- nrow(f.pols) - max(n) ## To have the same count of polygons
  n[length(n)] <- n[length(n)] + to_add
  # max(n) == nrow(f.pols) # > TRUE
  
  ## Add Class
  f.pols$classe <- NA
  start <- 0
  for(i in 2:(length(n))) {
    classe <- names(n)[i]
    message("Classe ", classe)
    from <- n[i-1]
    to <- n[i]
    f.pols$classe[from:to] <- rep(classe, length(from:to))
  }
  
  return(f.pols)
}

##=##=##=##=##=##=##=##=
## Get Land Usage Stats
##=##=##=##=##=##=##=##=

# https://www.statistiques.developpement-durable.gouv.fr/corine-land-cover-0
download.file("https://www.statistiques.developpement-durable.gouv.fr/sites/default/files/2019-07/clc_etat_com_n1.csv", "stats_occsol.csv")

## Read
f.stats <- read.table("stats_occsol.csv", sep=";", quote = "\"", header = TRUE)

## Skip the 3 first columns
f.stats <- f.stats[4:nrow(f.stats), ]

## Rename
names(f.stats)[4:8] <- c("artificialise", "agricole", "semi_naturel", "humide", "surface_eau")

## Convert to numeric
for(i in 4:8) {
  f.stats[[i]] <- f.stats[[i]] %>% as.character %>% as.numeric
}

## Calculate the stats
the_stats <- c("artificialise" = sum(f.stats$artificialise), 
                  "agricole" = sum(f.stats$agricole), 
                  "semi_naturel" = sum(f.stats$semi_naturel), 
                  "humide" = sum(f.stats$humide),
                  "surface_eau" = sum(f.stats$surface_eau)) # Sum areas


##=##=##=##=##=##=##=##=
# Process
##=##=##=##=##=##=##=##=

## Read Geo Data file
f <- st_read("C:/Users/mathieu/Documents/data/ADMIN-EXPRESS-COG_2-0__SHP__FRA_L93_2019-09-24/ADMIN-EXPRESS-COG_2-0__SHP__FRA_2019-09-24/ADMIN-EXPRESS-COG/1_DONNEES_LIVRAISON_2019-09-24/ADE-COG_2-0_SHP_LAMB93_FR/DEPARTEMENT.shp")

## Waffle !!
f.pols <- waffle(f, the_stats)

## Export
st_write(f.pols, "f.pols.gpkg", delete_layer = TRUE, delete_dsn = TRUE)