小贝子编程

R-使用Apply()和用户定义的功能在数据框中计算样本之间的距离



我有一个很大的样本数据集及其地理位置（> 500个样本）。我想使用两种不同的方法来计算所有样品之间的距离（即制作距离矩阵）：分离株之间的圆形距离，以及使用人类迁移的航点之间的圆形距离。前者使用地理包装很简单，我已经在样本数据集之后列出的代码实现了（n = 10）。

test <- structure(list(sample_lon = c(85.1, 101.65, 101.52, 100.5, 77.67, 
78.01, 41.8376999, 136.3070068, 43.0671997, 33.6925011), sample_lat = c(20.95, 
3.11, 3.07, 13.76, 27.49, 27.18, 56.2234001, 49.0937004, 52.3828011, 
-12.3848), Cairo = c(5482676.53004203, 7979572.75185404, 7969372.29645241, 
7300548.40307534, 4535773.29985337, 4577259.89703268, 3035531.48539288, 
8559356.13799981, 2675638.0088102, 4698651.89376045), Istanbul = c(5775286.47479145, 
8414621.15676231, 8406733.79230612, 7555850.43260897, 4744788.34257317, 
4791494.45027787, 1968286.36800993, 7814352.55494176, 1704767.09047958, 
5939243.14729151), PhnomPenh = c(2300013.05256412, 910323.277871997, 
918740.411068992, 487913.521378381, 3302731.94513236, 3256808.88281153, 
7291948.53468698, 5171445.2592504, 7089458.06789339, 8189260.33147658
), AddisAbada = c(5137690.28464086, 6987142.79928228, 6973508.88041575, 
6744246.81620126, 4567281.84480851, 4588137.16563698, 5241878.016939, 
9802623.85641921, 4823403.30909928, 2433156.16935115), UN = c("Southern-Asia", 
"South-Eastern-Asia", "South-Eastern-Asia", "South-Eastern-Asia", 
"Southern-Asia", "Southern-Asia", "Eastern-Europe", "Eastern-Europe", 
"Eastern-Europe", "Eastern-Africa"), continent = c("Asia", "Asia", 
"Asia", "Asia", "Asia", "Asia", "Europe", "Europe", "Europe", 
"Africa")), .Names = c("sample_lon", "sample_lat", "Cairo", "Istanbul", 
"PhnomPenh", "AddisAbada", "UN", "continent"), row.names = c(NA, 
10L), class = "data.frame")
require(geosphere)
test.dvse <- apply(test[c("sample_lon", "sample_lat")], 1, FUN=function(X) distVincentyEllipsoid(X, test[c("sample_lon", "sample_lat")]))

这成功返回了所有分离物之间距离的矩阵。

对于Waypoint分析，我创建了一个用户定义的功能，该功能着眼于样品来自该大洲，然后通过一系列规则来计算隔离物之间的距离。

dCI <- distVincentyEllipsoid(c(31,30), c(28,41)) #Cairo to Istanbul
dCP <- distVincentyEllipsoid(c(31,30), c(104,11)) #Cairo to Phnom Penh
dIP <- distVincentyEllipsoid(c(28,41), c(104,11)) #Istanbul Phnom Penh
calcWayDist <- function(lon1, lat1, con1, lon2, lat2, con2) {
  if (con1 == con2) {
    dd = distVincentyEllipsoid(c(lon1,lat1), c(lon2,lat2))
  } 
  else{
    if (setequal(c("Africa", "Europe"), c(con1, con2)) == TRUE) {
      if (which("Africa" == c(con1, con2)) == 1) {
        afr = distVincentyEllipsoid(c(lon1,lat1), c(31,30)) #dist from isolate to Cairo
        eur = distVincentyEllipsoid(c(lon2,lat2), c(28,41)) #dist from isolate to Istanbul
        dd = afr + dCI + eur
      }
      else {
        afr = distVincentyEllipsoid(c(lon2,lat2), c(31,30)) #dist from isolate to Cairo
        eur = distVincentyEllipsoid(c(lon1,lat1), c(28,41)) #dist from isolate to Istanbul
        dd = afr + dCI + eur
      }
    }
    else if (setequal(c("Africa", "Asia"), c(con1, con2)) == TRUE) {
      s1 = distVincentyEllipsoid(c(lon1,lat1), c(31,30)) #dist from isolate to Cairo
      s2 = distVincentyEllipsoid(c(lon2,lat2), c(31,30)) #dist from isolate to Cairo
      dd = s1 + s2
    }
    else if (setequal(c("Africa", "Melanesia"), c(con1, con2)) == TRUE) {
      if (which("Africa" == c(con1, con2)) == 1) {
        afr = distVincentyEllipsoid(c(lon1,lat1), c(31,30)) #dist from isolate to Cairo
        mel = distVincentyEllipsoid(c(lon2,lat2), c(104,11)) #dist from isolate to Phnom Penh
        dd = afr + dCP + mel
      }
      else {
        afr = distVincentyEllipsoid(c(lon2,lat2), c(31,30)) #dist from isolate to Cairo
        mel = distVincentyEllipsoid(c(lon1,lat1), c(104,11)) #dist from isolate to Phnom Penh
        dd = afr + dCP + mel
      }
    }
    else if (setequal(c("Europe", "Asia"), c(con1, con2)) == TRUE) {
      s1 = distVincentyEllipsoid(c(lon1,lat1), c(28,41)) #dist from isolate to Cairo
      s2 = distVincentyEllipsoid(c(lon2,lat2), c(28,41)) #dist from isolate to Cairo
      dd = s1 + s2
    }
    else if (setequal(c("Europe", "Melanesia"), c(con1, con2)) == TRUE) {
       if (which("Europe" == c(con1, con2)) == 1) {
          eur = distVincentyEllipsoid(c(lon1,lat1), c(28,41)) #dist from isolate to Istanbul
          mel = distVincentyEllipsoid(c(lon2,lat2), c(104,11)) #dist from isolate to Phnom Penh
          dd = eur + dIP + mel
       }
        else {
          eur = distVincentyEllipsoid(c(lon2,lat2), c(28,41)) #dist from isolate to Istanbul
          mel = distVincentyEllipsoid(c(lon1,lat1), c(104,11)) #dist from isolate to Phnom Penh
          dd = eur + dIP + mel 
        }
    }
    else if (setequal(c("Asia", "Melanesia"), c(con1, con2)) == TRUE) {
      s1 = distVincentyEllipsoid(c(lon1,lat1), c(104,11)) #dist from isolate to Phnom Penh
      s2 = distVincentyEllipsoid(c(lon2,lat2), c(104,11)) #dist from isolate to Phnom Penh
      dd = s1 + s2 
    }
  }
 return(dd) 
  }

手动输入值时，此功能可根据需要起作用。例如，

calcWayDist(3.6,15.6,"Africa",51.9,31.9,"Asia")

5225325

但是，当我尝试使用应用功能执行所有可能的分离对之间的计算时，我会收到错误。

test.waypoint <- apply(test, 1, FUN=function(X) calcWayDist(X["sample_lon"], X["sample_lat"], X["continent"], test["sample_lon"], test["sample_lat"], test["continent"]))

错误中的错误（p1）：经度> 360 另外：警告消息： 在if（con1 == con2）{：： 条件的长度> 1，仅使用第一个元素 从：.pointStomatrix（p1）

调用

i接下来尝试"矢量化"我的用户定义函数的输入参数：

calcWayDistbyVec <- function(s1, s2) {
  if (s1[3] == s2[3]) {
    dd = distVincentyEllipsoid(c(as.numeric(s1[1]),as.numeric(s1[2])), c(as.numeric(s2[1]),as.numeric(s2[2])))
  } 
  else{
    if (setequal(c("Africa", "Europe"), c(s1[3], s2[3])) == TRUE) {
      if (which("Africa" == c(s1[3], s2[3])) == 1) {
        afr = distVincentyEllipsoid(c(as.numeric(s1[1]),as.numeric(s1[2])), c(31,30)) #dist from isolate to Cairo
        eur = distVincentyEllipsoid(c(as.numeric(s2[1]),as.numeric(s2[2])), c(28,41)) #dist from isolate to Istanbul
        dd = afr + dCI + eur
      }
      else {
        afr = distVincentyEllipsoid(c(as.numeric(s2[1]),as.numeric(s2[2])), c(31,30)) #dist from isolate to Cairo
        eur = distVincentyEllipsoid(c(as.numeric(s1[1]),as.numeric(s1[2])), c(28,41)) #dist from isolate to Istanbul
        dd = afr + dCI + eur
      }
    }
    else if (setequal(c("Africa", "Asia"), c(s1[3], s2[3])) == TRUE) {
      d1 = distVincentyEllipsoid(c(as.numeric(s1[1]),as.numeric(s1[2])), c(31,30)) #dist from isolate to Cairo
      d2 = distVincentyEllipsoid(c(as.numeric(s2[1]),as.numeric(s2[2])), c(31,30)) #dist from isolate to Cairo
      dd = d1 + d2
    }
    else if (setequal(c("Africa", "Melanesia"), c(s1[3], s2[3])) == TRUE) {
      if (which("Africa" == c(s1[3], s2[3])) == 1) {
        afr = distVincentyEllipsoid(c(as.numeric(s1[1]),as.numeric(s1[2])), c(31,30)) #dist from isolate to Cairo
        mel = distVincentyEllipsoid(c(as.numeric(s2[1]),as.numeric(s2[2])), c(104,11)) #dist from isolate to Phnom Penh
        dd = afr + dCP + mel
      }
      else {
        afr = distVincentyEllipsoid(c(as.numeric(s2[1]),as.numeric(s2[2])), c(31,30)) #dist from isolate to Cairo
        mel = distVincentyEllipsoid(c(as.numeric(s1[1]),as.numeric(s1[2])), c(104,11)) #dist from isolate to Phnom Penh
        dd = afr + dCP + mel
      }
    }
    else if (setequal(c("Europe", "Asia"), c(s1[3], s2[3])) == TRUE) {
      d1 = distVincentyEllipsoid(c(as.numeric(s1[1]),as.numeric(s1[2])), c(28,41)) #dist from isolate to Cairo
      d2 = distVincentyEllipsoid(c(as.numeric(s2[1]),as.numeric(s2[2])), c(28,41)) #dist from isolate to Cairo
      dd = d1 + d2
    }
    else if (setequal(c("Europe", "Melanesia"), c(s1[3], s2[3])) == TRUE) {
       if (which("Europe" == c(s1[3], s2[3])) == 1) {
          eur = distVincentyEllipsoid(c(as.numeric(s1[1]),as.numeric(s1[2])), c(28,41)) #dist from isolate to Istanbul
          mel = distVincentyEllipsoid(c(as.numeric(s2[1]),as.numeric(s2[2])), c(104,11)) #dist from isolate to Phnom Penh
          dd = eur + dIP + mel
       }
        else {
          eur = distVincentyEllipsoid(c(as.numeric(s2[1]),as.numeric(s2[2])), c(28,41)) #dist from isolate to Istanbul
          mel = distVincentyEllipsoid(c(as.numeric(s1[1]),as.numeric(s1[2])), c(104,11)) #dist from isolate to Phnom Penh
          dd = eur + dIP + mel 
        }
    }
    else if (setequal(c("Asia", "Melanesia"), c(s1[3], s2[3])) == TRUE) {
      d1 = distVincentyEllipsoid(c(as.numeric(s1[1]),as.numeric(s1[2])), c(104,11)) #dist from isolate to Phnom Penh
      d2 = distVincentyEllipsoid(c(as.numeric(s2[1]),as.numeric(s2[2])), c(104,11)) #dist from isolate to Phnom Penh
      dd = d1 + d2 
    }
  }
 return(dd) 
  }

这适用于单个计算：

calcWayDistbyVec(c(3.6,15.6,"Africa"), c(41.8,56.2,"Europe"))

[1] 6435307

但是，再次尝试使用时，我会出现错误。

test.waypoint <- apply(test[c("sample_lon", "sample_lat", "continent")], 1, FUN=function(X) calcWayDistbyVec(X, test[c("sample_lon", "sample_lat", "continent")]))

继承中的错误（P，"空间点"）： （列表）无法强制键入" double" 另外：警告消息： 在if（s1 [3] == s2 [3]）{：： 条件的长度> 1，仅使用第一个元素 来自：继承（p，"空间点"）

如果我单点这样做。

test.waypoint.ind <- apply(test[c("sample_lon", "sample_lat", "continent")], 1, FUN=function(X) calcWayDistbyVec(X, c(3.6,15.6,"Africa")))

8702972 11199869 11189668 10520844 7756069 7797556 6438793 12284859 6175273 4535197

如果有人可以识别无效的方法，我将非常感谢它！我已经查看了帖子应用用户定义的功能，将用户定义的功能应用于特定的数据帧列，这帮助我得出了我的第一个距离矩阵计算，但似乎在这里不适用。

我将这个问题传递给了我的大学的R研究小组，而没有堆栈溢出帐户的人提供了解决方案。它不使用apply（）。

calcWayDistMODIFIED <- function(lon1, lat1, con1, lon2, lat2, con2) {
  if (con1 == con2) {
    dd = distVincentyEllipsoid(c(lon1,lat1), c(lon2,lat2))
  } else if (con1 == "Africa" & con2 == "Europe") {
    afr = distVincentyEllipsoid(c(lon1,lat1), c(31,30)) #dist from isolate to Cairo
    eur = distVincentyEllipsoid(c(lon2,lat2), c(28,41)) #dist from isolate to Istanbul
    dd = afr + dCI + eur
  } else if (con1 == "Europe" & con2 == "Africa") {
    afr = distVincentyEllipsoid(c(lon2,lat2), c(31,30)) #dist from isolate to Cairo
    eur = distVincentyEllipsoid(c(lon1,lat1), c(28,41)) #dist from isolate to Istanbul
    dd = afr + dCI + eur
  } else if (setequal(c("Africa", "Asia"), c(con1, con2))) {
    s1 = distVincentyEllipsoid(c(lon1,lat1), c(31,30)) #dist from isolate to Cairo
    s2 = distVincentyEllipsoid(c(lon2,lat2), c(31,30)) #dist from isolate to Cairo
    dd = s1 + s2
  } else if (con1 == "Africa" & con2 == "Melanesia") {
    afr = distVincentyEllipsoid(c(lon1,lat1), c(31,30)) #dist from isolate to Cairo
    mel = distVincentyEllipsoid(c(lon2,lat2), c(104,11)) #dist from isolate to Phnom Penh
    dd = afr + dCP + mel
  } else if (con1 == "Melanesia" & con2 == "Africa") {
    afr = distVincentyEllipsoid(c(lon2,lat2), c(31,30)) #dist from isolate to Cairo
    mel = distVincentyEllipsoid(c(lon1,lat1), c(104,11)) #dist from isolate to Phnom Penh
    dd = afr + dCP + mel
  } else if (setequal(c("Europe", "Asia"), c(con1, con2))) {
    s1 = distVincentyEllipsoid(c(lon1,lat1), c(28,41)) #dist from isolate to Cairo
    s2 = distVincentyEllipsoid(c(lon2,lat2), c(28,41)) #dist from isolate to Cairo
    dd = s1 + s2
  } else if (con1 == "Europe" & con2 == "Melanesia") {
    eur = distVincentyEllipsoid(c(lon1,lat1), c(28,41)) #dist from isolate to Istanbul
    mel = distVincentyEllipsoid(c(lon2,lat2), c(104,11)) #dist from isolate to Phnom Penh
    dd = eur + dIP + mel
  } else if (con1 == "Melanesia" & con2 == "Europe") {
    eur = distVincentyEllipsoid(c(lon2,lat2), c(28,41)) #dist from isolate to Istanbul
    mel = distVincentyEllipsoid(c(lon1,lat1), c(104,11)) #dist from isolate to Phnom Penh
    dd = eur + dIP + mel
  } else if (setequal(c("Asia", "Melanesia"), c(con1, con2))) {
    s1 = distVincentyEllipsoid(c(lon1,lat1), c(104,11)) #dist from isolate to Phnom Penh
    s2 = distVincentyEllipsoid(c(lon2,lat2), c(104,11)) #dist from isolate to Phnom Penh
    dd = s1 + s2
  }
  return(dd)
}
distance_function = function(d) {
  mat = matrix(0, ncol = nrow(d), nrow = nrow(d))
  for(i in 1:nrow(mat)) {
    for(j in 1:nrow(mat)) {
      mat[i,j] = calcWayDistMODIFIED(lon1 = d[i,'sample_lon'], lat1 = d[i,'sample_lat'], con1 = d[i,'continent'], lon2 = d[j,'sample_lon'], lat2 = d[j,'sample_lat'], con2 = d[j,'continent'])
    }
  }
  return(mat)
}

相关内容

最新更新


热门标签：

javascript python java c# php android html jquery c++ css ios sql mysql arrays asp.net json python-3.x ruby-on-rails .net sql-server django objective-c excel regex ruby linux ajax iphone xml vba spring asp.net-mvc database wordpress string postgresql wpf windows xcode bash git oracle list vb.net multithreading eclipse algorithm macos powershell visual-studio image forms numpy scala function api selenium