基于经纬度为DataFrame添加美国人口普查区（Census Tract）

你已经走对了第一步——用tidycensus获取普查区数据、用sf把经纬度转成点要素，接下来只需要通过空间连接就能把普查区ID匹配到你的数据里，具体步骤如下：

1. 修正普查区数据的坐标系转换 #

注意：不要用st_set_crs()修改已有坐标系的对象——它是给无CRS的数据设置坐标系，不是转换坐标系。tidycensus::get_acs()返回的几何数据默认是EPSG:4269（NAD83坐标系），要转换成和你的点数据一致的WGS84（EPSG:4326），应该用st_transform()，并且必须重新赋值：

library(sf)
library(tidycensus)

# 获取马里兰州普查区数据（带几何信息）
MD <- get_acs(state = "MD", geography = "tract",
              variables = "B19013_001", geometry = TRUE)

# 转换到WGS84坐标系（和你的点数据CRS一致）
MD <- st_transform(MD, crs = 4326)

# 把你的经纬度数据转成sf点要素
datsf <- st_as_sf(EIA_Data, coords = c("Longitude", "Latitude"), crs = 4326)

2. 空间左连接匹配普查区 #

用st_join()做空间左连接，指定join = st_within（判断点是否在普查区多边形内），这样就能把每个点对应的普查区信息关联上：

# 空间左连接，保留原数据所有行
matched_data <- st_join(datsf, MD, join = st_within, left = TRUE)

# 提取普查区ID（GEOID）作为"Tract"列
matched_data$Tract <- matched_data$GEOID

# 如果不需要几何列，可以转成普通data.frame
matched_data_df <- st_drop_geometry(matched_data)

3. 针对大型数据集的优化（可选） #

因为你的数据有大量重复经纬度，先对经纬度去重后再匹配，能大幅提升效率：

# 提取唯一经纬度行
unique_points <- datsf[!duplicated(st_coordinates(datsf)), ]

# 匹配唯一经纬度的普查区
unique_matched <- st_join(unique_points, MD, join = st_within, left = TRUE)
unique_matched$Tract <- unique_matched$GEOID

# 用经纬度把匹配结果合并回原数据（先转普通df）
unique_matched_df <- st_drop_geometry(unique_matched)[, c("Longitude", "Latitude", "Tract")]
final_data <- merge(EIA_Data, unique_matched_df, by = c("Longitude", "Latitude"), all.x = TRUE)

这样就能得到带"Tract"列的最终数据集，所有经纬度对应的普查区ID都会被正确匹配。

内容的提问来源于stack exchange，提问作者user17073706