Dissolved polygons using R not plotting correctly

Question

I'm trying to perform a dissolve in R. I've previously done this in QGIS but I want to achieve this in R to integrate with the rest of my workflow if possible.

I have an ESRI shapefile with small geographical polygons (output areas, if you're familiar with UK census geography). I also have a lookup table provided to me with a list of all OA codes with their associated aggregated geography code.

I can't provide the actual files I'm working on, but comparable files and a minimal reproducable example below:

https://www.dropbox.com/s/4puoof8u5btigxq/oa-soa.csv?dl=1 (130kb csv)
https://www.dropbox.com/s/xqbi7ub2122q14r/soa.zip?dl=1 (~4MB shp)

And code:

require("rgdal")  # for readOGR
require("rgeos")  # for gUnion
require("maptools")

unzip("soa.zip")
soa <- readOGR(dsn = "soa", "england_oac_2011")
proj4string(soa) <- CRS("+init=epsg:27700")  # British National Grid

lookup <- read.csv("oa-soa.csv")

slsoa <- gUnaryUnion(soa, id = lookup$LSOA11CD)

I've also tried:

slsoa <- unionSpatialPolygons(soa, lookup$$LSOA11CD)

but my understanding is that since I have (R)GEOS installed this uses the gUnion methods from the rgeos package anyway.

So, my problem is that the dissolve appears to work; I don't get an error message and the length() function suggests I now have fewer polygons:

length(soa@polygons)    # 1,817
length(slsoa@polygons)  # should be 338

but the plots appear to be the same (i.e. the internal dissolves haven't worked), as demonstrated by the following two plots:

plot(soa)
plot(slsoa)

I've looked around on the internet and stackoverflow to see if I can solve my issue and found several articles but without success.

problems when unioning and dissolving polygons in R (I don't think the quality of the shapefile is the problem because I'm using a lookup table to match geographies).
https://www.nceas.ucsb.edu/scicomp/usecases/PolygonDissolveOperationsR (uses two sp objects, not lookup table).
https://gis.stackexchange.com/questions/93441/problem-with-merging-and-dissolving-a-shapefile-in-r (as far as I can tell I've followed the relevant steps)

Does anyone have any idea what I'm doing wrong and why the plots aren't working correctly?

Thanks muchly.

jlhoward jlhoward · Accepted Answer · 2014-11-23T03:20:04

First, your soa shapefile has 1817 elements, each with a unique code (corresponding to lookup$OA11CD). But your lookup file has only 1667 rows. Obviously, lookup does not have "a list of all OA codes".

Second, unless lookup has the same codes as your shapefile in the same order, using gUnaryUnion(...) this way will yield garbage. You need to merge soa@data with lookup on the corresponding fields first.

Third, gUnaryUnion(...) cannot remove internal boundaries if the polygons are not contiguous (obviously).

This seems to work

soa   <- merge(soa,lookup,by.x="code",by.y="OA11CD",all.x=TRUE)
slsoa <- gUnaryUnion(soa,id=soa$LSOA11CD)
length(slsoa)
# [1] 338
par(mfrow=c(1,2),mar=c(0,0,0,0))
plot(soa);plot(slsoa)

Dissolved polygons using R not plotting correctly

1 Answers