使用 tigris 从 Lat/Lon 获取人口普查区
Get Census Tract from Lat/Lon using tigris
我有相对较多的坐标想要获取人口普查区(除了 FIPS 代码之外)。我知道我可以使用 call_geolocator_latlon
查找单个 lat/lon 对(就像 所做的那样),但这对我的目的来说似乎不切实际,因为该函数向人口普查局发出单一调用 API,我想在我的 ~200,000 双鞋上 运行 需要很长时间。
是否有更快的方法来执行此操作,也许是通过使用 block_groups
函数下载每个州的 shapefile 并从那里从 lat/lon 映射到人口普查区?
这不使用 tigris
,但使用 sf::st_within()
检查重叠区域点的数据框。
我在这里使用 tidycensus
将加利福尼亚地区的地图放入 R。
library(sf)
ca <- tidycensus::get_acs(state = "CA", geography = "tract",
variables = "B19013_001", geometry = TRUE)
现在模拟一些数据:
bbox <- st_bbox(ca)
my_points <- data.frame(
x = runif(100, bbox[1], bbox[3]),
y = runif(100, bbox[2], bbox[4])
) %>%
# convert the points to same CRS
st_as_sf(coords = c("x", "y"),
crs = st_crs(ca))
我在这里做 100 点是为了能够 ggplot()
结果,但是 1e6 的重叠计算很快,在我的笔记本电脑上只需要几秒钟。
my_points$tract <- as.numeric(st_within(my_points, ca)) # this is fast for 1e6 points
结果:
head(my_points) # tract is the row-index for overlapping census tract record in 'ca'
# but part would take forever with 1e6 points
library(ggplot2)
ggplot(ca) +
geom_sf() +
geom_sf(data = my_points, aes(color = is.na(tract)))
上面的答案很好。要获取人口普查区 ID,您还可以使用 st_join()
。区域 ID 的 NA 是位于加利福尼亚边界框内但不与州本身相交的那些点。
library(tigris)
library(tidyverse)
library(sf)
ca_tracts <- tracts("CA", class = "sf") %>%
select(GEOID, TRACTCE)
bbox <- st_bbox(ca_tracts)
my_points <- data.frame(
x = runif(200000, bbox[1], bbox[3]),
y = runif(200000, bbox[2], bbox[4])
) %>%
# convert the points to same CRS
st_as_sf(coords = c("x", "y"),
crs = st_crs(ca_tracts))
my_points_tract <- st_join(my_points, ca_tracts)
> my_points_tract
Simple feature collection with 200000 features and 2 fields
geometry type: POINT
dimension: XY
bbox: xmin: -124.4819 ymin: 32.52888 xmax: -114.1312 ymax: 42.0095
epsg (SRID): 4269
proj4string: +proj=longlat +datum=NAD83 +no_defs
First 10 features:
GEOID TRACTCE geometry
1 06025012400 012400 POINT (-114.6916 33.42711)
2 <NA> <NA> POINT (-118.4255 41.81896)
3 06053990000 990000 POINT (-121.8154 36.22736)
4 06045010200 010200 POINT (-123.6909 39.70572)
5 <NA> <NA> POINT (-116.9055 37.93532)
6 06019006405 006405 POINT (-119.511 37.09383)
7 06049000300 000300 POINT (-120.7215 41.3392)
8 <NA> <NA> POINT (-115.8916 39.32392)
9 06023990100 990100 POINT (-124.2737 40.14106)
10 06071008901 008901 POINT (-117.319 35.62759)
我有相对较多的坐标想要获取人口普查区(除了 FIPS 代码之外)。我知道我可以使用 call_geolocator_latlon
查找单个 lat/lon 对(就像
是否有更快的方法来执行此操作,也许是通过使用 block_groups
函数下载每个州的 shapefile 并从那里从 lat/lon 映射到人口普查区?
这不使用 tigris
,但使用 sf::st_within()
检查重叠区域点的数据框。
我在这里使用 tidycensus
将加利福尼亚地区的地图放入 R。
library(sf)
ca <- tidycensus::get_acs(state = "CA", geography = "tract",
variables = "B19013_001", geometry = TRUE)
现在模拟一些数据:
bbox <- st_bbox(ca)
my_points <- data.frame(
x = runif(100, bbox[1], bbox[3]),
y = runif(100, bbox[2], bbox[4])
) %>%
# convert the points to same CRS
st_as_sf(coords = c("x", "y"),
crs = st_crs(ca))
我在这里做 100 点是为了能够 ggplot()
结果,但是 1e6 的重叠计算很快,在我的笔记本电脑上只需要几秒钟。
my_points$tract <- as.numeric(st_within(my_points, ca)) # this is fast for 1e6 points
结果:
head(my_points) # tract is the row-index for overlapping census tract record in 'ca'
# but part would take forever with 1e6 points
library(ggplot2)
ggplot(ca) +
geom_sf() +
geom_sf(data = my_points, aes(color = is.na(tract)))
上面的答案很好。要获取人口普查区 ID,您还可以使用 st_join()
。区域 ID 的 NA 是位于加利福尼亚边界框内但不与州本身相交的那些点。
library(tigris)
library(tidyverse)
library(sf)
ca_tracts <- tracts("CA", class = "sf") %>%
select(GEOID, TRACTCE)
bbox <- st_bbox(ca_tracts)
my_points <- data.frame(
x = runif(200000, bbox[1], bbox[3]),
y = runif(200000, bbox[2], bbox[4])
) %>%
# convert the points to same CRS
st_as_sf(coords = c("x", "y"),
crs = st_crs(ca_tracts))
my_points_tract <- st_join(my_points, ca_tracts)
> my_points_tract
Simple feature collection with 200000 features and 2 fields
geometry type: POINT
dimension: XY
bbox: xmin: -124.4819 ymin: 32.52888 xmax: -114.1312 ymax: 42.0095
epsg (SRID): 4269
proj4string: +proj=longlat +datum=NAD83 +no_defs
First 10 features:
GEOID TRACTCE geometry
1 06025012400 012400 POINT (-114.6916 33.42711)
2 <NA> <NA> POINT (-118.4255 41.81896)
3 06053990000 990000 POINT (-121.8154 36.22736)
4 06045010200 010200 POINT (-123.6909 39.70572)
5 <NA> <NA> POINT (-116.9055 37.93532)
6 06019006405 006405 POINT (-119.511 37.09383)
7 06049000300 000300 POINT (-120.7215 41.3392)
8 <NA> <NA> POINT (-115.8916 39.32392)
9 06023990100 990100 POINT (-124.2737 40.14106)
10 06071008901 008901 POINT (-117.319 35.62759)