R：当点重叠/在距离内时追加数据；将缓冲区矩形添加到 set1，将半径添加到 set2

Posted 2023-03-23

【中文标题】R：当点重叠/在距离内时追加数据；将缓冲区矩形添加到 set1，将半径添加到 set2

【英文标题】：R: Append data when points overlap/within distance; add buffer rectangle to set1, add radius to set2

【发布时间】：2019-10-14 10:57:37

【问题描述】：

我有一个涡流数据库 (netcdf)，其中包含在海洋中移动的涡流，以及许多鱼的轨迹。当鱼在涡流范围内（（变化的）阈值距离内的位置和日期相同）时，我试图将涡流信息附加到鱼道。数据结构、复杂性和可重现的示例如下。我怀疑这不是超级难，我只是不断地打结，试图找出一个巧妙的解决方案。提前感谢您的任何想法！

数据结构/可重现示例：3 天的鱼和涡流轨迹

library(lubridate)
fish <- data.frame(lat = c(42.1, 42.6, 43.2), lon = c(-10, -10.1, -10.2), date = ymd(c("1990-01-01", "1990-01-02", "1990-01-02")))
eddy <- data.frame(lat = c(44, 42.3, 40), lon = c(-15, -10.1, -6), date = ymd(c("1990-01-01", "1990-01-02", "1990-01-02")), track = c(1,1,1), radius = c(81,82,83), vals = c(11,12,13))

这里有点棘手：

变半径环形涡流缓冲器

涡流的大小（从中计算重叠或距离）由半径列给出，所以：

sf_eddy <- st_as_sf(eddy)
st_buffer(sf_eddy, dist = radius)

对于 st_buffer “所有操作都基于每个功能”，所以我不确定这是否会将 sf_eddy 视为单个功能（可能是（多）线串）或多个单独的点。我想我可以玩弄它，直到我让它工作，这样每个涡流日（n = 151937）都是一个 sf 缓冲圈，我大概可以将其 vals 附加到，如果它们还没有包括在内。

设置距离 latlong 矩形鱼缓冲区

鱼的位置有一个固定的误差距离，+-1.9°lat，+-0.77°lon。所以从概念上讲，我想在每天的鱼位置周围创建一个又高又细的矩形多边形。 st_buffer 似乎不允许这样做，但 there may be other options 可能会根据 square 示例 here 在每个质心周围创建一个多边形。

如果鱼和任何涡流在空间和时间上重叠，则将涡流 val 附加到鱼

当 eddy$date==fish$date?虽然涡流计算是一次性的，但我可以保存输出。我怀疑鱼箱计算也不会特别费力。无论如何：

fish$vals <- rep(NA, nrow(fish))    
sf_fish <- sf_as_sf(fish)
for (i in seq(along = sf_fish))
if(st_intersects(sf_fish[i,],sf_eddy)) sf_fish[i,"vals"] <- sf_eddy$vals

显然，这可能是一种糟糕的做法。

但是，这是否是解决这个问题的一种模糊的合理方法，还是由于我对 sf、sp、rgeos 等的业余知识，我错过了更优雅的解决方案？我有 166 个鱼文件，平均 286 天的位置 = 47476 总查找。创建鱼/涡流日期匹配对的索引（同一天可能有多个涡流）然后仅在这些上测试交叉点会好/快吗？由于我需要将eddy$vals（实际上是2或3列）附加到fish，所以使用空间连接会更好，也许st_join？

编辑：评论和输出地图怪事： 600,000 条鱼/涡流 ggplot 总输出：

看起来我需要修复投影，因为涡流没有在 0 以西打印。它们现在肯定在同一个 CRS 中；也许原件是 0-360 而不是 -180:180。此外，虽然鱼箱打印得很好，但涡流太小：意味着 r=80km=~1deg@45N，而且它们的打印肯定比这小得多。我也会考虑的。另外：如果同一天有2个涡流，for()会有问题吗？

Edit2：我更正的原件是 0:360。涡流半径仍然太小： 根据代码讨论，带有 latlons 的原始 eddy 文件被转换为 st_as_sf 到 crs4326 (wgs84) 然后 st_transform(sf_eddy,6931) 这是北大西洋，如上图所示。结果文件中的几何列是 POLYGON，XY 尺寸，值是投影坐标而不是纬度，例如sf_eddy_buffered$geometry[1]bbox:xmin -4894790 ymin -2418555 xmax -4894648 ymax -2418413。 speed_radius 是 71，所以 xmin 和 xmax 之间的差异是 142，即 71*2，这是正确的。半径单位不等于投影距离单位吗？斯图尔特的例子是 81000。涡流半径单位：公里。投影/几何单位：https://epsg.io/6931 表示单位确实是米。将半径乘以 1000 并重试： 所以这一切都很好。刚刚尝试了更新的st_join 代码，认为它可以使用一些工作 - 可能因为 3639209 涡流行与 461 条鱼行而需要永远。它还在同一天为同一条鱼创建了多个涡流的多行。考虑到鱼周围的错误缓冲区框，可能是相对较大的重叠机会的函数。 我怀疑我可以通过标记重复项并删除最远的重复项来解决此问题。当前的执行方式（通过 v 和 rbind）比以前的方式（值的直接单元分配）IMO 更好，因为否则多个值将被发送到单个单元，要么崩溃，要么静默覆盖。

为了加快代码速度，我想我可以将 eddy 文件子集设置为仅当前鱼的日期范围内的重叠部分，并在连接循环中使用该子集。无论如何，连接循环已经按日期对鱼和涡流文件进行了子集化，所以这可能是多余的，或者可能减少需要处理的涡流数据的大小每次都会减少 CPU/RAM 的工作量？

Edit3：这大大加快了速度。为thisDate 添加一个简单的循环计数器显示它在 461 行中的 275 行完成。并制作了一个481长度的文件。可能有很多天鱼不在漩涡中（thisDate 循环导致没有v 和没有rbind），并且有很多天鱼在两个以上的漩涡中（@987654354 @ 循环导致多行 v 和 rbind)。仍然感觉它应该处理 461 循环，但是......

Edit4：一个小调整：

v <- st_join(thisFishRow, thisEddyRow, left = F, largest = TRUE)

largest 表示它只会加入具有最大重叠区域的涡流，导致最大 v 行长为 1，因此每天最多 1 个值。然而，路径的细微差别意味着鱼似乎可以在不同的持续漩涡之间来回穿梭。我的怀疑（基于类似的研究）是它们更有可能留在强大的反气旋涡流系统中。涡流控制的半径（以及面积）可能重叠超过位置。对于 1.5 个月的部分，我发现重叠度最高的涡流轨道编号在 3 个涡流之间反弹。涡道编号 217008 是最强大的反气旋。这是一个（废话）情节： 如您所见，217008 正好位于鱼错误框的中心，而更大、更分散、功率更低的涡流 217343 和 39625 位于边缘。然而，它们较大的尺寸经常会看到它们被撞到顶部，因为它们有更多的区域可以提供，并且不计算质心接近度。 所以：如果 fishbox 在同一天与 eddy 重叠，则将 eddy 包括在候选名单中 （thisFishRow & thisEddyRow 保持不变）。然后：根据与质心的最短距离（而不是st_join）从候选名单中选择涡流。待续！

编辑5：

fishNearEddies <- NULL
for (thisDate in sf_df_nona$date) 
 thisFishRow <- sf_df_nona[sf_df_nona$date == thisDate, ]
 thisEddyRow <- sebDateSub[sebDateSub$date == thisDate, ]
 overlapToday <- st_join(thisFishRow, thisEddyRow, left = F)
 if (nrow(overlapToday) > 0) 
  overlapEddies <- thisEddyRow[thisEddyRow$track %in% overlapToday$track,] 
  st_distance(x = thisFishRow, y = overlapEddies, by_element = TRUE) #0 0 due to overlap min dist = 0
  #close if
 #close for

st_distance 不起作用，因为空间特征（缓冲区和缓冲区圈）重叠，所以最小距离 = 0。我想我还需要质心来测试？

Edit6：使用工作代码进行最终编辑，并为 Stewart 的所有帮助提供答案。再次感谢先生。

# overlap join loop####
fishNearEddies <- NULL
EdCentroidDist <- NULL
counter <- 1
ofhowmany <- length(sf_df_nona$date)
for (thisDate in sf_df_nona$date) 
  thisFishRow <- sf_df_nona[sf_df_nona$date == thisDate, ] #will be 1 row per day
  thisEddyRow <- sebDateSub[sebDateSub$date == thisDate, ] #all eddies that day, could be multi rows
  overlapToday <- st_join(thisFishRow, thisEddyRow, left = F) #will join only if they spatial intersect and (already) time match
  if (nrow(overlapToday) > 0) 
    # now need to join based on closest centroid and NOT on highest overlap
    overlapEddies <- thisEddyRow[thisEddyRow$track %in% overlapToday$track,] #subset TER by overlap tracks
    fishcentroid <- st_centroid(thisFishRow)
    eddycentroid <- st_centroid(overlapEddies)
    fishEdDists <- st_distance(x = fishcentroid, y = eddycentroid, by_element = TRUE) #result vector corresponding to overlapEddies rownumbers
    fishNearEddies <- rbind(fishNearEddies, overlapEddies[which.min(fishEdDists),]) #row index for overlapEddies, no index but has date
    EdCentroidDist <- rbind(EdCentroidDist, as.numeric(min(fishEdDists)))
   #close if
  print(paste0(counter, " of ", ofhowmany, " fish days"))
  counter <- counter + 1
 #close for

if (!is.null(fishNearEddies))   #if there are overlaps, do processing. If not it'll fail

  fishNearEddies %<>% as.data.frame() %>% # convert to nonspatial so I can remove buffer column
    select(-geometry) %>% # remove geometry column. Attributes still remain. Whatever?
    cbind(EdCentroidDist) #add to FNE as column

【问题讨论】：

如何根据要素的质心计算距离（使用st_centroid），然后选择距离最小的要素？

【参考方案1】：

这里有一些东西可以帮助你：

library(lubridate)
library(sf)
library(ggplot2)

# sample data
fish <- data.frame(
    lat = c(41.1, 43.6, 44.2),
    lon = c(-7, -11, -15),
    date = ymd(c("1990-01-01", "1990-01-02", "1990-01-03"))
)

# Add a blank column for the eddy values we want
fish$vals <- rep(NA, nrow(fish))

eddy <- data.frame(
    lat = c(44, 42.3, 40),
    lon = c(-6, -10.1, -15),
    date = ymd(c("1990-01-01", "1990-01-02", "1990-01-03")),
    track = c(1,1,1),
    radius = c(81000,82000,83000),
    vals = c(11,12,13)
)

# Convert eddy to simple features, using WGS84 CRS
sf_eddy <- st_as_sf(eddy,  coords = c("lon", "lat"), crs=4326)

# Transform from geographical to projected so we can buffer correctly. You might need to pick a different CRS
sf_eddy <- st_transform(sf_eddy, 3035) # ETRS89 / LAEA Europe

# Buffer the eddy points based on their radii
sf_eddy_buffered <- st_buffer(sf_eddy, dist = sf_eddy$radius)

# Add error to fish position. There's probably a better way to do this.
fishErrLat <- 1.9
fishErrLon <- 0.77
fish$buffer <- paste('POLYGON((',
                    fish$lon - fishErrLon, ' ', fish$lat + fishErrLat, ',',
                    fish$lon + fishErrLon, ' ', fish$lat + fishErrLat, ',',
                    fish$lon + fishErrLon, ' ', fish$lat - fishErrLat, ',',
                    fish$lon - fishErrLon, ' ', fish$lat - fishErrLat, ',',

                    fish$lon - fishErrLon, ' ', fish$lat + fishErrLat,
                    '))',
                    sep=''
                )

# Convert fish to simple features, using WGS84 CRS
sf_fish <- st_as_sf(fish, wkt='buffer', crs=4326)

# Transform from geographical to projected
sf_fish <- st_transform(sf_fish, 3035)


#Plot what we've got so far    
g <- ggplot() + 
        geom_sf(data=sf_eddy_buffered, aes(fill=date)) +
        geom_sf(data=sf_fish, aes(fill=date))

print(g)

# Check for overlap
fishNearEddies <- NULL

for (thisDate in unique(sf_fish$date)) 

    thisFishRow <- sf_fish[sf_fish$date==thisDate, ]
    thisEddyRow <- sf_eddy_buffered[sf_eddy_buffered$date==thisDate, ]

    v <- st_join(thisFishRow, thisEddyRow, left=F)
    if (nrow(v) > 0) 
        fishNearEddies <- rbind(fishNearEddies, v)
    

并检查结果：

> fishNearEddies
Simple feature collection with 1 feature and 7 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 2588691 ymin: 2402830 xmax: 2703342 ymax: 2624501
epsg (SRID):    3035
proj4string:    +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
   lat lon     date.x     date.y track radius vals                         buffer
2 43.6 -11 1990-01-02 1990-01-02     1  82000   12 POLYGON ((2644792 2624501, ...

这只会给你那些与涡流相交的鱼。

【讨论】：

是的，for 循环需要进行编辑以处理多个匹配项。如果鱼与两个不同的vals 涡流相匹配，你如何决定选择哪个val？

我刚刚注意到您的误差为 ±1.7˚，因此我更新了代码以删除 1.7/2 中的分母。您的数据来自哪种跟踪技术？

我认为你使用sf_join 的想法比我原来使用st_intersects 的for 循环要好。我已经更新了我的答案。

你现在已经超出了我的知识水平。祝你好运！

谢谢斯图尔特；现在上传的答案，可能不是最优雅的，但似乎可以完成这项工作！干杯