Spatial operations: dissolve and overlay

snippets.R

setwd("/home/projects")

library(data.table)
library(rgdal)
library(rgeos)
library(maptools)
library(stringr)

# Load GAUL 2008 (release 2009) adm-2 from Darby
g <- readOGR("./maps/admin", "g2008_2")

# Dissolve multi-feature districts
# Note that this might kill the server, maybe best to run continent by continent or country by country
g.merged <- unionSpatialPolygons(g, ID=g$ADM2_CODE)

# Merge back all attributes from [g] into [g.merged]
g <- data.table(g@data)
setkey(g, ADM2_CODE)
g <- unique(g)
g.merged <- SpatialPolygonsDataFrame(g.merged, data.frame(g), match.ID=FALSE)

# Load all CRP layers into a list
f <- list.files("./2013-CRP/maps", glob2rx("*.shp"), recursive=TRUE)
l <- data.frame(dir=dirname(f), base=str_replace(basename(f), ".shp", ""), stringsAsFactors=F)
f <- lapply(1:dim(l)[1], function(i) readOGR(paste0("./2013-CRP/maps/", l[i, "dir"]), l[i, "base"]))

# Make sure the coordinate systems are identical before overlaying
# If they don't match then reproject the layers
sapply(f, slot, "proj4string")
# [[1]]
# CRS arguments:
#     +proj=longlat +a=6378137 +b=6378137 +no_defs

# If no projection is available, then assume same as others or best guess, e.g.
f[[4]]@proj4string <- f[[1]]@proj4string

g.merged@proj4string
# CRS arguments:
#     +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

# They don't match, so reproject all CRP layers to match GAUL
f <- lapply(f, spTransform, g.merged@proj4string)

# Define generic overlay function (att is a vector of CRP attributes we want to keep)
# e.g. att <- c("title", "description", "resolution", etc..)
gaulOverlay <- function(x, att) {
  tmp <- over(x, g.merged)
  tmp <- data.table(featureID=row.names(x),
      g.merged@data[tmp, c("ADM0_CODE", "ADM0_NAME", "ADM1_CODE", "ADM1_NAME", "ADM2_CODE", "ADM2_NAME")])
  tmp[, att] <- x@data[, att]
  return(tmp)  
}

# Overlay CRP layer by CRP layer
att <- c("title", "description", "resolution", "source")
out <- lapply(f, gaulOverlay, att)

# Combine all
out <- do.call(rbind, out)

# For a quick visual check
plot(g.merged[g.merged$ADM2_CODE %in% out[[1]]$ADM2_CODE,], col="red")
plot(f[[1]][out[[1]]$featureID,], col="#D9D9D920", add=TRUE)
# This shows that all intersected districts are returned by over() even if only a tiny area intersects
# We can tweak that behavior to discard small sleeves

# Or check interactively through Github's new geojson feature
writeOGR(g.merged[g.merged$ADM2_CODE %in% out[[1]]$ADM2_CODE,],
    dsn="./test.geojson", layer="lyr1", driver="GeoJSON")
writeOGR(f[[1]][out[[1]]$featureID,],
    dsn="./test.geojson", layer="lyr2", driver="GeoJSON")
# View on github => https://github.com/mbacou/config/blob/master/test.geojson
# (only showing Ethiopia districts in this example)