使用长/纬度坐标查找相应的人口普查区

Question

我正在尝试使用长/纬度坐标来查找他们所属的人口普查区域。这就是我的数据目前的样子。

Latitude    Longitude
38.040357   -84.490224
39.8366017  -86.2430572
38.1952533  -84.3974234
38.1485095  -84.5342804
38.1485095  -84.5342804
...

我需要使用什么类型的包、实际代码或其他任何东西来获得我正在寻找的结果？我有 tigris 包但是不知道怎么用

Answer 1

这里是

sf

和一些坐标的快速介绍。为了好玩，我将使用美国人口普查局的一些制图边界文件，可在 https://www.census.gov/geographies/mapping-files/time-series/geo/carto-boundary-file 找到.html。我下载了“States”，特别是

cb_2018_us_state_500k.zip

并解压到临时目录。

### shapefile
US <- sf::read_sf("~/tmp/cb/cb_2018_us_state_500k.shp")
### your data
dat <- structure(list(Latitude = c(38.040357, 39.8366017, 38.1952533, 38.1485095, 38.1485095), Longitude = c(-84.490224, -86.2430572, -84.3974234, -84.5342804, -84.5342804)), class = "data.frame", row.names = c(NA, -5L))

我们需要将您的数据转换为

sf

对象并应用适当的 CRS。该 shapefile 中的数据是在 NAD83 上定义的。

sf::st_crs(US)
# Coordinate Reference System:
#   User input: NAD83 
#   wkt:
# GEOGCRS["NAD83",
#     DATUM["North American Datum 1983",
#         ELLIPSOID["GRS 1980",6378137,298.257222101,
#             LENGTHUNIT["metre",1]]],
#     PRIMEM["Greenwich",0,
#         ANGLEUNIT["degree",0.0174532925199433]],
#     CS[ellipsoidal,2],
#         AXIS["latitude",north,
#             ORDER[1],
#             ANGLEUNIT["degree",0.0174532925199433]],
#         AXIS["longitude",east,
#             ORDER[2],
#             ANGLEUNIT["degree",0.0174532925199433]],
#     ID["EPSG",4269]]

如果您不知道“CRS”是什么意思，那么您可能会遇到很多问题，因为……更改坐标基准可能很小，或者可能意味着跨越边界。如果坐标是通过某种现代收集设备导出的，那么可以安全地假设您的坐标是 WGS84 (EPSG:4326)。为了这个答案，我会盲目地假设这一点，但请尽职调查并核实。如果您不这样做，人口普查结果可能会默默地错误。

sf::st_set_crs(US, 4326)
datsf <- sf::st_as_sf(dat, coords=c("Longitude","Latitude"), crs = 4326L)

我们可以展示情节只是为了好玩。我会使用

ggplot2

，但这不是唯一的方法。

library(ggplot2)
ggplot(US) +
  geom_sf() +
  geom_sf(data = datsf, color = "red")

让我们检查一下您的积分落在哪里。

ind <- sf::st_intersects(datsf, US, sparse=T)
ind
# Sparse geometry binary predicate list of length 5, where the predicate was `intersects'
#  1: 48
#  2: 54
#  3: 48
#  4: 48
#  5: 48

我们可以通过

更好地形象化这一点

ggplot(US[unique(unlist(ind)),]) +
  geom_sf() +
  geom_sf(data = datsf, color = "red")

那么你的观点所在的区域是什么呢？第 1 行和第 3-5 行位于“48”（肯塔基州），第 2 行位于“54”（印第安纳州），

US[unique(unlist(ind)),]
# # A tibble: 2 × 10
#   STATEFP STATENS  AFFGEOID    GEOID STUSPS NAME     LSAD         ALAND     AWATER                                                                    geometry
#   <chr>   <chr>    <chr>       <chr> <chr>  <chr>    <chr>        <dbl>      <dbl>                                                          <MULTIPOLYGON [°]>
# 1 21      01779786 0400000US21 21    KY     Kentucky 00    102279490672 2375337755 (((-89.40565 36.52817, -89.39869 36.54233, -89.39681 36.55197, -89.39663 3…
# 2 18      00448508 0400000US18 18    IN     Indiana  00     92789302676 1538002829 (((-88.09776 37.90403, -88.09448 37.90596, -88.08624 37.90571, -88.07145 3…

你也许可以用它来形式化

dat$State <- US$NAME[unlist(ind)]
dat
#   Latitude Longitude    State
# 1 38.04036 -84.49022 Kentucky
# 2 39.83660 -86.24306  Indiana
# 3 38.19525 -84.39742 Kentucky
# 4 38.14851 -84.53428 Kentucky
# 5 38.14851 -84.53428 Kentucky