R code for cartograms (see Dougenik et al. 1985)

cartodougenik.R

library(tmap)
library(sp)
library(rgeos)
library(maptools)


# Cartogram 
# algorithm from  Dougenik, Chrisman, Niemeyer (1985): An Algorithm To Construct Continuous Area Cartograms. In: Professional Geographer, 37(1), 75-81.
cartogram <- function(shp, weight, itermax=15, maxSizeError=1.0001) {

  # keep row names
  rown <- rownames(shp@data)
  # identify multi polygons
  multipol <- sapply(seq_len(nrow(shp)), function(i) length(shp@polygons[[i]]@Polygons))
  ## save boundaries in lists  
  tmpcoords <- lapply(seq_len(nrow(shp)), function(i) lapply(seq_len(multipol[i]), function(j) shp@polygons[[i]]@Polygons[[j]]@coords))

  # sum up total value
  value <- shp@data[,weight]
  
  if(any(is.na(value)))
     stop("NA not allowed in weight vector")

  valueTotal <- sum(value, na.rm=T)

  # set meanSizeError
  meanSizeError <- 100

  shp.iter <- shp

  # iterate until itermax is reached
  for(z in 1:itermax) {
    # break if mean Sizer Error is less than maxSizeError
    if(meanSizeError < maxSizeError) break

    # polygon centroids (centroids for multipart polygons)
    centroids <- coordinates(gCentroid(shp.iter, byid=T))

    # area for polygons and total area
    area <- gArea(shp.iter, byid=T)
    area[area <0 ] <- 0
    areaTotal <- gArea(shp.iter)

    # prepare force field calculations
    desired <- areaTotal*value/valueTotal
    radius <- sqrt(area/pi)
    mass <- sqrt(desired/pi) - sqrt(area/pi)
    sizeError <- apply(cbind(area,desired), 1, max)/apply(cbind(area,desired), 1, min)
    meanSizeError <- mean(sizeError, na.rm=T)
    forceReductionFactor <- 1/(1+meanSizeError)

    message(paste0("Mean size error for iteration ", z ,": ", meanSizeError))



    for(i in seq_along(shp.iter)) {
      for(k in seq_len(multipol[i])) {
      newpts <- shp.iter@polygons[[i]]@Polygons[[k]]@coords

      #distance 
      for(j in  seq_len(nrow(centroids))) {
        
        # distance to centroid j        
        distance <- spDistsN1(newpts, centroids[j,])

        # calculate force vector        
        Fij <- mass[j] * radius[j] / distance
        Fbij <- mass[j] * (distance/radius[j])^2 * (4 - 3*(distance/radius[j]))
        Fij[distance <= radius[j]] <- Fbij[distance <= radius[j]]
        Fij <- Fij * forceReductionFactor / distance

        # calculate new border coordinates
        newpts <- newpts + cbind(Fij,Fij) * (newpts - centroids[rep(j,nrow(newpts)),])    

      }

      # save final coordinates from this iteration to coordinate list
      tmpcoords[[i]][[k]] <- newpts
     }
    }
    
    # construct sp-object for area and centroid calculation
    shp.iter <- SpatialPolygons(lapply(seq_along(tmpcoords), function(x) checkPolygonsHoles(Polygons(lapply(tmpcoords[[x]], Polygon), rown[x]))),
                          proj4string = CRS(proj4string(shp)))

  }

  # construct final shape  
  shp.carto <- SpatialPolygons(lapply(seq_along(tmpcoords), function(x) checkPolygonsHoles(Polygons(lapply(tmpcoords[[x]], Polygon),rown[x]))),
                          proj4string = CRS(proj4string(shp)))

  # add data
  shp.cartodf <- SpatialPolygonsDataFrame(shp.carto, shp@data)
  return(shp.cartodf)
}

# Example
data(NLD_prov)
NLD_prov.carto <- cartogram(NLD_prov, "population", 5)
tm_shape(NLD_prov) + tm_borders(gray(0.8)) + tm_shape(NLD_prov.carto) + tm_borders()


data(Europe)
eu <- gSimplify(Europe, 0.1)
eu <- SpatialPolygonsDataFrame(eu, Europe@data, match.ID=F)

eu$pop_est[is.na(eu$pop_est)] <- quantile(eu$pop_est, probs=0.25, na.rm=T)
eu.carto <- cartogram(eu, "pop_est", 15)
tm_shape(eu) + tm_borders(gray(0.8)) + tm_shape(eu.carto) + tm_borders()

library(lineprof)
l <- lineprof(cartogram(NLD_prov, "population", 5))