Come combinare gli shapefile a livello di stato dall'ufficio del censimento degli Stati Uniti in una forma nazionale

Question

L'ufficio di censimento non fornisce uno shapefile nazionale di aree di microdati di uso pubblico (la più piccola geografia disponibile nella American Community Survey). Ho provato a combinarli tutti con alcuni metodi diversi, ma anche quello che de-duplica gli identificatori si interrompe una volta che colpisce la California. Sto facendo qualcosa di sciocco o questo richiede una soluzione difficile? Ecco il codice per riprodurre fino al punto in cui le cose si rompono.

library(taRifx.geo) 
library(maptools) 

td <- tempdir() ; tf <- tempfile() 
setInternet2(TRUE) 
download.file("ftp://ftp2.census.gov/geo/tiger/TIGER2014/PUMA/" , tf) 

al <- readLines(tf) 
tl <- al[ grep("geo/tiger/TIGER2014/PUMA/tl_2014_" , al) ] 
fp <- gsub("(.*)geo/tiger/TIGER2014/PUMA/tl_2014_([0-9]*)_puma10\\.zip(.*)" , "\\2" , tl) 

# get rid of alaska 
fp <- fp[ fp != '02' ] 

af <- paste0("ftp://ftp2.census.gov/geo/tiger/TIGER2014/PUMA/tl_2014_" , fp , "_puma10.zip") 

d <- NULL 
for (i in af){ 
    try(file.remove(z) , silent = TRUE) 
    download.file(i , tf , mode = 'wb') 
    z <- unzip(tf , exdir = td) 
    b <- readShapePoly(z[ grep('shp$' , z) ]) 
    if (is.null(d)) d <- b else d <- taRifx.geo:::rbind.SpatialPolygonsDataFrame(d , b , fix.duplicated.IDs = TRUE) 
} 

# Error in `row.names<-.data.frame`(`*tmp*`, value = c("d.0", "d.1", "d.2", : 
    # duplicate 'row.names' are not allowed 
# In addition: Warning message: 
# non-unique values when setting 'row.names': ‘d.0’, ‘d.1’, ‘d.10’, ‘d.11’, ‘d.12’, ‘d.13’, ‘d.14’, ‘d.15’, ‘d.16’, ‘d.17’, ‘d.18’, ‘d.19’, ‘d.2’, ‘d.3’, ‘d.4’, ‘d.5’, ‘d.6’, ‘d.7’, ‘d.8’, ‘d.9’ 

Answer 1

Ecco un altro approccio, che include una scorciatoia per ottenere l'elenco delle directory FTP. Come menzionato @Pop, la chiave è garantire che gli ID siano tutti unici.

library(RCurl) 
library(rgdal) 

# get the directory listing 
u <- 'ftp://ftp2.census.gov/geo/tiger/TIGER2014/PUMA/' 
f <- paste0(u, strsplit(getURL(u, ftp.use.epsv = FALSE, ftplistonly = TRUE), 
         '\\s+')[[1]]) 

# download and extract to tempdir/shps 
invisible(sapply(f, function(x) { 
    path <- file.path(tempdir(), basename(x)) 
    download.file(x, destfile=path, mode = 'wb') 
    unzip(path, exdir=file.path(tempdir(), 'shps')) 
})) 

# read in all shps, and prepend shapefile name to IDs 
shps <- lapply(sub('\\.zip', '', basename(f)), function(x) { 
    shp <- readOGR(file.path(tempdir(), 'shps'), x) 
    shp <- spChFIDs(shp, paste0(x, '_', sapply(slot(shp, "polygons"), slot, "ID"))) 
    shp 
}) 

# rbind to a single object 
shp <- do.call(rbind, as.list(shps)) 

# plot (note: clipping to contiguous states for display purposes) 
plot(shp, axes=T, xlim=c(-130, -60), ylim=c(20, 50), las=1) 

# write out to wd/USA.shp 
writeOGR(shp, '.', 'USA', 'ESRI Shapefile') 

Answer 2

Il tuo problema, come si dovrebbe avere indovinato è dovuto al fatto che ci sono gli ID dei poligoni duplicati nel vostro oggetto d.

Infatti, tutti gli ID poligono nei file "shp" sono "0". Pertanto, hai utilizzato fix.duplicated.IDs = TRUE per renderli diversi.

Questo è strano perché lo taRifx.geo:::rbind.SpatialPolygonsDataFrame avrebbe dovuto risolvere il problema impostando fix.duplicated.IDs = TRUE. Più precisamente, le informazioni vengono trasmesse a sp::rbind.SpatialPolygons che chiama la funzione "interna" sp:::makeUniqueIDs, che utilizza infine la funzione base::make.unique.

Non volevo vedere cosa è andato storto in questa catena. In alternativa, , ti consiglio di configurare gli ID dei tuoi poligoni anziché utilizzare l'opzione fix.duplicated.IDs.

Per risolvere il problema da soli, sostituire il per-loop dal codice seguente:

d <- NULL 
count <- 0 
for (i in af){ 
    try(file.remove(z) , silent = TRUE) 
    download.file(i , tf , mode = 'wb') 
    z <- unzip(tf , exdir = td) 
    b <- readShapePoly(z[ grep('shp$' , z) ]) 

    for (j in 1:length(b@polygons)) 
     b@polygons[[j]]@ID <- as.character(j + count) 
    count <- count + length(b@polygons) 

    if (is.null(d)) 
     d <- b 
    else 
     d <- taRifx.geo:::rbind.SpatialPolygonsDataFrame(d , b) 
} 

La semplice ciclo for su j cambia solo l'ID di ogni poligono nell'oggetto b prima di aspettare che a d.