Point-in-Polygon 但坐标系关闭，即使设置为相同的 CRS

Question

我有一个关于坐标的小问题，我想看看哪些观测值属于哪些县。我正在使用 NHGIS 县边界：

> counties = st_read("path/US_county_2012.shp", stringsAsFactors = FALSE)

当我使用“st_crs()”提取 CRS 时，我得到的输出似乎与解释 R 的 sf 包的论坛和教程中其他人的输出不匹配：

> st_crs(counties)
Coordinate Reference System:
  User input: USA_Contiguous_Albers_Equal_Area_Conic 
  wkt:
PROJCRS["USA_Contiguous_Albers_Equal_Area_Conic",
    BASEGEOGCRS["NAD83",
        DATUM["North American Datum 1983",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["USA_Contiguous_Albers_Equal_Area_Conic",
        METHOD["Albers Equal Area",
            ID["EPSG",9822]],
        PARAMETER["Latitude of false origin",37.5,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",-96,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",29.5,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",45.5,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["USA - CONUS - onshore"],
        BBOX[24.41,-124.79,49.38,-66.91]],
    ID["ESRI",102003]]

所以我使用相同的 CRS 将我的 tibble 转换为 sf 对象：

> head(dupes_tibble)
index       lon      lat
1  7911 -84.60410 33.44512
2  5211 -85.57854 42.88454
3  7075 -85.53756 42.86731
4  6600 -85.53756 42.86731
5  2042 -95.71289 37.09024
6  2553 -77.44137 38.30777
> dupes_sf = st_as_sf(dupes_tibble, coords = c("lon", "lat"), crs = st_crs(counties))

然后我看看哪些观测值属于哪些县：

dupes_county = st_join(dupes_sf, counties, join = st_within)

我没有收到错误，但是当我将 ggplot 中的地图与坐标应该是什么进行比较时，很明显我的 tibble 点有偏差。我猜测由于某种原因将我的点设置为相同的 CRS 无法正确转换坐标，可能是因为 st_crs() 的奇怪输出。可能有人知道我可能做错了什么吗？

Answer 1

当您从您的点创建 sf 对象时，您首先需要指定 lat/lon 坐标在 WGS84 coordinate system (EPSG 4326) 中。然后从那里开始，下一步是将点转换为与多边形相同的 CRS。

这里是一个例子，我使用tigris package下载县界，但它应该和你的NHGIS shp文件一样。

library(tidyverse)
library(sf)
library(tigris)

# get county boundaries
counties <- counties(cb = TRUE, class = "sf") %>% 
  filter(!as.numeric(STATEFP) %in% c(2, 15, 60, 66, 69, 72, 78)) # lower 48 only

# create points data 
dupes_tibble <- tribble(
  ~index,   ~lon,  ~lat,
  7911, -84.60410, 33.44512,
  5211, -85.57854, 42.88454,
  7075, -85.53756, 42.86731,
  6600, -85.53756, 42.86731,
  2042, -95.71289, 37.09024,
  2553, -77.44137, 38.30777
  ) 

# convert to sf and transform crs
dupes_sf <- dupes_tibble %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>% 
  st_transform(st_crs(counties))

# plot
ggplot() +
  geom_sf(data = counties) +
  geom_sf(data = dupes_sf)

Point-in-Polygon 但坐标系关闭，即使设置为相同的 CRS

Point-in-Polygon but coordinate system off even when set to same CRS

gis

r

latitude-longitude

point-in-polygon

sf