Background

When using the getcartr package you discover that the output cartogram is not easily plotted on top of the original set of polygons for comparison. This problem can be remedied using the 2d Cartesian transformation functions provided by German Carrillo in his vec2dtransf package. I’ll illustrate this by creating and then transforming a cartogram of Irish Administrative County populations.

Getting started

First we load the libraries. We need rgdal for the readOGR() function that reads an ESRI shapefile into a spatial polygons (in our case) data frame. The getcartr library provides the quick.carto() function for creating a Gastner-Newman style cartogram. We also need rgeos since we’ll need to compare the areas in both cartograms to determine that no untoward distortions have taken place as a result of our antics.

library(rgdal)
library(getcartr)
library(vec2dtransf)
library(rgeos)

Load the external shapefile into an SPDF. The DataDir variable points to the location on my computer’s hard disk where the shapefiles are located.

cty <- readOGR(DataDir,"Census2011_Admin_Counties_generalised20m")

Now we compute the cartogram, using the total population in 2011 as the stock variable. Some of our other experiments with the Gastner-Newman algorithm suggest that areas of very high population density are a little too small in the output cartogram, and areas of low density are a little too large. This is particularly true for Dublin City.

cty_c <- quick.carto(cty,cty$Total2011)

The next step is to extract the bounding boxes for the original SPDF, cty.bbox and for the cartogram, cty_c.bbox. You’ll see that the coordinates are rather different.

cty.bbox <- bbox(cty)
cty_c.bbox <- bbox(cty_c)
print(cty.bbox)
       min      max
x 17491.14 334558.6
y 19589.93 466919.3
print(cty_c.bbox)
        min      max
x 20.862654 114.8385
y  3.781609 123.5291

Co-location details

We can use these to set an an affine transformation to map the locations of the cartogram bounding box corners onto the corresponding locations in the original SPDF This requires the use of an affine transformation matrix which will allow us to carry out the requisite scaling and transformation. An affine transformation will also deal with rotation of planar coordinates, although this will not be needed in this example.

The corresponding control point coordinates are organised a data frame; those for the origin system in the first two columns, and those for the destination system in the next two columns. Each row represents a separate control point.

xfrom <- c(cty_c.bbox[1,1],cty_c.bbox[1,2],cty_c.bbox[1,1],cty_c.bbox[1,2])
yfrom <- c(cty_c.bbox[2,1],cty_c.bbox[2,2],cty_c.bbox[2,2],cty_c.bbox[2,1])
Xto   <- c(cty.bbox[1,1],  cty.bbox[1,2],  cty.bbox[1,1],  cty.bbox[1,2])
Yto   <- c(cty.bbox[2,1],  cty.bbox[2,2],  cty.bbox[2,2],  cty.bbox[2,1])
CPs <- data.frame(xfrom,yfrom,Xto,Yto)
print(CPs)
      xfrom      yfrom       Xto       Yto
1  20.86265   3.781609  17491.14  19589.93
2 114.83846 123.529078 334558.59 466919.31
3  20.86265 123.529078  17491.14 466919.31
4 114.83846   3.781609 334558.59  19589.93

An affine transformation matrix contains the coefficients of the following equations.

X = ax + by + c
Y = dx + ey + f

… where (x,y) are the coordinates of a point in the origin coordinate system (in our case, the cartogram), and (X,Y) are the coordinates of the corresponding point in the destination coordinate system (that of the original SPDF).

Three steps are required: (i) create an object of class AffineTransformation from the control points, (ii) compute the parameters of the transformation matrix and (iii) apply the transformation to the cartogram. We can print the parameters, either with their identifying coefficient letters or as a matrix. The first line scales and translates the eastings, and the second line scale and translates the northings.

atM <- AffineTransformation(controlPoints=CPs)
calculateParameters(atM)
cty_ct <- applyTransformation(atM, cty_c)
getParameters(atM)
            a             b             c             d             e 
 3.373926e+03  2.916521e-13 -5.289792e+04  5.648917e-13  3.735606e+03 
            f 
 5.463328e+03 
print(round(matrix(getParameters(atM),2,3,byrow=TRUE),3))
         [,1]     [,2]       [,3]
[1,] 3373.926    0.000 -52897.916
[2,]    0.000 3735.606   5463.328

There is a slight difference in the scaling of the eastings and northings. The eastings are scaled by a factor of 3373.926 and the northings by 3735.606. This results in an east-west enlargement of slightly over 10%. The translation of the eastings is -52897.92m, and that of the northings 5463.33m. Notice that the other values, corresponding c and f, are almost zero: there is no rotation to be applied to the coordinates.

And if we plot the scaled and translated cartogram onto the original county boundaries we can see that there’s a reasonable correspondence - we also get some idea of the rescaling that’s been done as a results of the application of the cartogram algorithm.

plot(cty,   border="lightgrey")
plot(cty_ct,border="black",add=T)
box()

Confirmation

As we have applied slightly different scalings to each of the x and y coordinates, the raises the question of whether this any distortion into relative areas of the polygons. If we compute the areas of the cartogram from quick.carto() and the transformed cartogram from applyTransformation() and plot them together, all the points should lie on a straight line which passes through the origin.

N <- length(cty)
cty_cAreas  <- rep(0,N)
cty_ctAreas <- rep(0,N)
for (i in 1:N) {
   cty_cAreas[i]  <- gArea(cty_c[i,])
   cty_ctAreas[i] <- gArea(cty_ct[i,])
}
plot(cty_cAreas,cty_ctAreas,pch=16,cex=0.75,
     xlab="Original Cartogram",ylab="Scaled & Translated Cartogram",
     xaxt="n",yaxt="n",main="Polygon Area Comparison")
axis(1,cex.axis=0.75)
axis(2,cex.axis=0.75)
abline(lm(cty_ctAreas~cty_cAreas+0),lty=2)

… which rather reassuringly, they do.

A feature of the getcartr library is the provision of the carto.transform() function. This returns a function, say to.carto(), which maps data in the original coordinate system into the cartogram coordinate system. This allows additional topographic data such as a road network, or river network to be transformed into the cartogram space.

To apply the affine transformation outlined in this note, you create the cartogram boundaries first, and compute its affine transformation matrix. Then you can use this in AffineTransformation() on each of the additional outputs from to.carto().

All this can be put together as a function.

carto.relocate <- function(cart,dest){
   dest.bbox <- bbox(dest)
   cart.bbox <- bbox(cart)
   xfrom <- c(cart.bbox[1,1],cart.bbox[1,2],cart.bbox[1,1],cart.bbox[1,2])
   yfrom <- c(cart.bbox[2,1],cart.bbox[2,2],cart.bbox[2,2],cart.bbox[2,1])
   Xto   <- c(dest.bbox[1,1],dest.bbox[1,2],dest.bbox[1,1],dest.bbox[1,2])
   Yto   <- c(dest.bbox[2,1],dest.bbox[2,2],dest.bbox[2,2],dest.bbox[2,1])
   CPs <- data.frame(xfrom,yfrom,Xto,Yto)
   atM <- AffineTransformation(controlPoints=CPs)
   calculateParameters(atM)
   return(applyTransformation(atM, cart))
}

new_cart <- carto.relocate(cty_c, cty)
plot(cty,border="lightgrey")
plot(new_cart,border="black",add=T)
box()