如何将美国人口普查局的州级形状文件合并为全国性形状

人口普查局不提供全国范围内公共使用微数据区域的形状文件（美国社区调查中可用的最小地理区域）。我尝试用几种不同的方法将它们结合起来，但即使是消除重复标识符的方法一旦到达加利福尼亚州也会崩溃。我是在做一些愚蠢的事情还是需要一个困难的解决方法？下面是可重现直到出现问题的代码。

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’ 

正如您应该猜到的那样，您的问题是由于对象中存在重复的多边形 IDd.

事实上，“shp”文件中的所有多边形 ID 都是"0"。因此，您使用了fix.duplicated.IDs = TRUE让他们与众不同。

这很奇怪，因为taRifx.geo:::rbind.SpatialPolygonsDataFrame应该按照你的设置修复它fix.duplicated.IDs = TRUE。更准确地说，信息被传送到sp::rbind.SpatialPolygons它调用“内部”函数sp:::makeUniqueIDs，最终使用该函数base::make.unique.

我不想看到这个链条出了什么问题。或者，我建议你自己设置多边形的 ID，而不是使用fix.duplicated.IDs option.

要自行修复此问题，请将 for 循环替换为以下代码：

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

简单的 for 循环j只改变对象中每个多边形的IDb在将其绑定到之前d.