Расчет с помощью вложенных циклов в R

Question

У меня есть набор геокодов с тремя столбцами: широта, долгота и кластер. Я вычислил средний центр кластеров и сохранил результаты в двух списках Center_lat и Center_lon.

Теперь я хочу рассчитать расстояние от каждого наблюдения (3000+) до каждого центра кластера (30) с помощью формулы Хаверсина. Чтобы получить матрицу размером 3000 на 30.

Я попытался использовать вложенный цикл, но я получил одинаковое расстояние для всех кластеров. Вот код.

for (i in 1:n){ 
    for (k in 1:c){ 
    lat1=radians(Geocode[i,1]) 
    lon1=radians(Geocode[i,2]) 
    lat2=radians(Center_lat[k,2]) 
    lon2=radians(Center_lon[k,2]) 
    } 
    R <- 3958.756 # Earth mean radius [miles] 
    dist_mat[i,] <- acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) * cos(lon2-lon1)) * R 
    } 

Я также подумываю использовать Lapply для замены вложенного цикла. Но я не уверен, как использовать функцию ... Любая помощь приветствуется.

# Convert to radian 
radians = function(theta=0){return(theta * pi/180)} 

# Calculates the geodesic distance from each property to the center of it's current cluster using the 
# Spherical Law of Cosines (slc) 
get_dist <- function(lat1, lon1, lat2, lon2) { 
    R <- 3958.756 # Earth mean radius [miles] 
    d <- acos(sin(radians(lat1))*sin(radians(lat2)) + 
       cos(radians(lat1))*cos(radians(lat2)) * cos(radians(lon2)-radians(lon1))) * R 
    return(d) # Distance in miles 
} 

dist_mat<-lapply() 

Answer 1

Это тип вычислений, который вы хотите векторизовать в R. Здесь мы используем outer, чтобы сгенерировать все возможные комбинации индексов строк из ваших данных Geocode и Center_x, а затем применить функцию расстояния одним махом.

Во-первых, получить данные проще использовать форму (одну матрицу для местоположений, другой для центров, первые лат столбцов, второй Лон):

# See Data section below for actual data used 

# G <- radians(Geocode) 
# C <- radians(cbind(Center_lat[, 2], Center_lon[, 2]))  
R <- 3958.756 # Earth mean radius [miles] 

Определите функцию, обратите внимание, как мы используем индексы, чтобы посмотреть до реальных координат в G и C, и как функция векторизованную (т.е. нам нужно только назвать это один раз со всеми данными):

my_dist <- function(xind, yind) 
    acos(
    sin(G[xind, 1]) * sin(C[yind, 1]) + 
    cos(G[xind, 1]) * cos(C[yind, 1]) * cos(C[yind, 2] - G[xind, 2]) 
) * R 

и применять его с outer:

DISTS <- outer(seq.int(nrow(G)), seq.int(nrow(C)), my_dist) 

str(DISTS) 
# num [1:3000, 1:30] 4208 6500 8623 7303 3864 ... 

quantile(DISTS) # to make sure stuff is reasonable: 
#  0%  25%  50%  75%  100% 
# 0.000 4107.574 6204.799 8333.155 12422.059  

Это работает около 30 мс в моей системе.

---

данных:

set.seed(1) 
lats <- runif(10000, -60, 60) * pi/180 
lons <- runif(10000, -179, 180) * pi/180 

G.ind <- sample(10000, 3000) 
C.ind <- sample(10000, 30) 

G <- cbind(lats[G.ind], lons[G.ind]) 
C <- cbind(lats[C.ind], lons[C.ind]) 

Answer 2

Похоже, вы хотите записать в матрицу один раз в строке и в столбце, так что вы хотели бы изменить матрицу внутри как для петель, например:

for (i in 1:n){ 
    for (k in 1:c){ 
     lat1=radians(Geocode[i,1]) 
     lon1=radians(Geocode[i,2]) 
     lat2=radians(Center_lat[k,2]) 
     lon2=radians(Center_lon[k,2]) 
     R <- 3958.756 # Earth mean radius [miles] 
     dist_mat[i,k] <- acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2) * cos(lon2-lon1)) * R 
    } 
}