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

问题描述

我有一个涡流在海洋中移动的数据库（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看起来不允许这样做，但可能还有其他选项，可能是根据此处square的示例在每个质心周围创建一个多边形。

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

当 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 的原始涡流文件被转换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 更好，因为否则多个值将被发送到单个单元格，要么崩溃，要么默默地覆盖。

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

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

Edit4：一个小调整：

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

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

标签： rbufferintersectionspsf

解决方案

Here's something that should get you going:

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)
    }

}

And check the results:

> 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, ...

This will give you only those fish that intersect with an eddy.

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

问题描述

解决方案

推荐阅读