weighted centroid function

weightedCentres.R

## Function for returning weighted centroid location
## currently hard-coded to return mean, q10, q20, q80 & q90.
weightedCentre <- function(x, y, z) {
  require(matrixStats); require(Hmisc)
  if (anyNA(c(x, y, z))) {
    stop("There are missing values present in x, y, or z")
  }
  if (length(z) != length(x)) {
    stop("Number of weights supplied not equal to number of coordinates")
  }
  x_y = cbind(x, y)
  w = z
  mean = matrixStats::colWeightedMeans(x_y, w = w, na.rm = FALSE)
  q10 = cbind(Hmisc::wtd.quantile(x = x_y[,1], weights = w, probs = 0.10), ## x
              Hmisc::wtd.quantile(x = x_y[,2], weights = w, probs = 0.10)) ## y
  q20 = cbind(Hmisc::wtd.quantile(x = x_y[,1], weights = w, probs = 0.20),
              Hmisc::wtd.quantile(x = x_y[,2], weights = w, probs = 0.20))
  q80 = cbind(Hmisc::wtd.quantile(x = x_y[,1], weights = w, probs = 0.80),
              Hmisc::wtd.quantile(x = x_y[,2], weights = w, probs = 0.80))
  q90 = cbind(Hmisc::wtd.quantile(x = x_y[,1], weights = w, probs = 0.90),
              Hmisc::wtd.quantile(x = x_y[,2], weights = w, probs = 0.90))
  colnames(q10)<-colnames(q20)<-colnames(q80)<-colnames(q90) <- c("x", "y")
  wght_cent <- matrix(rbind(mean, q10, q20, q80, q90), ncol = 2)
  rownames(wght_cent) <- c("mean", "q10", "q20", "q80", "q90")
  return(wght_cent)
}
### EXAMPLE ###
library(sp)
data(meuse)
# x-coords, y-coords, "weights"
df_x <- meuse$x; df_y <- meuse$y; df_w = meuse$cadmium
plot(df_x, df_y, pch = 19, col = heat.colors(20))
points(mean(df_x), mean(df_y), col = "blue", pch = 12) ## geographic centroid
new_cents <- weightedCentre(x = df_x, y = df_y, z = df_w)
## add weighted centroids to map
points(new_cents, col = c("#0000FF", "#8A2BE2", "#008B00", "#4F4F4F", "#8B4513"),
       pch = 16)
legend(x = 178800, y = 333500, legend = c("geo cent", rownames(new_cents)),
       col = c("blue", "#0000FF", "#8A2BE2", "#008B00", "#4F4F4F", "#8B4513"),
       pch = c(12, rep(16,5)))