I have a very simple operation I would like to perform: combining two shape files. Specifically, I have census tract shape files for each state in the US that I would like to combine into one shape file. Ultimately, I want to take the combined shape file and perform an overlay on a set of latitude and longitude coordinates, to determine which census tracts my coordinates fall into.
I have seen much discussion on this (Combining bordering shapefiles in R). However, all of the discussion are outdated and I am hoping that the packages have been improved in the meantime.
The files I am using are located here: ftp://ftp2.census.gov/geo/tiger/TIGER2010/TRACT/2010/
The code below can be used to recreate the files I am working with. It requires downloading two files, approximately 11 megabytes total. The run time should be only a minute.
Any help is greatly appreciated. This seems like a trivial operation to perform. Perhaps if I had more experience with geo-spatial mapping I could make better use of the documentation available.
Here are a few things I have tried:
### 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 = '
', 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)
Thank you for your time.
See Question&Answers more detail:
os