1

I am starting in R and trying to get this loop to execute. I am trying to get the loop to calculate consecutive distances between coordinates using a function (Vincenty's formula). 'Distfunc' is the file to the function. The function is then called up by 'x' below. All I want then is a data frame or a list of the distances between coordinates. Greatful of any help!

Distfunc <- source("F://Distfunc.R")
for (i in length(Radians)) {
    LatRad1 <- Radians[i,1]
    LongRad1 <- Radians[i,2]
    LatRad2 <- Radians[i+1,1]
    LongRad2 <- Radians[i+1,2]
    x <- gcd.vif(LongRad1, LatRad1, LongRad2, LatRad2) 
    print(data.frame(x[i]))
}
tguzella
  • 1,441
  • 1
  • 12
  • 15
PharmR
  • 260
  • 1
  • 3
  • 12
  • 3
    It may not solve everything, but the most obvious thing to me is that you should use `for (i in 1:length(Radians)){` or `for (1 in seq_along(Radians)){` – Benjamin Sep 15 '15 at 12:07
  • 3
    @Benjamin: Actually, `for (i in seq_along(Radians)){` – tguzella Sep 15 '15 at 12:08
  • Yeah, what @tguzella said. (Big difference between `i` and `1`) – Benjamin Sep 15 '15 at 12:09
  • 4
    @Benjamin, Well, actually that is not a good idea, since the code needs the value at row `i + 1`. Therefore, the for loop should be: `for (i in seq(1, nrow(Radians) - 1)) {` – tguzella Sep 15 '15 at 12:14

1 Answers1

1

Well, without a good description of the problem you are facing and a proper reproducible example it is very difficult to provide any good insight. To start off, see How to make a great R reproducible example?.

There are many things that are not clear in the way you are doing things. First of all, why assign the results of source(...) to the variable Distfunc?

Anyways, here is some code that I put together in trying to understand this; it runs without problems, but it is not clear that it accomplishes what you expect (since you don't provide much information). In particular, the codet uses the implementation for function gcd.vif by Mario Pineda-Krch (http://www.r-bloggers.com/great-circle-distance-calculations-in-r/). The code below is aimed at clarity, since you mention that you are starting in R.

# Calculates the geodesic distance between two points specified by radian latitude/longitude using
# Vincenty inverse formula for ellipsoids (vif)
# By Mario Pineda-Krch (http://www.r-bloggers.com/great-circle-distance-calculations-in-r/)
gcd.vif <- function(long1, lat1, long2, lat2) {

    # WGS-84 ellipsoid parameters
    a <- 6378137 # length of major axis of the ellipsoid (radius at equator)
    b <- 6356752.314245 # ength of minor axis of the ellipsoid (radius at the poles)
    f <- 1/298.257223563 # flattening of the ellipsoid

    L <- long2-long1 # difference in longitude
    U1 <- atan((1-f) * tan(lat1)) # reduced latitude
    U2 <- atan((1-f) * tan(lat2)) # reduced latitude
    sinU1 <- sin(U1)
    cosU1 <- cos(U1)
    sinU2 <- sin(U2)
    cosU2 <- cos(U2)

    cosSqAlpha <- NULL
    sinSigma <- NULL
    cosSigma <- NULL
    cos2SigmaM <- NULL
    sigma <- NULL

    lambda <- L
    lambdaP <- 0
    iterLimit <- 100
    while (abs(lambda-lambdaP) > 1e-12 & iterLimit>0) {
        sinLambda <- sin(lambda)
        cosLambda <- cos(lambda)
        sinSigma <- sqrt( (cosU2*sinLambda) * (cosU2*sinLambda) +
            (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda) )
        if (sinSigma==0) return(0)  # Co-incident points
        cosSigma <- sinU1*sinU2 + cosU1*cosU2*cosLambda
        sigma <- atan2(sinSigma, cosSigma)
        sinAlpha <- cosU1 * cosU2 * sinLambda / sinSigma
        cosSqAlpha <- 1 - sinAlpha*sinAlpha
        cos2SigmaM <- cosSigma - 2*sinU1*sinU2/cosSqAlpha
        if (is.na(cos2SigmaM)) cos2SigmaM <- 0  # Equatorial line: cosSqAlpha=0
        C <- f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha))
        lambdaP <- lambda
        lambda <- L + (1-C) * f * sinAlpha *
            (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)))
        iterLimit <- iterLimit - 1
    }
    if (iterLimit==0) return(NA)  # formula failed to converge
    uSq <- cosSqAlpha * (a*a - b*b) / (b*b)
    A <- 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)))
    B <- uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)))
    deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM^2) -
        B/6*cos2SigmaM*(-3+4*sinSigma^2)*(-3+4*cos2SigmaM^2)))
    s <- b*A*(sigma-deltaSigma) / 1000
    return(s) # Distance in km
}

# Initialize the variable 'Radians' with random data
Radians <- matrix(runif(20, min = 0, max = 2 * pi), ncol = 2)

lst <- list() # temporary list to store the results
for (i in seq(1, nrow(Radians) - 1)) { # loop through each row of the 'Radians' matrix
    LatRad1 <- Radians[i, 1]
    LongRad1 <- Radians[i, 2]
    LatRad2 <- Radians[i + 1, 1]
    LongRad2 <- Radians[i + 1, 2]
    gcd_vif <- gcd.vif(LongRad1, LatRad1, LongRad2, LatRad2)

    # Store the input data and the results
    lst[[i]] <- c(
        latitude_position_1 = LatRad1, 
        longtude_position_1 = LongRad1, 
        latitude_position_2 = LatRad2, 
        longtude_position_2 = LongRad2, 
        GCD = gcd_vif
    )
}
Results <- as.data.frame(do.call(rbind, lst)) # store the input data and the results in a data frame
Community
  • 1
  • 1
tguzella
  • 1,441
  • 1
  • 12
  • 15