Получите переписной тракт от широты/долготы, используя tigris

У меня относительно большое количество координат, по которым я хотел бы получить переписной участок (помимо кода ФИПС). Я знаю, что могу искать отдельные пары широта/долгота, используя call_geolocator_latlon (как это сделано здесь ), но для моих целей это кажется нецелесообразным, поскольку функция выполняет один вызов API бюро переписи, и я полагаю, что на мои ~ 200 000 пар уйдет очень много времени.

Есть ли более быстрый способ сделать это, возможно, загрузив шейп-файлы для каждого штата с помощью функции block_groups и оттуда сопоставив широту/долготу с переписным массивом?


r tigris census
person mlinegar    schedule 09.09.2018    source источник
comment
Вы можете проверить cenpy и pysal   -  person akrun    schedule 09.09.2018

Ответы (2)


arrow_upward
4
arrow_downward

Это не использует 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)))

демонстрация карты ca

person Nate    schedule 09.09.2018
comment
Это фантастика! Кроме того, вы случайно не знаете, можно ли получить код FIPS/название округа в том же запросе, что и переписной участок? - person mlinegar; 10.09.2018
comment
tidycensus поставляется со встроенным фреймом данных fips_codes, так что вы можете сделать merge/join по названию округа. - person Nate; 10.09.2018

arrow_upward
3
arrow_downward

Отличный ответ выше. Чтобы получить идентификаторы переписных участков, вы также можете использовать st_join(). 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)
person kwalkertcu    schedule 10.09.2018