Avoiding multiple for-loops in R to calculate a matrix

Question

So in the course of generating some fake data to answer a map question, I found myself writing the following:

# 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

Now obviously those two central for-loops are not as R'ish as I would like. It seems like I should be able to get an mapply or something to do the job, but sadly that returns a list, and does not really do what I want. The other applys don't seem to do the right thing either.

What am I missing here?

Answer 1

You should try to use matrix algebra. No need to use any functions from the apply family:

lat <- seq(-90, 90, by = 5)
lon <- seq(-180, 180, by = 10)
1500 * tcrossprod(sin(pi * lat / 180), cos(pi * lon / 180))^2 + 4500

Answer 2

you can use 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's answer seems to be the faster than the other solution. Here some benchmark

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