cartograms.R

library(rgeos)
library(maptools)
 
load("F:/DATAS/DEPARTEMENT.rda")
 
foret <- read.csv2("IN/foret.csv")
 
 
######################
## SCALED CARTOGRAM ##
######################
 
shrinkExactlyThat <- function(pols, scl) {
  
  out = vector(mode="list", length=nrow(pols))
  
  for (i in seq(along=out)) {
    pol <- pols[i, ]
    scld <- elide(pol, scale=diff(bbox(g)[1, ])*scl[i]) # scale data
    shft <- coordinates(pol) - coordinates(scld) 
    out[[i]] <- elide(scld, shift=shft) # shift to original coordinates
    row.names(out[[i]]) <- as.character(i)
  }
  res <- do.call("rbind", out)
  resdf <- SpatialPolygonsDataFrame(res, data=pols@data)
  return(res)
}
 
 
######################
## ERODED CARTOGRAM ##
######################
 
shrinkThat <- function(pols, ratios, inc) {
  
  out <- vector(mode="list", length = nrow(pols))
  for (j in seq(along=out)) {
    print(paste("traitement du polygone:", j))
    
    pol <- pols[j, ]
    
    goOn <- TRUE
    i <- 1
    diff <- list() 
    
    while(goOn) {
      diff[[i]] <- abs(gArea(pol)*ratios[j] - gArea(gBuffer(pol, width=-i*inc)))
      if (i > 1) {
        if (diff[[i]] > diff[[i-1]]) {
          amount <- -((i-1)*inc)
          break
        }
      }
      i <- i+1
    }  
    out[[j]] <- gBuffer(pol, width=amount)
    row.names(out[[j]]) <- row.names(pol)
  }
  res <- do.call("rbind", out)
  resdf <- SpatialPolygonsDataFrame(res, data=pols@data)
  return(resdf)
}
 
 
#############################
## INTEGRATION DES DONNEES ##
#############################
 
m <- match(deps$CODE_DEPT, foret$X)
deps$Taux <- foret$Taux.de.boisement..en...de.la.surface.totale[m]
ratios <- deps$Taux/100
 
 
###########################
## NON CONTIG CARTOGRAMS ##
###########################
 
shrinked <- shrinkThat(deps, ratios, 1000)
rescaled <- shrinkExactlyThat(deps, ratios)
 
plot(deps, col=gray(.5), border="white") # base
plot(shrinked, col=rgb(0,.4,0), border=NA,add=T)
plot(rescaled, col="green", border=NA, add=T)
 
 
writeOGR(shrinked, "OUT", "shrinked", "ESRI Shapefile")
writeOGR(rescaled, "OUT", "rescaled", "ESRI Shapefile")