Use Rcartogram on a SpatialPolygonsDataFrame object

Question

I'm trying to do the same thing asked in this question, Cartogram + choropleth map in R, but starting from a SpatialPolygonsDataFrame and hoping to end up with the same type of object.

I could save the object as a shapefile, use scapetoad, reopen it and convert back, but I'd rather have it all within R so that the procedure is fully reproducible, and so that I can code dozens of variations automatically.

I've forked the Rcartogram code on github and added my efforts so far here.

Essentially what this demo does is create a SpatialGrid over the map, look up the population density at each point of the grid and convert this to a density matrix in the format required for cartogram() to work on. So far so good.

But, how to interpolate the original map points based on the output of cartogram()?

There are two problems here. The first is to get the map and grid into the same units to allow interpolation. The second is to access every point of every polygon, interpolate it, and keep them all in right order.

The grid is in grid units and the map is in projected units (in the case of the example longlat). Either the grid must be projected into longlat, or the map into grid units. My thought is to make a fake CRS and use this along with the spTransform() function in package(rgdal), since this handles every point in the object with minimal fuss.

Accessing every point is difficult because they are several layers down into the SpPDF object: object>polygons>Polygons>lines>coords I think. Any ideas how to access these while keeping the structure of the overall map intact?

Answer 1

This problem can be solved with the getcartr package, available on Chris Brunsdon's GitHub, as beautifully explicated in this blog post.

The quick.carto function does exactly what you want -- takes a SpatialPolygonsDataFrame as input and has a SpatialPolygonsDataFrame as output.

Reproducing the essence of the example in the blog post here in case the link goes dead, with my own style mixed in & typos fixed:

(Shapefile; World Bank population data)

library(getcartr)
library(maptools)
library(data.table)

world <- readShapePoly("TM_WORLD_BORDERS-0.3.shp")
#I use data.table, see blog post if you want a base approach;
#  data.table wonks may be struck by the following step as seeming odd;
#  see here: http://stackoverflow.com/questions/32380338
#  and here: https://github.com/Rdatatable/data.table/issues/1310
#  for some background on what's going on.
world@data <- setDT(world@data)

world.pop <- fread("sp.pop.totl_Indicator_en_csv_v2.csv",
                   select = c("Country Code", "2013"),
                   col.names = c("ISO3", "pop"))

world@data[world.pop, Population := as.numeric(i.pop), on = "ISO3"]

#calling quick.carto has internal calls to the
#  necessary functions from Rcartogram
world.carto <- quick.carto(world, world$Population, blur = 0)

#plotting with a color scale
x <- world@data[!is.na(Population), log10(Population)]
ramp <- colorRampPalette(c("navy", "deepskyblue"))(21L)
xseq <- seq(from = min(x), to = max(x), length.out = 21L)
#annoying to deal with NAs...
cols <- ramp[sapply(x, function(y)
  if (length(z <- which.min(abs(xseq - y)))) z else NA)]

plot(world.carto, col = cols,
     main = paste0("Cartogram of the World's",
                   " Population by Country (2013)"))