2015-09-15 22 views
1

我開始在R並試圖讓這個循環來執行。我試圖讓循環計算使用函數(Vincenty的公式)在座標之間的連續距離。 'Distfunc'是該函數的文件。該函數然後由下面的'x'調用。我想要的只是一個數據框或座標之間的距離列表。偉大的任何幫助!不能得到一個R循環來執行

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])) 
} 
+3

它可能不會解決所有問題,但對我來說最明顯的一點是,你應該使用'的(我在1:長度(Radians)){'或'for(1 in seq_along(Radians)){' – Benjamin

+3

@Benjamin:其實,對於(我在seq_along(Radians)){' – tguzella

+0

是的,@tguzella說。 ('i'和'1'之間有很大的區別) – Benjamin

回答

1

好了,沒有你所面臨的問題的一個很好的說明和適當的重複的例子,它是很難提供什麼好的見解。首先,請參閱How to make a great R reproducible example?

有很多事情在你做事的方式上並不清晰。首先,爲什麼將source(...)的結果分配給變量Distfunc

無論如何,這裏是一些代碼,我試圖理解這一點,它運行沒有問題,但它不完全符合你的期望(因爲你沒有提供太多的信息)。特別是,代碼使用Mario Pineda-Krch的函數gcd.vifhttp://www.r-bloggers.com/great-circle-distance-calculations-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 
+0

非常感謝您的幫助! – PharmR