你能帮助找出解决dotsInPolys引发的长度不匹配错误的最佳方法吗？我认为这是因为多边形数据中有 NA 或 NULL 或一些放克，使其太长。下面是重现错误的代码。最终，我想使用 Leaflet 绘制多个种族，但此时我无法生成随机点所需的纬度/纬度。
require(maptools)
require(tidycensus)
person.number.divider <- 1000
census_api_key("ENTER KEY HERE", install = TRUE)
racevars <- c(White = "B02001_002", #"P005003"
Black = "B02001_003", #Black or African American alone
Latinx = "B03001_003"
)
nj.county <- get_acs(geography = "county", #tract
year = 2015,
variables = racevars,
state = "NJ", #county = "Harris County",
geometry = TRUE,
summary_var = "B02001_001")
library(sf)
st_write(nj.county, "nj.county.shp", delete_layer = TRUE)
nj <- rgdal::readOGR(dsn = "nj.county.shp") %>%
spTransform(CRS("+proj=longlat +datum=WGS84"))
nj@data <- nj@data %>%
tidyr::separate(NAME,
sep =",",
into = c("county", "state")) %>%
dplyr::select(estimat,variabl, GEOID, county) %>%
spread(key = variabl, value = estimat) %>%
mutate(county = trimws(county))
black.dots <- dplyr::select(nj@data, Black) / person.number.divider #%>%
black.dots <- dotsInPolys(nj, as.integer(black.dots$Black), f="random")
# Error in dotsInPolys(nj, as.integer(black.dots$Black), f = "random") :
# different lengths
length(nj) # 63 This seems too many, because I believe NJ has 21 counties.
length(black.dots$Black) # 21
这篇文章(关于解决dotsInPolys错误(maptools(的建议(几乎帮助了我，但我看不到如何将其应用于我的情况。
我可以通过删除黑色流行音乐大于 0 的 NA 和县来更改 nj 空间多边形数据帧的长度，但随后地图不会绘制多个县(也许人口普查下载有问题？
看起来您可能已经弄清楚了这一点，但我想分享另一种使用sf::st_sample()而不是maptools::dotsInPolys()的方法。这样做的一个优点是，您无需将从tidycensus获得的sf对象转换为sp对象。
在下面的示例中，我按种族将人口普查数据拆分为三个sf对象列表，然后对列表的每个元素(每个种族(执行st_sample()。接下来，我将采样点重新组合到一个sf对象中，每个点都有一个新的竞赛变量。最后，我使用tmap来制作地图，尽管您也可以使用ggplot2或leaflet来绘制地图。
library(tidyverse)
library(tidycensus)
library(sf)
library(tmap)
person.number.divider <- 1000
racevars <- c(White = "B02001_002", #"P005003"
Black = "B02001_003", #Black or African American alone
Latinx = "B03001_003"
)
# get acs data with geography in "tidy" form
nj.county <- get_acs(geography = "county", #tract
year = 2015,
variables = racevars,
state = "NJ", #county = "Harris County",
geometry = TRUE,
summary_var = "B02001_001"
)
# split by race
county.split <- nj.county %>%
split(.$variable)
# randomly sample points in polygons based on population
points.list <- map(county.split, ~ st_sample(., .$estimate / person.number.divider))
# combine points into sf collections and add race variable
points <- imap(points.list, ~ st_sf(tibble(race = rep(.y, length(.x))), geometry = .x)) %>%
reduce(rbind)
# map!
tm_shape(nj.county) +
tm_borders(col = "darkgray", lwd = 0.5) +
tm_shape(points) +
tm_dots(col = "race", size = 0.01, pal = "Set2")
我没有足够的代表直接发布地图图像，但它就在这里。