Итак, в процессе генерации некоторых поддельных данных для ответа на вопрос о карте я обнаружил, что писал следующее:
# Generate some fake data
lat <- seq(-90, 90, by = 5)
lon <- seq(-180, 180, by = 10)
phi <- matrix(0, nrow = length(lat), ncol = length(lon))
i <- 1
for (l1 in lat) {
j <- 1
for (l2 in lon) {
phi[i, j] <- (sin(pi * l1 / 180) * cos(pi * l2 / 180))^2
j <- j+1
}
i <- i+1
}
phi <- 1500*phi + 4500 # scale it properly
Теперь, очевидно, эти две центральные петли не такие, как R'ish, как хотелось бы. Похоже, я должен получить mapply или что-то сделать для этой работы, но, к сожалению, это возвращает список и не делает то, что я хочу. Другие применяют, похоже, тоже не делают правильных решений.
Что мне здесь не хватает?
Вы должны попытаться использовать матричную алгебру. Нет необходимости использовать какие-либо функции из семейства приложений:
lat <- seq(-90, 90, by = 5)
lon <- seq(-180, 180, by = 10)
1500 * tcrossprod(sin(pi * lat / 180), cos(pi * lon / 180))^2 + 4500
вы можете использовать outer
x = outer(lat, lon, FUN = function(x,y) {(sin(pi * x/180) * cos(pi * y /180))^2})
identical(x * 1500 + 4500, phi)
# [1] TRUE
Ответ NBATrends кажется более быстрым, чем другое решение. Здесь некоторый ориентир
library(microbenchmark)
microbenchmark(within(df, {
phi <- (sin(pi * lat / 180) * cos(pi * lon / 180))^2
phi <- 1500*phi + 4500
}), 1500 * tcrossprod(sin(pi * lat / 180), cos(pi * lon / 180))^2 + 4500, outer(lat, lon, FUN = function(x,y) {(sin(pi * x/180) * cos(pi * y /180))^2}),
((as.matrix(l1)%*%t(as.matrix(l2)))^2) * 1500 + 4500)
Unit: microseconds
expr min lq mean median uq max neval
within(df, { phi <- (sin(pi * lat/180) * cos(pi * lon/180))^2 phi <- 1500 * phi + 4500 }) 255.670 262.0095 270.50948 266.6880 277.7060 385.467 100
1500 * tcrossprod(sin(pi * lat/180), cos(pi * lon/180))^2 + 4500 11.471 12.3770 22.30177 12.9805 13.5850 868.130 100
outer(lat, lon, FUN = function(x, y) { (sin(pi * x/180) * cos(pi * y/180))^2 }) 137.645 139.7590 144.39520 141.5700 145.1925 179.905 100
((as.matrix(l1) %*% t(as.matrix(l2)))^2) * 1500 + 4500 16.301 17.6595 20.20390 19.6215 20.5270 80.294 100
Линейная алгебра может быть проще для вашего приложения, потому что вы просто умножаете элементарные два вектора, что можно сделать через v * u ^ T. В R матричное умножение %*%.
lat <- seq(-90, 90, by = 5)
lon <- seq(-180, 180, by = 10)
l1 <- sin(pi * lat / 180)
l2 <- s(pi * lon/ 180)
# compute the matrix
phi <- as.matrix(l1)%*%t(as.matrix(l2))
# square each element of the matrix
phi <- phi^2
# scale properly
# square each element of the matrix
phi <- 1500*phi + 4500
Зачем прикрепляться к матричной структуре и использовать ее, когда вы можете ее обозначить?
df <- expand.grid(lat = seq(-90, 90, by = 5),
lon = seq(-180, 180, by = 10))
df <- within(df, {
phi <- (sin(pi * lat / 180) * cos(pi * lon / 180))^2
phi <- 1500*phi + 4500
})
Вы всегда можете конвертировать обратно с помощью инструкций здесь.
Используя sapply(), но я бы предпочел решение outer():
#using sapply
phi_1 <-
t(
sapply(lat, function(l1)
sapply(lon, function(l2)(sin(pi * l1 / 180) * cos(pi * l2 / 180))^2))
) * 1500 + 4500
#compare result
identical(phi_1, phi)
# [1] TRUE