r语言 - 如何在第二个 df 中查找最接近坐标的值/匹配值

Question

我在海上有一系列的地理位置，我正在尝试获取地质沉积物类型信息。我正在使用英国国家地质沉积物数据库（df1）的导出，这是一个坐标和沉积物信息的大型数据集。目前，我一直在对BGS导出文件（df1）中的坐标进行四舍五入，并平均/重新计算这些坐标方块的沉积物类型，然后我在（df2）中对坐标进行四舍五入，并将其与这些方块匹配以获得沉积物分类。

BGS 导出如下所示 （df1）;

    NUM     X       Y           GRAV    SAND    MUD
1   228     1.93656 52.31307    1.07    98.83   0.10
2   142     1.84667 52.45333    0.00    52.60   47.40
3   182     1.91950 52.17750    9.48    90.38   0.14
4   124     1.88333 52.70833    0.00    98.80   1.20
5   2807    1.91050 51.45000    2.05    97.91   0.05
6   2787    1.74683 51.99382    41.32   52.08   6.60
7   2776    1.66117 51.63550    9.83    87.36   2.81
8   2763    1.82467 51.71767    43.92   47.25   8.83
9   2753    1.76867 51.96349    57.66   39.18   3.15
10  68      2.86967 52.96333    0.30    98.90   0.80
11  2912    1.70083 51.77783    26.90   64.87   8.22
12  2914    1.59750 51.88882    32.00   65.02   2.97
13  2886    1.98833 51.34267    1.05    98.91   0.04
14  2891    1.87817 51.31549    68.57   31.34   0.08
15  2898    1.37433 51.41249    35.93   61.48   2.59
16  45      2.06667 51.82500    9.70    88.10   2.20
17  2904    1.63617 51.45999    16.28   66.67   17.05

我在海上的位置看起来像这样（df2）;

haul    DecStartLat DecStartLong
1993H_2 55.23983    -5.512830
2794H_1 55.26670    -5.516700
1993H_1 55.27183    -5.521330
0709A_71    55.26569    -5.519730
0396H_2 55.44120    -5.917800
0299H_2 55.44015    -5.917310
0514A_26    55.46897    -5.912167
0411A_64    55.47289    -5.911820
0410A_65    55.46869    -5.911930
0514A_24    55.63585    -5.783500
0295H_4 55.57250    -5.754300
0410A_62    55.63656    -6.041870
0413A_53    55.73280    -6.020600
0396H_13    55.66470    -6.002300
2794H_8 55.83330    -5.883300
0612A_15    55.84025    -5.912130
0410A_74    55.84311    -5.910180
0299H_16    55.90568    -5.732490
0200H_18    55.88600    -5.742900
0612A_18    55.90450    -5.835880

这是我的剧本...

get.Sed.type <- function(x,y) {
  x$Y2 <- round(x$Y, digits=1)
  x$X2 <- round(x$X, digits=1)
  x$BGSQ <- paste(x$Y2,x$X2,sep="_")
  x$RATIO <- x$SAND/x$MUD
  x <- aggregate(cbind(GRAV,RATIO)~BGSQ,data=x,FUN=mean)
  FOLK <- (x$GRAV)
  FOLK[(FOLK)<1] <- 0
  FOLK[(FOLK)>=1&(FOLK)<5] <- 1
  FOLK[(FOLK)>=5&(FOLK)<30] <- 5
  FOLK[(FOLK)>=30&(FOLK)<80] <- 30
  FOLK[(FOLK)>=80] <- 80
  R_CLASS <- (x$RATIO)
  R_CLASS[(R_CLASS)<1/9] <- 0
  R_CLASS[(R_CLASS)>=1/9&(R_CLASS)<1] <- 0.1
  R_CLASS[(R_CLASS)>=1&(R_CLASS)<9] <- 1
  R_CLASS[(R_CLASS)>=9] <- 9
  x$FOLK_CLASS <- NULL
  x$FOLK_CLASS[(R_CLASS)==0&(FOLK)==0] <- "M"
  x$FOLK_CLASS[(R_CLASS)%in%c(0,0.1)&(FOLK)==5] <- "gM"
  x$FOLK_CLASS[(R_CLASS)==0.1&(FOLK)==0] <- "sM"
  x$FOLK_CLASS[(R_CLASS)==0&(FOLK)==1] <- "(g)M"
  x$FOLK_CLASS[(R_CLASS)==0.1&(FOLK)==1] <- "(g)sM"
  x$FOLK_CLASS[(R_CLASS)==9&(FOLK)==0] <- "S"
  x$FOLK_CLASS[(R_CLASS)==1&(FOLK)==0] <- "mS"
  x$FOLK_CLASS[(R_CLASS)==9&(FOLK)==1] <- "(g)S"
  x$FOLK_CLASS[(R_CLASS)==1&(FOLK)==1] <- "(g)sM"
  x$FOLK_CLASS[(R_CLASS)==1&(FOLK)==5] <- "gmS"
  x$FOLK_CLASS[(R_CLASS)==9&(FOLK)==5] <- "gS"
  x$FOLK_CLASS[(FOLK)==80] <- "G"
  x$FOLK_CLASS[(R_CLASS)%in%c(0,0.1)&(FOLK)==30] <- "mG"
  x$FOLK_CLASS[(R_CLASS)==1&(FOLK)==30] <- "msG"
  x$FOLK_CLASS[(R_CLASS)==9&(FOLK)==30] <- "sG"
  y$Lat <- round(y$DecStartLat, digits=1)
  y$Long <- round(y$DecStartLong, digits=1)
  y$LATLONG100_sq <- paste(y$Lat,y$Long,sep="_")
  y <- merge(y, x[,c(1,4)],all.x=TRUE,by.x="LATLONG100_sq",by.y="BGSQ")
  #Delete unwanted columns
  y <- y[, !(colnames(y) %in% c("Lat","Long","LATLONG100_sq"))]
  #Name column something logical
  colnames(y)[colnames(y) == 'FOLK_CLASS'] <- 'BGS_class'
  return(y)
}

但是，我在 db2 中有十几个位置，BGS 导出 （db1） 中没有相应的值，我想知道如何要求它为相应正方形周围的正方形做另一个平均值（即四舍五入到更大的数字并重复该过程）或者要求它在 BGS 导出文件中找到最接近的坐标并采用现有的价值。

Answer 1

对于问题中提到的第二个选项，我建议将问题框定如下：

假设您有一组来自 db1 的 m 坐标和来自 db2 的 n 个坐标，m <=n，并且当前这些集合的交集为空。

您希望将 db1 中的每个点与 db2 中的点进行匹配，以便将匹配的"误差"（例如距离总和）最小化。

解决这个问题的一种简单贪婪的方法可能是生成一个 m x n 矩阵，其中包含每对坐标之间的距离，并按顺序为每个点选择最接近的匹配项。当然，如果有很多点需要匹配，或者如果您追求最优解决方案，您可能需要考虑更复杂的匹配算法（例如匈牙利算法）。

法典：

  #generate some data (this data will generate sub-optimal matching with greedy matching)
  db1 <- data.frame(id=c("a1","a2","a3","a4"), x=c(1,5,10,20), y=c(1,5,10,20))
  db2 <- data.frame(id=c("b1","b2","b3","b4"),x=c(1.1,2.1,8.1,14.1), y=c(1.1,1.1,8.1,14.1))
  #create cartesian product
  product <- merge(db1, db2, by=NULL)
  #calculate auclidean distances for each possible matching
  product$d <- sqrt((product$x.x - product$x.y)^2 + (product$y.x - product$y.y)^2)
  #(naively & greedily) find the best match for each point
  sorted <- product[ order(product[,"d"]), ]
  found <- vector()
  res <- vector() #this vector will hold the result
  for (i in 1:nrow(db1)) {
    for (j in 1:nrow(sorted)) {
      db2_val <- as.character(sorted[j,"id.y"])
      if (sorted[j,"id.x"] == db1[i, "id"] && length(grep(db2_val, found)) == 0) {    
        #print(paste("matching ", db1[i, "id"], " with ", db2_val))
        res[i] <- db2_val      
        found <- c(found, db2_val)
        break
      }
    }
  }

请注意，我相信通过使用循环以外的方法可以改进代码并使其更加优雅。

Answer 2

希望我没有误解，但就我从标题中得到的，您需要根据最小距离进行匹配。如果允许这个距离是欧几里得距离，那么可以使用快速 RANN 包，如果不是，那么就需要计算大圆距离。

提供的一些数据

BGS_df <- 
  read.table(text = 
               "    NUM     X       Y           GRAV    SAND    MUD
                1   228     1.93656 52.31307    1.07    98.83   0.10
                2   142     1.84667 52.45333    0.00    52.60   47.40
                3   182     1.91950 52.17750    9.48    90.38   0.14
                4   124     1.88333 52.70833    0.00    98.80   1.20
                5   2807    1.91050 51.45000    2.05    97.91   0.05",
             header = TRUE)
my_positions <-
  read.table(text = 
               "haul    DecStartLat DecStartLong
                1993H_2 55.23983    -5.512830
                2794H_1 55.26670    -5.516700
                1993H_1 55.27183    -5.521330",
             header = TRUE)

欧氏距离（使用RANN包）

library(RANN)
# For each point in my_positions, find the nearest neighbor from BGS_df:
# Give X and then Y (longtitude and then latitude)
# Note that argument k sets the number of nearest neighbours, here 1 (the closest)
closest_RANN <- RANN::nn2(data = BGS_df[, c("X", "Y")], 
                          query = my_positions[, c("DecStartLong", "DecStartLat")], 
                          k = 1)
results_RANN <- cbind(my_positions[, c("haul", "DecStartLong", "DecStartLat")],
                      BGS_df[closest_RANN$nn.idx, ])
results_RANN
#        haul DecStartLong DecStartLat NUM       X        Y GRAV SAND MUD
# 4   1993H_2     -5.51283    55.23983 124 1.88333 52.70833    0 98.8 1.2
# 4.1 2794H_1     -5.51670    55.26670 124 1.88333 52.70833    0 98.8 1.2
# 4.2 1993H_1     -5.52133    55.27183 124 1.88333 52.70833    0 98.8 1.2

大圆距离（使用geosphere包装）

library(geosphere)
# Compute matrix of great circle distances
dist_mat <- geosphere::distm(x = BGS_df[, c("X", "Y")],
                             y = my_positions[, c("DecStartLong", "DecStartLat")],
                             fun = distHaversine) # can try other distances
# For each column (point in my_positions) get the index of row of min dist
# (corresponds to row index in BGS_df) 
BGS_idx <- apply(dist_mat, 2, which.min)
results_geo <- cbind(my_positions[, c("haul", "DecStartLong", "DecStartLat")],
                     BGS_df[BGS_idx, ])
identical(results_geo, results_RANN) # here TRUE, but not always expected