使用sf::st_join空间连接时，科罗拉多点集错误匹配至堪萨斯普查区

我尝试用sf::st_join()将点sf对象与科罗拉多州、堪萨斯州的普查区面数据进行空间连接，但通过leaflet确认点确实位于科罗拉多州区域后，执行空间连接的结果却显示所有点都匹配到了堪萨斯州（state_code=20）。

我的代码： #

数据生成与普查区获取：

library(tidycensus)
library(sf)
library(dplyr)
library(tidyverse)

# Set seed for reproducibility
set.seed(42)

# Generate dummy data for points in New York
points <- data.frame(
  longitude = runif(300, min = -109, max = -102),  # Approximate longitude boundaries of Colorado
  latitude = runif(300, min = 36.993076, max = 41)  # Approximate latitude boundaries of Colorado
)

# Print the first few rows of the dummy data
points <- st_as_sf(points, coords = c("longitude", "latitude"), crs = "ESRI:102003")

tract2010 <- get_decennial(geography = "tract", variables = "P001001", year = 2010,
                          state = as.list(c("Colorado", "Kansas")), geometry = TRUE)

tract2010$state_code <- substr(tract2010$GEOID, 1, 2)
table(tract2010$state_code)

# make same CRS
tract2010 <- st_transform(tract2010, st_crs(points))

Leaflet验证点位置：

# test where it is
library(leaflet)
leaflet() %>%
  addTiles() %>%
  addMarkers(data = points)

（配图：点位于科罗拉多州区域）

空间连接代码：

#spatial join
points <- st_join(points, tract2010, join = st_within)
table(points$state_code, useNA = "always")

问题原因与解决方法 #

核心问题：点数据的CRS设置错误

你生成的是WGS84经纬度格式的坐标（longitude/latitude），但却错误地指定了CRS为ESRI:102003（北美阿尔伯斯等积投影，单位为米）。这会导致经纬度数值被当作米单位解析，点的实际位置完全偏移，最终全部落在堪萨斯州对应的投影区域内。

修正步骤：

1. 正确设置点数据的初始CRS：
   因为输入的是经纬度，所以要指定为WGS84（EPSG:4326）：

   points <- st_as_sf(points, coords = c("longitude", "latitude"), crs = "EPSG:4326")
   

2. 统一坐标系：
   将点数据转换到与普查区一致的坐标系（建议使用投影坐标系，避免经纬度的精度问题）：

   # 将点转换为普查区的坐标系（get_decennial默认返回EPSG:4326，也可手动指定投影CRS）
   points <- st_transform(points, st_crs(tract2010))
   # 或者直接指定投影CRS：points <- st_transform(points, "ESRI:102003")
   

修正后的完整代码示例：

library(tidycensus)
library(sf)
library(dplyr)
library(leaflet)

set.seed(42)

# 生成科罗拉多州范围的经纬度点
points <- data.frame(
  longitude = runif(300, min = -109, max = -102),
  latitude = runif(300, min = 36.993076, max = 41)
)

# 正确设置初始CRS为WGS84
points <- st_as_sf(points, coords = c("longitude", "latitude"), crs = "EPSG:4326")

# 获取科罗拉多和堪萨斯的普查区数据
tract2010 <- get_decennial(geography = "tract", variables = "P001001", year = 2010,
                            state = c("Colorado", "Kansas"), geometry = TRUE)

tract2010$state_code <- substr(tract2010$GEOID, 1, 2)

# 统一坐标系
points <- st_transform(points, st_crs(tract2010))

# 执行空间连接
points_joined <- st_join(points, tract2010, join = st_within)
table(points_joined$state_code, useNA = "always")

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