How can I change the colour of my points on my db-RDA triplot in R?

QUESTION: I am building a triplot for the results of my distance-based RDA in R, library(vegan). I can get a triplot to build, but can't figure out how to make the colours of my sites different based on their location. Code below.

#running the db-RDA 
spe.rda.signif=capscale(species~canopy+gmpatch+site+year+Condition(pair), data=env, dist="bray")

#extract % explained by first 2 axes
perc <- round(100*(summary(spe.rda.signif)$cont$importance[2, 1:2]), 2)

#extract scores (coordinates in RDA space)
sc_si <- scores(spe.rda.signif, display="sites", choices=c(1,2), scaling=1)
sc_sp <- scores(spe.rda.signif, display="species", choices=c(1,2), scaling=1)
sc_bp <- scores(spe.rda.signif, display="bp", choices=c(1, 2), scaling=1)

#These are my location or site names that I want to use to define the colours of my points
site_names <-env$site
site_names

#set up blank plot with scaling, axes, and labels
plot(spe.rda.signif,
scaling = 1, 
type = "none", 
frame = FALSE,
xlim = c(-1,1),
ylim = c(-1,1),
main = "Triplot db-RDA - scaling 1",
xlab = paste0("db-RDA1 (", perc[1], "%)"),
ylab = paste0("db-RDA2 (", perc[2], "%)")
)

#add points for site scores - these are the ones that I want to be two different colours based on the labels in the original data, i.e., env$site or site_names defined above. I have copied the current state of the graph
points(sc_si,
pch = 21, # set shape (here, circle with a fill colour)
col = "black", # outline colour
bg = "steelblue", # fill colour
cex = 1.2) # size

Current graph

I am able to add species names and arrows for environmental predictors, but am just stuck on how to change the colour of the site points to reflect their location (I have two locations defined in my original data). I can get them labelled with text, but that is messy.

Any help appreciated!

I have tried separating shape or colour of point by site_name, but no luck.

If you only have a few groups (in your case, two), you could make the group a factor (within the plot call). In R, factors are represented as an integer "behind the scenes" - you can represent up to 8 colors in base R using a simple integer:

set.seed(123)
df <- data.frame(xvals = runif(100),
                 yvals = runif(100),
                 group = sample(c("A", "B"), 100, replace = TRUE))

plot(df[1:2], pch = 21, bg = as.factor(df$group), 
     bty = "n", xlim = c(-1, 2), ylim = c(-1, 2))
legend("topright", unique(df$group), pch = 21, 
       pt.bg = unique(as.factor(df$group)), bty = "n")

If you have more than 8 groups, or if you would like to define your own colors, you can simply create a vector of colors the length of your groups and still use the same factor method, though with a few slight tweaks:

# data with 10 groups
set.seed(123)
df <- data.frame(xvals = runif(100),
                 yvals = runif(100),
                 group = sample(LETTERS[1:10], 100, replace = TRUE))

# 10 group colors
ccols <- c("red", "orange", "blue", "steelblue", "maroon",
                "purple", "green", "lightgreen", "salmon", "yellow")

plot(df[1:2], pch = 21, bg = ccols[as.factor(df$group)], 
     bty = "n", xlim = c(-1, 2), ylim = c(-1, 2))
legend("topright", unique(df$group), pch = 21, 
       pt.bg = ccols[unique(as.factor(df$group))], bty = "n")

For pch just a slight tweak to wrap it in as.numeric:

pchh <- c(21, 22)
ccols <- c("slateblue", "maroon")

plot(df[1:2], pch = pchh[as.numeric(as.factor(df$group))], bg = ccols[as.factor(df$group)], 
     bty = "n", xlim = c(-1, 2), ylim = c(-1, 2))
legend("topright", unique(df$group), 
       pch = pchh[unique(as.numeric(as.factor(df$group)))], 
       pt.bg = ccols[unique(as.factor(df$group))], bty = "n")