求R中数据点的某一半径内的点数

Question

我有两个数据集，一个用于医院，另一个用于程序。每个数据集都有纬度和经度坐标。程序要么在医院里，要么在医院外面，尽管如果在医院里给出的坐标不一定精确。我试图在每个医院周围形成一个一定大小的半径，并确定平均有多少个手术点在这个半径内。举个例子，如果我有100家医院和3000家医院，我想在所有医院周围形成一个半径，看看平均有多少家医院属于这个指定的半径。我的初始代码在下面，但我知道这可以更快地完成。编码为R.谢谢！

for(i in 1:NROW(hospitals)){
  hospital <- hospitals[i,]
  radius <- .016

  # find all the procedures that lie in the .016 sized radius from this hospital

  hospital$latitude_low <- hospital$lat - radius
  hospital$longitude_low <- hospital$long - radius
  hospital$latitude_high <- hospital$lat + radius
  hospital$longitude_high <- hospital$long + radius

  in_rad <- procedures[(procedures$long >= hospital$longitude_low & procedures$long <= 
  hospital$longitude_high & procedures$lat <= hospital$latitude_high & procedures$lat >= 
  hospital$latitude_low),]

  num <- NROW(in_rad)
  hospitals[i,]$number_of_procedures <- num
}

Answer 1

当你问一个问题时，你应该总是包括一些例子数据。像这样

lat <- c(-23.8, -25.8)
lon <- c(-49.6, -44.6)
hosp <- cbind(lon, lat)


lat <- c(-22.8, -24.8, -29.1, -28, -20)
lon <- c(-46.4, -46.3, -45.3, -40, -30)
procedures <- cbind(lon, lat)

你的数据是经度/纬度吗？如果是这样，则需要使用适当的方法来计算距离。例如

 library(geosphere)
 dm <- distm(procedures, hosp)

或

 library(raster)
 d <- pointDistance(procedures, hosp, lonlat=TRUE)

两者都计算出从所有程序到所有医院的距离。对于非常大的数据集，这将失败，但是按照您的描述，它应该工作得很好。现在你可以用一个阈值(这里是400,000米)来找出每个医院在这个距离内的程序。

apply(d < 400000, 2, which)
#[[1]]
#[1] 1 2

#[[2]]
#[1] 1 2 3

所以手术1，2和3就在第2医院的距离之内。

如果数据不是经度/纬度，则可以使用

 d <- pointDistance(procedures, hosp, lonlat=FALSE)

Answer 2

这里有几样东西可以改进。首先，你实际上并不是在医院半径0.16的范围内计算程序，而是在0.32 * 0.32单位正方形内以医院为中心计算程序。对于具体的问题，可能不是什么大问题，但实际上，在特定的距离内计算点会更快，就像您实际想要的那样。

其次，您有一个倾向于存储您计算过的任何变量，即使您只使用它们一次。这可以帮助理解代码，但有时效率较低，而且肯定会使代码更长，特别是如果您喜欢使用long_descriptive_variable_names。

第三，最后是子集procedures，然后度量行数，而不仅仅是使用子集本身的长度。

最后(但不太重要)，将结果一次写到一个新列中。您可以一次性使用sapply来完成这一切。

因此，您的代码可以被一些简单得多的代码替换，比如：

hospitals$number_of_procedures <- sapply(1:NROW(hospitals), function(i)
  {
    d <- (procedures$long - hospitals[i,]$long)^2 + (procedures$lat - hospitals[i,]$lat)^2
    length(which(d < 0.16^2))
  })