Dealing with NAs with ordisurf / gg_ordisurf in NMDS analysis

Question

I am trying to plot something like the image below using gg_ordisurf from the ggordiplots package (which seems to work the same as ordisurf of the vegan package`):

However I have NAs values in my continuous variable. And I don't want to remove from the NMDS the sites where I have NAs, or neither predict them based on the others. Isn't there a way to "na.omit/na.rm" them in the function computing the smoothing surface ?

library(vegan)
library(ggordiplots)

# Example data
community_matrix <- matrix(rpois(50, lambda = 5), nrow = 10)
rownames(community_matrix) <- paste0("Site", 1:10)
colnames(community_matrix) <- paste0("Sp", 1:5)

NMDS = metaMDS(community_matrix, k=2, trace = 0, distance = "bray")

coords= as.data.frame(scores(NMDS, "sites")) # coordinates on the NMDS axis
coords$continuous <- c(3, 5, NA, 2, 4, NA, 6, 7, 5, 3) #continuous variable with NAs

# Plot graph
gg_ordisurf(NMDS,
            env.var = coords$continuous, 
            var.label = "Continuous variable")

Answer 1

After a bit of research on gg_ordisurf github source code, the problem was coming from this part of the code function:

  # Calculate default binwidth
  # Can change the binwidth depending on how many contours you want
  if(missing(binwidth)) {
    r <- range(env.var)
    binwidth <- (r[2]-r[1])/15
  } else {
    binwidth = binwidth
  }

so I added na.rm=TRUE in the range function:

  # MODIFED gg_ordisurf function
  if(missing(binwidth)) {
    r <- range(env.var, na.rm = TRUE)
    binwidth <- (r[2]-r[1])/15
  } else {
    binwidth = binwidth
  }

and made my own gg_ordisurf function by copying the source code.

Note: the vegan::ordisurf function was working perfectly fine with this issue.