Script for teaching, part of the book Geocompuation with R, the original of which can be found here: https://github.com/Robinlovelace/geocompr/blob/master/code/10-centroid-alg.R

10-centroid-alg.R

# Aim: take a matrix representing a convex polygon, return its centroid,
# demonstrate how algorithms work

# Pre-requisite: an input object named poly_mat with 2 columns representing
# vertices of a polygon, with 1st and last rows identical

# Step 1: create sub-triangles, set-up ------------------------------------

O = poly_mat[1, ] # create a point representing the origin
i = 2:(nrow(poly_mat) - 2)
T_all = lapply(i, function(x) {
  rbind(O, poly_mat[x:(x + 1), ], O)
})


# Step 2: calculate triangle centroids ------------------------------------

C_list = lapply(T_all,  function(x) (x[1, ] + x[2, ] + x[3, ]) / 3)
C = do.call(rbind, C_list)

# Step 3: calculate triangle areas ----------------------------------------

A = vapply(T_all, function(x) {
  abs(x[1, 1] * (x[2, 2] - x[3, 2]) +
        x[2, 1] * (x[3, 2] - x[1, 2]) +
        x[3, 1] * (x[1, 2] - x[2, 2]) ) / 2
}, FUN.VALUE = double(1))

# Step 4: calculate area-weighted centroid average ------------------------

poly_area = sum(A)
print(paste0("The area is: ", poly_area))
poly_centroid = c(weighted.mean(C[, 1], A), weighted.mean(C[, 2], A))

# Step 5: output results --------------------------------------------------

print(paste0(
  "The coordinates of the centroid are: ",
  round(poly_centroid[1], 2),
  ", ",
  round(poly_centroid[2], 2)
))