在 R 中使用纬度和经度计算观测值答案

【问题标题】：Counting Observations Using Latitude and Longitude in R在 R 中使用纬度和经度计算观测值
【发布时间】：2014-03-08 05:27:04
【问题描述】：

菜鸟问题。我无法弄清楚这段代码有什么问题。我试图找出在 400 米半径圆中发生的观察次数。我有每个观察的纬度和经度。我正在尝试创建一个新列，该列将显示 400 米半径范围内的竞争餐厅数量。我包含了我正在使用的代码的数据样本，以及数据帧的 STR。提前致谢。

for (i in seq(nrow(expandedDataFrame2)))
{
  # circle's centre
  xcentre <- df[i,'latitude']
  ycentre <- df[i,'longitude']

  # checking how many restaurants lie within 400 m of the above centre, noofcloserest   column will contain this value
  expandedDataFrame2[i,'noofcloserest'] <- sum(
(expandedDataFrame2[,'latitude'] - xcentre)^2 + 
  (expandedDataFrame2[,'longitude'] - ycentre)^2 
<= 400^2
) - 1

# logging part for deeper analysis
cat(i,': ')

 cat((expandedDataFrame2[,'latitude'] - xcentre)^2 + 
    (expandedDataFrame2[,'longitude'] - ycentre)^2 
  <= 400^2)

cat('\n')

}

示例：

              business_id             restaurantType                                           full_address open       city
1 --5jkZ3-nUPZxUvtcbr8Uw                Greek             1336 N Scottsdale Rd\nScottsdale, AZ 85257    1 Scottsdale
2 --BlvDO_RG2yElKu9XA1_g           Sushi Bars 14870 N Northsight Blvd\nSte 103\nScottsdale, AZ 85260    1 Scottsdale
3 -_Ke8q969OAwEE_-U0qUjw Beer, Wine & Spirits                   18555 N 59th Ave\nGlendale, AZ 85308    0   Glendale
4 -_npP9XdyzILAjtFfX8UAQ           Vietnamese          6025 N 27th Avenue\nSte 24\nPhoenix, AZ 85073    1    Phoenix
5 -2xCV0XGD9NxfWaVwA1-DQ                Pizza                      9008 N 99th Ave\nPeoria, AZ 85345    1     Peoria
6 -3WVw1TNQbPBzaKCaQQ1AQ              Chinese                     302 E Flower St\nPhoenix, AZ 85012    1    Phoenix
   review_count                       name longitude state stars latitude     type      categories1          categories2
1           11 George's Gyros Greek Grill -111.9269    AZ   4.5 33.46337 business       Greek                 <NA>
2           37               Asian Island -111.8983    AZ   4.0 33.62146 business  Sushi Bars             Hawaiian
3            6    Jug 'n Barrel Wine Shop -112.1863    AZ   4.5 33.65387 business        <NA> Beer, Wine & Spirits
4           15          Thao's Sandwiches -112.0739    AZ   3.0 33.44990 business  Vietnamese           Sandwiches
5            4          Nino's Pizzeria 2 -112.2766    AZ   4.0 33.56626 business       Pizza                 <NA>
6          145                China Chili -112.0692    AZ   3.5 33.48585 business     Chinese                 <NA>
  categories3 categories4 categories5 categories6 categories7 categories8 categories9 categories10 isRestaurant Freq
1        <NA>        <NA>        <NA>        <NA>        <NA>        <NA>        <NA>         <NA>         TRUE   66
2     Chinese        <NA>        <NA>        <NA>        <NA>        <NA>        <NA>         <NA>         TRUE   58
3        <NA>        <NA>        <NA>        <NA>        <NA>        <NA>        <NA>         <NA>         TRUE    8
4        <NA>        <NA>        <NA>        <NA>        <NA>        <NA>        <NA>         <NA>         TRUE   44
5        <NA>        <NA>        <NA>        <NA>        <NA>        <NA>        <NA>         <NA>         TRUE  166
6        <NA>        <NA>        <NA>        <NA>        <NA>        <NA>        <NA>         <NA>         TRUE  166
     avgRev  avgStar  duration delta
1 31.32836 3.694030  381 days     0
2 68.62712 3.661017  690 days     0
3 34.33333 3.555556  604 days     1
4 63.22222 3.577778 1916 days     0
5 30.84431 3.482036  226 days     0
6 23.79042 3.535928 2190 days     0

数据的结构是；

str(expandeddataframe2)

'data.frame':   2833 obs. of  28 variables:
 $ business_id   : chr  "--5jkZ3-nUPZxUvtcbr8Uw" "--BlvDO_RG2yElKu9XA1_g" "-_Ke8q969OAwEE_-U0qUjw" "-_npP9XdyzILAjtFfX8UAQ" ...
 $ restaurantType: chr  "Greek" "Sushi Bars" "Beer, Wine & Spirits" "Vietnamese" ...
 $ full_address  : chr  "1336 N Scottsdale Rd\nScottsdale, AZ 85257" "14870 N Northsight Blvd\nSte 103\nScottsdale, AZ 85260" "18555 N 59th Ave\nGlendale, AZ 85308" "6025 N 27th Avenue\nSte 24\nPhoenix, AZ 85073" ...
 $ open          : Factor w/ 2 levels "0","1": 2 2 1 2 2 2 2 2 2 2 ...
 $ city          : chr  "Scottsdale" "Scottsdale" "Glendale" "Phoenix" ...
 $ review_count  : num  11 37 6 15 4 145 255 35 7 7 ...
 $ name          : chr  "George's Gyros Greek Grill" "Asian Island" "Jug 'n Barrel Wine Shop" "Thao's Sandwiches" ...
 $ longitude     : num  -112 -112 -112 -112 -112 ...
 $ state         : chr  "AZ" "AZ" "AZ" "AZ" ...
 $ stars         : num  4.5 4 4.5 3 4 3.5 4.5 4 2.5 4.5 ...
 $ latitude      : num  33.5 33.6 33.7 33.4 33.6 ...
 $ type          : chr  "business" "business" "business" "business" ...
 $ categories1   : chr  "Greek" "Sushi Bars" NA "Vietnamese" ...
 $ categories2   : chr  NA "Hawaiian" "Beer, Wine & Spirits" "Sandwiches" ...
 $ categories3   : chr  NA "Chinese" NA NA ...
 $ categories4   : chr  NA NA NA NA ...
 $ categories5   : chr  NA NA NA NA ...
 $ categories6   : chr  NA NA NA NA ...
 $ categories7   : chr  NA NA NA NA ...
 $ categories8   : chr  NA NA NA NA ...
 $ categories9   : chr  NA NA NA NA ...
 $ categories10  : chr  NA NA NA NA ...
 $ isRestaurant  : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...
 $ Freq          : num  66 58 8 44 166 166 98 35 45 166 ...
 $ avgRev        : num [1:2833(1d)] 31.3 68.6 34.3 63.2 30.8 ...
  ..- attr(*, "dimnames")=List of 1
  .. ..$ : chr  "Greek" "Sushi Bars" "Beer, Wine & Spirits" "Vietnamese" ...
 $ avgStar       : num [1:2833(1d)] 3.69 3.66 3.56 3.58 3.48 ...
  ..- attr(*, "dimnames")=List of 1
  .. ..$ : chr  "Greek" "Sushi Bars" "Beer, Wine & Spirits" "Vietnamese" ...
 $ duration      :Class 'difftime'  atomic [1:2833] 381 690 604 1916 226 ...
  .. ..- attr(*, "units")= chr "days"
 $ delta         : num  0 0 1 0 0 0 0 0 0 0 ...

【问题讨论】：

你能发布一个（最小的）可重现的例子吗？我会在你的中心周围画一个 400 m 的圆圈，并使用 sp 包中的函数（如 points.in.polygon）来计算这个圆圈内的点数。
您不能认真地期望在经纬度上使用毕达哥拉斯距离公式会得到以米为单位的距离，是吗？使用具有特定地球半径的大圆距离函数，或者如果您的点位于一个小区域中，则投影到 UTM 网格系统。

标签： r

【解决方案1】：

所以这是一种方法，它使用包sp 中的函数spDistsN1(...)。调用你的数据框df，

library(sp)

get.dists <- function(i) {
  ref.pt <- with(df[i,],c(longitude,latitude))
  points <- as.matrix(with(df[-i,],cbind(longitude,latitude)))
  dists  <- spDistsN1(points, ref.pt, longlat=T)
  return(length(which(dists<0.4)))
}
df$count <- sapply(1:nrow(df),get.dists)

spDistsN1(points, ref.pt) 计算从ref.pt 到points 中每个点的大圆距离。如果longlat=T 以公里为单位返回距离。所以函数get.dists 生成一个从参考行到每隔一行的距离向量，然后使用length(which(dists<0.4)) 计算有多少小于0.4km。使用sapply(...) 为df 中的每一行调用此函数。

请注意，在您的示例数据集中，没有一家餐厅位于 400m 范围内。

【讨论】：

太棒了！感谢您还解释了函数中发生了什么。帮助像我这样的菜鸟大有帮助。