Why do my fitted values from ordisurf not match actual variable values?

Question

Here is my matrix for the NMDS ordination (PA.mm.5p.rel).

structure(c(0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 
1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 
1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 
0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 
1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 
0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 
0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 
0, 0, 0, 0, 0, 0, 1), .Dim = c(45L, 8L), .Dimnames = list(c("3", 
"4", "5", "7", "8", "10", "11", "12", "13", "15", "16", "18", 
"19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", 
"30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "40", 
"41", "42", "43", "44", "45", "46", "47", "48", "49", "50", "51"
), c("Onisimus", "Amphipods", "Krill", "Copepods", "Fish", "Gammarus", 
"Arctic.cod", "Sand.Lance")))

Here is my environmental/individual data (PA_Info).

structure(list(Date = structure(c(3L, 3L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 9L, 
10L, 12L, 14L, 15L, 16L, 16L, 16L, 16L, 17L, 18L, 18L, 19L, 19L, 
19L, 20L, 20L, 20L, 21L, 21L, 22L, 23L, 25L), .Label = c("2017-07-19", 
"2017-07-25", "2017-08-07", "2017-08-08", "2017-08-09", "2017-08-12", 
"2017-08-14", "2017-08-15", "2017-08-22", "2017-09-01", "2017-09-04", 
"2018-07-20", "2018-07-23", "2018-07-26", "2018-07-29", "2018-07-30", 
"2018-08-01", "2018-08-02", "2018-08-04", "2018-08-06", "2018-08-08", 
"2018-08-14", "2018-08-19", "2018-08-20", "2018-08-28", "2019-08-08", 
"2019-08-10", "2019-08-11", "2019-08-12", "2019-08-13", "2019-08-14"
), class = "factor"), Floy = c(NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 2019L, 
2014L, NA, 2060L, NA, 2575L, NA), FishID = structure(c(16L, 15L, 
25L, 21L, 23L, 22L, 20L, 19L, 18L, 24L, 28L, 26L, 29L, 27L, 30L, 
8L, 5L, 7L, 4L, 10L, 3L, 11L, 12L, 13L, 31L, 33L, 34L, 35L, 36L, 
37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L, 46L, 47L, 48L, 49L, 
50L, 51L, 53L), .Label = c("Char 08", "Char 09", "Char 112", 
"Char 116", "Char 121", "Char 123", "Char 126", "Char 129", "Char 130", 
"Char 132", "Char 138", "Char 150", "Char 154", "Char 156", "Char 54", 
"Char 56", "Char 58", "Char 60", "Char 61", "Char 62", "Char 63", 
"Char 65", "Char 66", "Char 68", "Char 69", "Char 86", "Char 87", 
"Char 88", "Char 89", "Char 90", "SSS001", "SSS002", "SSS003", 
"SSS004", "SSS005", "SSS006", "SSS007", "SSS008", "SSS009", "SSS010", 
"SSS011", "SSS012", "SSS013", "SSS014", "SSS015", "SSS016", "SSS017A", 
"SSS017B", "SSS018", "SSS019", "SSS020", "SSS021", "SSS022", 
"SSS23", "SSS24", "SSS25", "SSS26", "SSS27", "SSS28", "SSS29", 
"SSS30", "SSS31", "SSS32", "SSS33", "SSS34", "SSS35", "SSS36", 
"SSS37", "SSS38", "SSS39", "SSS40", "SSS41", "SSS42", "SSS43", 
"SSS44", "SSS45", "SSS46", "SSS47", "SSS48", "SSS49", "SSS50", 
"SSS51", "SSS52", "SSS53", "TCH01", "TCH02", "TCH03", "TCH04", 
"TCH05", "TCH06", "TCH07"), class = "factor"), TL = c(71, 33.2, 
35.5, 59.5, 53.7, 65.3, 54.5, 57.5, 61.5, 60, 56.5, 53.5, 55, 
57, 61, 38.5, 65.2, NA, NA, 59, 29.5, 17, 79.5, 81.3, 90, NA, 
67, 73.2, 86, 85, 63.2, 72.4, 65.4, 66.8, 67, 59, 63.2, 79, 75, 
90.8, 76, 41.4, 82, 83, 77.6), FL = structure(c(55L, 6L, 7L, 
25L, 16L, 40L, 19L, 22L, 31L, 27L, 21L, 18L, 20L, 22L, 29L, 8L, 
43L, NA, NA, 25L, 5L, 3L, 67L, 72L, 63L, NA, 44L, 57L, 62L, 71L, 
35L, 53L, 45L, 47L, 42L, 24L, 39L, 65L, 60L, 73L, 61L, 9L, 69L, 
70L, 62L), .Label = c("13.5", "15.8", "16", "16.5", "28.3", "31.5", 
"34.5", "38", "39.4", "40.6", "42", "44", "46.8", "48.6", "50", 
"51.5", "52", "53", "53.2", "53.5", "55.4", "56", "57", "57.2", 
"57.5", "58", "58.3", "58.4", "58.8", "59", "59.2", "59.6", "60", 
"61", "61.1", "61.2", "61.4", "61.8", "62", "62.5", "62.6", "63", 
"63.5", "64", "65", "65.4", "65.6", "66", "66.4", "66.8", "67", 
"67.6", "68", "69", "69.9", "70", "70.1", "70.4", "72.6", "72.8", 
"75", "76", "76 (pcl)", "76.4", "76.8", "77.4", "77.5", "78", 
"79.8", "81.4", "83", "83.8", "89.8"), class = "factor"), Girth = c(NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, 29, NA, NA, NA, NA, 32.4, 
31, 30, NA, NA, 31, 36, 38.6, 43, NA, 20.6, 37, 46, 35.4), Mass = c(3400, 
8.2, 450, 1800, 1500, 2900, 1300, 1300, 2100, 1900, 1650, 1450, 
1625, 1600, 1850, 550, 2700, 3600, NA, 1300, 237.2, 41.2, 4850, 
5100, 6000, NA, 2350, 2700, 4550, 5600, 1750, 2400, 2020, 1800, 
NA, 1940, 2250, 4000, 3720, 5550, 3150, 620, 4350, 6600, 3900
), Sex = structure(c(5L, 3L, 5L, 1L, 1L, 5L, 5L, 5L, 1L, 5L, 
5L, 5L, 5L, 1L, 5L, 5L, 5L, 5L, 1L, 1L, 1L, 3L, 1L, 5L, NA, NA, 
NA, 5L, 5L, 5L, 1L, 5L, 5L, 1L, 5L, 1L, 1L, 5L, 5L, 5L, 1L, 1L, 
1L, 5L, 5L), .Label = c("F", "F?", "I", "J", "M"), class = "factor"), 
    Year = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L), .Label = c("2017", "2018"), class = "factor"), 
    MassCat = c("Medium", "Small", "Small", "Medium", "Small", 
    "Medium", "Small", "Small", "Medium", "Medium", "Medium", 
    "Small", "Medium", "Medium", "Medium", "Small", "Medium", 
    "Medium", NA, "Small", "Small", "Small", "Large", "Large", 
    "Large", NA, "Medium", "Medium", "Large", "Large", "Medium", 
    "Medium", "Medium", "Medium", NA, "Medium", "Medium", "Large", 
    "Large", "Large", "Medium", "Small", "Large", "Large", "Large"
    ), Datect = structure(c(1502078400, 1502078400, 1502251200, 
    1502251200, 1502251200, 1502251200, 1502251200, 1502251200, 
    1502251200, 1502251200, 1502510400, 1502510400, 1502510400, 
    1502510400, 1502510400, 1502683200, 1502683200, 1502683200, 
    1502683200, 1502683200, 1502683200, 1502769600, 1503374400, 
    1504238400, 1532059200, 1532577600, 1532836800, 1532923200, 
    1532923200, 1532923200, 1532923200, 1533096000, 1533182400, 
    1533182400, 1533355200, 1533355200, 1533355200, 1533528000, 
    1533528000, 1533528000, 1533700800, 1533700800, 1534219200, 
    1534651200, 1535428800), class = c("POSIXct", "POSIXt"), tzone = ""), 
    week = structure(c(4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
    4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 7L, 1L, 
    2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 
    4L, 4L, 5L, 5L, 7L), .Label = c("29", "30", "31", "32", "33", 
    "34", "35"), class = "factor"), yday = c(219, 219, 221, 221, 
    221, 221, 221, 221, 221, 221, 224, 224, 224, 224, 224, 226, 
    226, 226, 226, 226, 226, 227, 234, 244, 201, 207, 210, 211, 
    211, 211, 211, 213, 214, 214, 216, 216, 216, 218, 218, 218, 
    220, 220, 226, 231, 240)), row.names = c(10L, 11L, 13L, 14L, 
15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 
28L, 29L, 30L, 33L, 34L, 35L, 36L, 38L, 40L, 41L, 42L, 43L, 44L, 
45L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L, 55L, 56L, 57L, 
58L, 60L), class = "data.frame")

Creating the NMDS and a dataframe of the points.

PA.mm.5p.rel.nmds2<- metaMDS(PA.mm.5p.rel, distance="bray", k = 2, autotransform=FALSE)

allNMDS_DF<-as.data.frame(PA.mm.5p.rel.nmds2$points)

Running an ordisurf and using a function to extract x, y, and z coordinates from the sp.sf object.

sp.sf <- ordisurf(PA.mm.5p.rel.nmds2 ~ PA_Info$TL)

extract.xyz <- function(obj) {
  xy <- expand.grid(x = obj$grid$x, y = obj$grid$y)
  xyz <- cbind(xy, c(obj$grid$z))
  names(xyz) <- c("x", "y", "z")
  return(xyz)
}

This is where I extract the coordinates, however the z values in contour.vals do not match the range of values observed in PA_Info$TL. You can also see this within the legend when you plot and compare the scales in "level" to "total length". Because there are differences I'm worried that I cannot trust this output from the ordisurf or the related p-value (summary(sp.sf)). Does anyone know what could be going wrong?

contour.vals <- extract.xyz(obj = sp.sf, plot = F)
head(contour.vals)
contour.vals <- na.omit(contour.vals)

allNMDS_DFTL <- allNMDS_DF[which(!is.na(PA_Info$TL)),]

p <- ggplot(data = contour.vals, aes(x, y, z = z)) + stat_contour(aes(colour = ..level..)) + 
  coord_cartesian(xlim = c(-2, 2), ylim = c(-1, 1.5)) + theme_bw()+
  geom_point(data = allNMDS_DFTL, aes(MDS1, MDS2, z = 0, shape = PA_Info$Year[which(!is.na(PA_Info$TL))], fill = PA_Info$TL[which(!is.na(PA_Info$TL))]), size=8,shape=21)+
  labs(fill = 'Total Length (cm)') + xlab("Axis 1 (38.2%)")+ylab("Axis 2 (27.9%)")