dancefloor.R

# code used for the following map : http://datagistips.blogspot.fr/

library(raster)
library(rgdal)
library(rgeos)

clc = raster("DATAS/CLC12_RIDF_RGF.tif") # here is the 200m Corine Land Cover GeoTiff
reg = readOGR("DATAS", "idf_geofla") # ile de france GéoFla Departments

# RECLASS
print("calculate")
clc2 = clc
values(clc2)[grep("^1.*$", values(clc))] = 1 # all raster values beginning with 1
values(clc2)[grep("^2.*$", values(clc))] = 2 # all raster values beginning with 2
values(clc2)[grep("^3.*$", values(clc))] = 3 # all raster values beginning with 3
values(clc2)[grep("^4.*$", values(clc))] = 4 # all raster values beginning with 4
values(clc2)[grep("^5.*$", values(clc))] = 5 # all raster values beginning with 5

# SIZES
width = diff(bbox(reg)[1,]) # width of idf
sizes = round(exp(seq(log(1000), log(diff(bbox(reg)[1,])), length.out=100))); plot(sizes) # the square sizes we choose

# CALCULATE
print("calculate")

print(Sys.time())
for (size in sizes) {
  
  output = file.path("OUT", paste(paste("grille", size, sep="_"), "shp", sep="."))
  
  if (!file.exists(output)) {
    
    print(Sys.time())
    print(size)
    
    grille = creerGrille(reg, size) # we create the grid
    
    counts = extractFrequencies(grille, clc2) # we get the number of pixels for each category
    
    RGB = round(prop.table(counts[, 1:3], 1)*255) # proportions * 255 : RGB code
    alpha = round(prop.table(counts[,4] + counts[,5])*255) # alpha (optional)
    
    grille$R = RGB[,1] # URBAN
    grille$G = RGB[,3] # NATURAL
    grille$B = RGB[,2] # AGRI
    grille$alpha = alpha # WATER
    
    writeOGR(grille, "OUT", paste("grille", size, sep="_"), "ESRI Shapefile", overwrite=T) # exporting data. The size is in the layer name.
    print(Sys.time())
  } 
}
print(Sys.time())

# DUPLICATE DATA
synoptic = makeAtlasLayer(reg, "size", sizes)
writeOGR(synoptic, "OUT", "synoptic", "ESRI Shapefile", overwrite=T)

# JUST A DUMMY ANIMATED RAINBOW COLORED (the one in the upper left side of the image)
par(bg = "white")

size=100
endSize = size/2

xmin = 0; ymin = 0; xmax = size; ymax = size

incs = log(seq(1, endSize, length.out=length(sizes))) # increments following log distribution
incs = endSize/max(incs) * incs # increments

cols = rainbow(35, start = 0, end = 1)
for (i in 1:length(incs)) {
  id = sizes[order(-sizes)][i]
  png(file.path("IMAGES/STRIP/", paste(id, ".png", sep="")), width=1000, height=1000) # the image path
  par(bg = NA)
  plot(c(xmin, xmax), c(ymin, ymax), type = "n", xlab = "", ylab = "", main = "", axes=FALSE)
  rect(xmin+incs[1:i], ymin+incs[1:i], xmax-incs[1:i], ymax-incs[1:i], col=cols, border=NA, lwd=.04) # decreasing size rectangles.
  dev.off()
}