将点和多边形数据与 R 中的 sf 相交时结果不正确

问题描述 投票：0回答：1

当我使用 R 中的 

sf
包将一些点数据与多边形相交时，我得到了奇怪的结果，虽然我认为我已经找到了解决方案，但我对这个问题感到非常困惑，我真的很感激它是否有人可以帮助我了解正在发生的事情以及我是否遗漏了任何潜在的问题。

我有一个只有坐标的数据集；我需要知道每个点位于哪个 2000 个区域中，因此我想将其与 2000 个区域 shapefile 相交，以便分配正确的区域 ID。我的真实数据是保密的，因此举个例子，我目前正在使用 DC 地铁地图，可在此处获取。我的真实数据只是坐标，而不是形状文件，因此最初缺乏投影。我将从地铁地图中删除投影以模仿我的实际数据（并复制数据，以便我可以展示我的解决方案）：

library(sf)

# open metro shapefile using file paths not shown here
metro <- read_sf(dsn = metro.fp, layer = metro.filename)

dim(metro)
[1] 98 15

# remove metro projection to mirror initial state of real data
st_crs(metro) <- ""

# duplicate this data for later example of my solution
metro2 <- metro

我使用的 shapefile 包含从 NHGIS.org 下载的 2000 个美国人口普查区（2010 年 TIGER/Line 几何）的多边形。 （如果有一种方法可以更轻松地分享这个示例，请告诉我；我可以使用 

dput()
，但即使与 
head()
结合使用，输出也是巨大的。）它使用的是 ESRI 投影，我认为是与当前的问题相关。当我查看它的投影时，我不会显示完整的输出，因为它很长，但是当您使用 
st_crs()
时，它会显示该地图正在使用 
USA_Contiguous_Albers_Equal_Area_Conic
，底部有 
ID["ESRI",102003]
。

当我将点数据集转换为与多边形相同的投影时 - 使用 

st_crs()
而不是 
st_transform()
因为我的点数据缺乏开始的投影 - 然后将两者相交，我最终遇到了这个奇怪的问题我的所有积分据说都在堪萨斯州（州 fips 为 20），尽管这些都是整个华盛顿特区、马里兰州和弗吉尼亚州的地铁站：

# open tract map using file paths not shown here
tract00 <- read_sf(dsn = fp, layer = file.name)

# give metro same projection as tract map
st_crs(metro) <- st_crs(tract00)

# intersect the two
metro.tract <- st_intersection(x = metro,
                               y = tract00)

# now all states are kansas (state fips of 20):
table(metro.tract$STATEFP00)

20 
98

如果我给我的初始地铁数据集一个随机投影，然后使用 

st_transform
将其放入与区域 shapefile 相同的投影中，那么我正确地知道车站存在于 DC、马里兰州和弗吉尼亚州（DC 为 11， MD 为 24，VA 为 51）：

# now try again, first giving metro2 a random projection so I can use st_transform
st_crs(metro2) <- "WGS84" 

# give metro2 same projection as tract map
metro2 <- st_transform(metro2, crs = st_crs(tract00))

# intersect the two
metro.tract2 <- st_intersection(x = metro2,
                                y = tract00)

# now the states are correct:
table(metro.tract2$STATEFP00)

11 24 51 
40 26 32

有人知道发生了什么事吗？虽然我的解决方案生成的结果看起来正确，但因为我不明白问题到底是什么，我担心我可能会错过真实数据的其他潜在问题。

我最好的猜测是，这与区域形状文件的 ESRI 投影有关，因为 st_transform 的 sf 文档说：“ESRI Shapefile 格式需要额外小心，因为 WKT1 不会明确存储轴顺序。 ”但这基本上只是我的猜测，我不确定为什么我的解决方案在这种情况下有效。

基本上，如果有人可以告诉我发生了什么事或者我是否遗漏了任何其他潜在问题，我将不胜感激。

r gis shapefile esri sf-package
1个回答
0
投票

我将依次阐述您的观点以澄清这个问题。首先，加载所需的包并创建一些代表性数据以供使用：

library(tigris) # For US census tracts
library(sf)
library(dplyr)
library(ggplot2)
library(usmap)

# Read tract data for DC. Project to ESRI:102003 as default CRS is NAD83:EPSG4269 
tract00 <- tracts(state = "dc", year = 2000) %>%
  st_transform("ESRI:102003")

tract00$geometry
# Geometry set for 188 features 
# Geometry type: MULTIPOLYGON
# Dimension:     XY
# Bounding box:  xmin: 1610777 ymin: 306994.4 xmax: 1629412 ymax: 329361
# Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic
# First 5 geometries:
# MULTIPOLYGON (((1617395 319019.9, 1617392 31902...
# MULTIPOLYGON (((1617767 319815.2, 1617772 31981...
# MULTIPOLYGON (((1616807 318539.2, 1616794 31854...
# MULTIPOLYGON (((1617619 318298.1, 1617601 31829...
# MULTIPOLYGON (((1618010 318832.6, 1618021 31883...

# DC metro data shp previously unzipped to working directory, downloaded from 
# https://opendata.dc.gov/datasets/DCGIS::metro-stations-regional/explore
metro <- st_read("Metro_Stations_Regional.shp") 

metro$geometry
# Geometry set for 98 features 
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: -77.49154 ymin: 38.76653 xmax: -76.84455 ymax: 39.11994
# Geodetic CRS:  WGS 84
# First 5 geometries:
# POINT (-76.91147 38.82645)
# POINT (-77.05367 38.81415)
# POINT (-77.06081 38.80659)
# POINT (-77.07088 38.80043)
# POINT (-77.07521 38.79392)

注意metro 和tract00 之间的值差异。

当我将点数据集转换为与多边形相同的投影时......我最终遇到了这个奇怪的问题，其中我的所有点都应该在堪萨斯州......

您链接到的 DC Metro 数据集的坐标单位是十进制度，而区域数据的坐标单位是米。通过运行：

st_crs(metro) <- st_crs(tract00)
# Warning message:
#   st_crs<- : replacing crs does not reproject data; use st_transform for that

你是说那些十进制度值实际上是米。这正是您这样做时会收到警告的原因。结果是您的地铁数据在纵向上距离应有的位置 > 1.6e6，在纬度上距离应有 > 3.06e5，例如多萝西和托托回到堪萨斯州;)

如果我给我的初始地铁数据集一个随机投影，然后使用

st_transform
将其放在与区域形状文件相同的投影中，那么我正确地知道车站存在于DC...

事实证明，在这个假设中，您通过将 WGS84:EPSG4326 分配给 Metro2 数据而变得“幸运”。这是因为我们知道原始 Metro 文件的 CRS 是 WGS84。所以当你跑步时：

st_crs(metro2) <- "WGS84" 

metro2 <- st_transform(metro2, crs = st_crs(tract00))


st_transform()
函数“我看到metro的单位是十进制，并且我根据CRS知道它的位置，所以我将这些十进制度坐标值转换为米，并使用tract00的CRS来定位它们。”

但是，正如我所说，你很“幸运”。对于您的实际数据集，您尚未指明哪个 CRS 构成坐标值的基础。这是一个问题，因为地理坐标系 (GCS) 值看起来相似，但可能差异很大。考虑以下几点：

# Reload metro dataset
metro <- st_read("Metro_Stations_Regional.shp")

# Convert metro to df with 'unknown' coordinates to replicate your actual data
metro2 <- data.frame(id = 1:nrow(metro),
                     st_coordinates(metro))

head(metro2)
#   id         X        Y
# 1  1 -76.91147 38.82645
# 2  2 -77.05367 38.81415
# 3  3 -77.06081 38.80659
# 4  4 -77.07088 38.80043
# 5  5 -77.07521 38.79392
# 6  6 -76.99537 38.86297

# Coordinates correctly assumed to be WGS84
p_4236 <- st_as_sf(metro2, coords = c("X", "Y"), crs = 4326)

# Coordinates incorrectly assumed to be NAD27/EPSG:4267
p_4267 <- st_as_sf(metro2, coords = c("X", "Y"), crs = 4326) %>% # lon/lat
  st_transform(4267) %>%
  st_set_crs(4326)
# Warning message:
#   st_crs<- : replacing crs does not reproject data; use st_transform for that

p_4267$geometry
# Geometry set for 98 features 
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: -77.49185 ymin: 38.76649 xmax: -76.84488 ymax: 39.11991
# Geodetic CRS:  WGS 84
# First 5 geometries:
# POINT (-76.91179 38.82642)
# POINT (-77.05399 38.81412)
# POINT (-77.06114 38.80656)
# POINT (-77.0712 38.8004)
# POINT (-77.07553 38.79389)

# Return distance between original and incorrectly projected location
head(st_distance(metro, p_4267 , by_element = TRUE))
# Units: [m]
# [1] 28.61364 28.23506 28.21754 28.19191 28.18171 28.38108

尽管 Metro 与 p_4267 具有相同的 CRS，并且坐标值看起来相似，但实际上它们相距约 28m。总之，在预测之前，您需要知道哪个 CRS 构成了您的实际数据的基础。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.