关于解决 dotsInPolys 错误的建议 (maptools)

Question

我是 R 的新手，并且仍在学习解决我遇到的问题的一些方法。我运行陷入了困境，想知道是否有人有建议。

我正在尝试构建点密度图，但运行遇到了 dotsInPolys 函数的错误。该行：

scc.rand <- dotsInPolys(sccpolys, as.integer(plotvar), f="random")

这给了我错误：

> sccdots.rand <- dotsInPolys(sccpolys, as.integer(plotvar), f="random")
Error in dotsInPolys(sccpolys, as.integer(plotvar), f = "random") : 
  different lengths

sccpolys 和 plotvar 的 documentation indicates 需要相同的长度，但我不确定如何仔细检查，或者更重要的是，更正问题。有人对我如何检查问题有什么建议吗？提前致谢。

这是我正在处理的整套代码：

library(maptools)

# Population data
sccpop <- read.csv("nhgis0010_ds98_1970_tract.csv", stringsAsFactors = FALSE)
sccpop.sub <- sccpop[sccpop$COUNTY=="Santa Clara",c(1,3,21,22,23)]

# Shapefile for Census tracts
scctract.shp <- readShapePoly("1970-ca-tracts.shp")
sccpolys <- SpatialPolygonsDataFrame(scctract.shp, data=as(scctract.shp, "data.frame"))

# Merge datasets
sccdata <- merge(sccpolys@data, sccpop.sub, sort=FALSE)
plotvar <- sccdata$C0X001 / 1000 # one dot per 1,000 people
head(sccpolys@data)
head(sccpop.sub)

# Generate random dots in polygons
sccdots.rand <- dotsInPolys(sccpolys, as.integer(plotvar), f="random")

# County boundaries
baycounties.shp <- readShapePoly("ca-counties-1970.shp")
baycounties <- SpatialPolygonsDataFrame(baycounties.shp, data=as(baycounties.shp, "data.frame"))

par(mar=c(0,0,0,0))
plot(baycounties, lwd=0.1)

# Add dots
plot(sccdots.rand, add=TRUE, pch=19, cex=0.1, col="#00880030")

Answer 1

@LincolnMullen 是对的。合并后你有：

> length(sccpolys)
[1] 787

和

> length(plotvar)
[1] 210

为了解决这个问题，您可以替换

sccdots.rand <- dotsInPolys(sccpolys, as.integer(plotvar), f="random")

和

sccdots.rand <- dotsInPolys(sccpolys[sccpolys$GISJOIN %in% sccdata$GISJOIN,], as.integer(plotvar), f="random")

Answer 2

问题是您拥有的区域（即 shapefile 中的多边形）多于您希望绘制的区域。圣克拉拉有 787 个传单，而只有 210 个传单。此外，SpatialPolygonsDataFrame 中的 @data 插槽有一些不必要的操作。这是清理合并代码的解决方案。

library(maptools)

shp <- readShapePoly("1970-ca-tracts.shp")
sccpop <- read.csv("nhgis0010_ds98_1970_tract.csv", stringsAsFactors = FALSE)
sccpop.sub <- sccpop[sccpop$COUNTY=="Santa Clara",c(1,3,21,22,23)]

shp <- merge(shp, sccpop.sub)

现在我们有了 SpatialPolygonsDataFrame 数据，但所有非圣克拉拉县都缺少值。我们想像你上面那样改变人口。下面的第一行进行转换，向数据框添加一列。最简单的方法是将其保存在数据框中而不是作为外部变量。第二行过滤掉所有与人口无关的多边形，即非圣克拉拉县。

shp@data$plotvar <- as.integer(shp@data$C0X001 / 1000)
shp <- shp[!is.na(shp@data$plotvar), ]

现在我们可以像您之前那样进行了。

sccdots.rand <- dotsInPolys(shp, shp@data$plotvar, f="random")

baycounties.shp <- readShapePoly("ca-counties-1970.shp")

par(mar=c(0,0,0,0))
plot(baycounties.shp, lwd=0.1)
plot(sccdots.rand, add=TRUE, pch=19, cex=0.1, col="#00880030")

FWIW，我使用 rgdal::readOGR() 加载 shapefile 获得了更好的结果，但 maptools::readShapePoly() 在这里工作正常。

关于解决 dotsInPolys 错误的建议 (maptools)

Advice on troubleshooting dotsInPolys error (maptools)

r

rstudio