使用 over(), R 连接 SpatialPointsDataFrame 和 SpatialLinesDataFrame

Join SpatialPointsDataFrame and SpatialLinesDataFrame using over(), R

我遇到了一些看起来应该很简单的事情。抱歉,我不熟悉在 R 中使用空间数据。

我正在尝试将城市数据映射到世界海岸线地图上。我从自然地球数据集 (https://www.naturalearthdata.com/downloads/) 1:110m 数据中提取了海岸线并生成了空间线数据框:

coast_rough_sldf
class       : SpatialLinesDataFrame 
features    : 134 
extent      : -180, 180, -85.60904, 83.64513  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 3
names       : scalerank, featurecla, min_zoom 
min values  :         0,  Coastline,      0.0 
max values  :         1,    Country,      1.5 

我还有一个城市数据集,样本如下:

city_coast <- data.frame(Latitude = c(-34.60842, -34.47083, -34.55848, -34.76200, -34.79658, -34.66850), 
              Longitude = c(-58.37316, -58.52861, -58.73540, -58.21130, -58.27601, -58.72825), 
              Name1 = c("Buenos Aires", "San Isidro", "San Miguel", "Berazategui", "Florencio Varela", "Merlo"), 
              distance = c(7970.091,  5313.518, 26156.700, 11670.274, 18409.738, 33880.259))
city_coast

Latitude Longitude            Name1  distance
1 -34.60842 -58.37316     Buenos Aires  7970.091
2 -34.47083 -58.52861       San Isidro  5313.518
3 -34.55848 -58.73540       San Miguel 26156.700
4 -34.76200 -58.21130      Berazategui 11670.274
5 -34.79658 -58.27601 Florencio Varela 18409.738
6 -34.66850 -58.72825            Merlo 33880.259

然后我成功创建了空间点数据框:

city_spdf <- SpatialPointsDataFrame(coords = select(city_coast, c("Longitude", "Latitude")),
                                    proj4string = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84"),
                                    data = select(city_coast, c("Name1", "distance")))

city_spdf

class       : SpatialPointsDataFrame 
features    : 6 
extent      : -58.7354, -58.2113, -34.79658, -34.47083  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 2
names       :       Name1,  distance 
min values  : Berazategui,  5313.518 
max values  :  San Miguel, 33880.259 

现在我想加入 city_spdf 和 coast_sldf,这样我就可以使用 tmap 绘制它们。查看教程似乎我应该使用 over():

city_coast_shp <- over(coast_rough_sldf, city_spdf)

city_coast_shp

Name1 distance
1  <NA>       NA

这显然是错误的。切换对象的顺序会改变一些事情,但仍然无法满足我的需求。

谁能告诉我这个 over 函数有什么不对的地方?我见过的每个例子都只是让人们加入两个空间对象。如果我遗漏了一些非常简单的东西,我们深表歉意。

就像@elmuertefurioso 在评论中指出的那样,我认为这没有按照您预期的方式工作的一个原因是几何类型的混淆。

由于 coastline 数据是线,而不是像 tmap 中的 data(World) 这样的多边形,因此您在使用 [=17= 进行计算和比较时会受到一些限制],也就是积分。

sf方式读取数据:

library(sf)

# downloaded from https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/physical/ne_110m_coastline.zip
coastline <- read_sf("~/Downloads/ne_110m_coastline/ne_110m_coastline.shp")

cities <- data.frame(
  Latitude = c(-34.60842, -34.47083, -34.55848, -34.76200, -34.79658, -34.66850), 
  Longitude = c(-58.37316, -58.52861, -58.73540, -58.21130, -58.27601, -58.72825), 
  Name1 = c("Buenos Aires", "San Isidro", "San Miguel", "Berazategui", "Florencio Varela", "Merlo"), 
  distance = c(7970.091,  5313.518, 26156.700, 11670.274, 18409.738, 33880.259)
  )

为了在 sf 对象之间进行任何比较,它们必须具有相同的坐标参考系统。因此,当我们读入 cities 时,我们会将 CRS 设置为 coastline.

的 CRS
cities <- st_as_sf(
  cities,
  coords = c("Longitude", "Latitude"), # must be x, y order
  crs = st_crs(coastline) # must be equivilant between objects
  )

现在您可以使用 st_{comparison}() 系列函数进行比较。

函数 over() 及其 sf 对应函数 st_intersects() 可以处理一组点和多边形,但我们这里没有。我们可以对点和线使用 st_nearest_feature() 等距离函数,为每个城市获取距 coastline 最近的几何图形。

st_nearest_feature(cities, coastline)

它 return 是 coastlines 中最近的几何图形的行索引,这里的所有城市恰好是相同的,因为它们都在阿根廷。顺序在函数中很重要,因为它定义了被问到的问题如果我们将其翻转为 st_nearest_feature(coastline, cities),它将 return 是 coastline 中每个几何图形最近的城市,因此 return将有 134 个元素。

也就是说,您实际上不必进行任何连接或比较即可将您的点绘制在同一个 tmap 上。

library(tmap)

tmap_mode("view")

tm_shape(coastline) +
  tm_lines() +
tm_shape(cities) +
  tm_bubbles("distance")

我不是 tmap 用户,但我只是放大了这个屏幕截图以显示它的工作原理。