Skip to content

Instantly share code, notes, and snippets.

@datagistips
Last active April 19, 2017 22:01
Show Gist options
  • Save datagistips/7110caa1b9c15de769d2 to your computer and use it in GitHub Desktop.
Save datagistips/7110caa1b9c15de769d2 to your computer and use it in GitHub Desktop.
# 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()
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment