Как совместить шейп-файлы государственного уровня от бюро переписи Соединенных Штатов в общенациональной форме

Question

Бюро переписи не предоставляет общенациональный шейп-файл областей микроданных общественного пользования (самая маленькая география, доступная в Обзоре Американского сообщества). Я попытался объединить их все с помощью нескольких разных методов, но даже тот, который де-дублирует идентификаторы, ломается, когда он попадает в Калифорнию. Я делаю что-то глупое или это требует сложного обходного пути? Вот код, который можно воспроизвести до того момента, когда все сломается.

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

Вот еще один подход, который включает короткое сокращение для получения списка каталогов FTP. Как упоминалось в @Pop, ключ должен гарантировать, что идентификаторы уникальны.

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

Ваша проблема, как вы должны были догадаться это связано с тем, что существуют повторяющиеся идентификаторы полигонов в вашем объекте d.

Действительно, все идентификаторы многоугольника в ваших файлах «shp»: "0". Таким образом, вы использовали fix.duplicated.IDs = TRUE, чтобы сделать их разными.

Это странно, потому что taRifx.geo:::rbind.SpatialPolygonsDataFrame должен был исправить это, когда вы установили fix.duplicated.IDs = TRUE. Точнее, информация передается на sp::rbind.SpatialPolygons, которая вызывает «внутреннюю» функцию sp:::makeUniqueIDs, которая, наконец, использует функцию base::make.unique.

Я не хотел видеть, что пошло не так в этой цепочке. Альтернативно, Я советую вам установить идентификаторы ваших полигонов, вместо использования опции fix.duplicated.IDs.

Чтобы исправить это самостоятельно, замените ваш для цикла по следующему коду:

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) 
} 

Простой для петли на j только изменяет идентификатор каждого многоугольника в объекте b перед выжидать его до d.