Last active
July 25, 2019 19:31
-
-
Save benfasoli/b6d973cdf42439c177580dfd4ff87fe8 to your computer and use it in GitHub Desktop.
Normalized inverse distance weighting applied over a grid
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/usr/bin/env Rscript | |
n <- 5 | |
# fill matrix of distances | |
d <- matrix(nrow = n, ncol = n) | |
d[1:n, 1] <- 0:(n-1) | |
d[1, 1:n] <- 0:(n-1) | |
d <- d / (n - 1) | |
for (i in 2:n) { | |
for (j in 2:n) { | |
d[i, j] = sqrt(d[i, 1]^2 + d[1, j]^2) | |
} | |
} | |
# fill matrix of weights | |
w <- d | |
w[ , ] <- NA | |
for (i in 2:(n-1)) { | |
for (j in 2:(n-1)) { | |
r1 <- d[i, j] | |
r2 <- d[(n+1)-i, j] | |
r3 <- d[(n+1)-i, (n+1)-j] | |
r4 <- d[i, (n+1)-j] | |
w[i, j] = (1/r1) / ((1/r1) + (1/r2) + (1/r3) + (1/r4)) | |
} | |
} | |
# ensure the sum of applied weights is always 1 | |
ws <- w | |
ws[ , ] <- NA | |
for (i in 2:(n-1)) { | |
for (j in 2:(n-1)) { | |
w1 <- w[i, j] | |
w2 <- w[(n+1)-i, j] | |
w3 <- w[(n+1)-i, (n+1)-j] | |
w4 <- w[i, (n+1)-j] | |
ws[i, j] = w1 + w2 + w3 + w4 | |
} | |
} | |
# distance from (0, 0) | |
d | |
# [,1] [,2] [,3] [,4] [,5] | |
# [1,] 0.00 0.2500000 0.5000000 0.7500000 1.000000 | |
# [2,] 0.25 0.3535534 0.5590170 0.7905694 1.030776 | |
# [3,] 0.50 0.5590170 0.7071068 0.9013878 1.118034 | |
# [4,] 0.75 0.7905694 0.9013878 1.0606602 1.250000 | |
# [5,] 1.00 1.0307764 1.1180340 1.2500000 1.414214 | |
# weight applied | |
w | |
# [,1] [,2] [,3] [,4] [,5] | |
# [1,] NA NA NA NA NA | |
# [2,] NA 0.4488813 0.3086089 0.2007458 NA | |
# [3,] NA 0.3086089 0.2500000 0.1913911 NA | |
# [4,] NA 0.2007458 0.1913911 0.1496271 NA | |
# [5,] NA NA NA NA NA | |
# sum of applied weights (should be 1) | |
ws | |
# [,1] [,2] [,3] [,4] [,5] | |
# [1,] NA NA NA NA NA | |
# [2,] NA 1 1 1 NA | |
# [3,] NA 1 1 1 NA | |
# [4,] NA 1 1 1 NA | |
# [5,] NA NA NA NA NA | |
df <- data.frame(d = c(d), w = c(w)) | |
with(na.omit(df), plot(d, w, type = 'p')) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment