Tengo una operación muy simple que me gustaría realizar: combinar dos archivos de formas. Específicamente, tengo archivos de formas del censo para cada estado de los Estados Unidos que me gustaría combinar en un archivo de formas. En última instancia, quiero tomar el archivo de formas combinadas y realizar una superposición en un conjunto de coordenadas de latitud y longitud, para determinar en qué censos se encuentran mis coordenadas.

He visto mucha discusión sobre esto (Combinando shapefiles limítrofes en R). Sin embargo, toda la discusión está desactualizada y espero que los paquetes se hayan mejorado mientras tanto.

Los archivos que estoy usando se encuentran aquí:ftp://ftp2.census.gov/geo/tiger/TIGER2010/TRACT/2010/

El siguiente código se puede usar para recrear los archivos con los que estoy trabajando. Requiere descargar dos archivos, aproximadamente 11 megabytes en total. El tiempo de ejecución debe ser de sólo un minuto.

Cualquier ayuda es muy apreciada. Esto parece una operación trivial a realizar. Quizás si tuviera más experiencia con la cartografía geoespacial podría hacer un mejor uso de la documentación disponible.

Aquí hay algunas cosas que he intentado:

### Insert your file path here
FPATH <- './data'

### Set up library
require(rgeos)
require(maptools)
require(RCurl)
require(parallel)
cl <- makeCluster(detectCores())

### Download files... (~11,000 KB total for this example)
ftp <- 'ftp://ftp2.census.gov/geo/tiger/TIGER2010/TRACT/2010/'
files <- getURLContent(ftp, dirlistonly = T)
files <- unlist(strsplit(files, split = '\r\n', fixed = T))
files <- grep('2010_[[:digit:]]{2}_', files, value = T)[1:2]  # Only grab two files for this example
clusterMap(cl, download.file, url = paste0(ftp, files), destfile = paste0(FPATH, files))

### Unzip shape files
files <- list.files(FPATH, full.names = T)
clusterMap(cl, unzip, zipfile = files, exdir = FPATH)

### Read in shape files
files <- list.files(FPATH, pattern = "shp$", full.names = T)
a <- readShapePoly(fn = files[1])
b <- readShapePoly(fn = files[2])

### Attempt to join two shape files
spRbind(a, b)  # Error in spRbind(as(obj, "SpatialPolygons"), as(x, "SpatialPolygons")) : non-unique polygon IDs
gUnion(a, b)   # Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, "rgeos_union")   : geos_geospolygon2SpatialPolygons: ng > length(IDs)

Gracias por tu tiempo.