title: How do I customise the colours of the groups and overlay a contour with values over an NMDS ordination spider graph plot using ggplot2?
body: <p>I've ran an NMDS using <code>vegan</code> on some ecological data showing species frequency and habitat types and produced a spider graph plot showing the differences between the groups(habitats) using some code I found.
  I'd like to be able to change the colours of the habitats to represent a gradient from one main habitat to the next and add a contour over the top showing the change in Shannon diversity using <code>ordisurf</code> with figures on the contours.</p>
  <p>Here is some example data (apologies if this is too large of a subset) named <code>dataNMDS</code>.</p>
  <pre><code>structure(list(Plot_ID = c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 
  1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), 
      Plot_gradient = c(&quot;Grassland&quot;, &quot;Grassland&quot;, &quot;Grassland&quot;, 
      &quot;EGA&quot;, &quot;EGA&quot;, &quot;EGA&quot;, &quot;EGB&quot;, &quot;EGB&quot;, &quot;EGB&quot;, &quot;EGC&quot;, &quot;EGC&quot;, &quot;EGC&quot;, 
      &quot;EHC&quot;, &quot;EHC&quot;, &quot;EHC&quot;, &quot;EHB&quot;, &quot;EHB&quot;, &quot;EHB&quot;, &quot;EHA&quot;, &quot;EHA&quot;, &quot;EHA&quot;, 
      &quot;Heath&quot;, &quot;Heath&quot;, &quot;Heath&quot;), SSD_stdev = c(6.475059845, 7.730873172, 
      5.451752012, 1.854723699, 3.847076812, 1.939071943, 3.382306905, 
      2.332380758, 1.673320053, 2.481934729, 4.955804677, 1.624807681, 
      3.2, 6.337191807, 0.748331477, 0.632455532, 3.521363372, 
      0.8, 1.264911064, 9.583318841, 0.748331477, 1.977473135, 
      6.66765326, 3.404937591), log_mean = c(3.77, 3.77, 3.68, 
      3.65, 3.68, 3.41, 3.64, 3.8, 3.4, 3.65, 3.72, 3.43, 3.65, 
      3.73, 3.36, 3.57, 3.88, 3.38, 3.6, 3.79, 3.38, 3.84, 3.53, 
      3.3), conductivity = c(770L, 555L, 508L, 508L, 553L, 532L, 
      605L, 443L, 454L, 561L, 533L, 554L, 502L, 435L, 523L, 452L, 
      384L, 494L, 502L, 411L, 378L, 461L, 507L, 383L), percentage_c = c(32.144, 
      28.316, 31.498, 9.612, 21.184, 20.764, 14.819, 21.187, 16.619, 
      19.554, 28.565, 23.694, 7.539, 14.12, 23.769, 9.696, 13.022, 
      18.966, 10.053, 15.969, 12.395, 6.013, 17.848, 10.716), percentage_n = c(2.398, 
      2.224, 2.321, 0.799, 1.798, 1.497, 1.192, 1.875, 1.091, 1.618, 
      2.149, 1.679, 0.557, 1.172, 1.361, 0.679, 1.121, 1.325, 0.715, 
      1.368, 0.83, 0.465, 1.394, 0.734), percentage_ic = c(0.047, 
      0.049, 0.046, 0.048, 0.049, 0.048, 0.049, 0.048, 0.05, 0.051, 
      0.048, 0.047, 0.049, 0.051, 0.048, 0.05, 0.048, 0.047, 0.051, 
      0.046, 0.05, 0.05, 0.046, 0.05), Racomitrium_lanuginosum = c(2.4, 
      11.2, 3.2, 52, 32, 24, 16, 40, 0, 48, 0, 16, 88, 8, 100, 
      96, 76, 64, 100, 40, 68, 71, 57, 74), Vaccinium_myrtillus = c(2.4, 
      22.4, 9.6, 28, 24, 20, 32, 20, 28, 16, 44, 60, 52, 52, 44, 
      28, 16, 16, 16, 52, 48, 34.5, 18.5, 48), Agrostis_capillaris = c(37.6, 
      54.4, 11.2, 60, 64, 16, 96, 36, 16, 68, 44, 24, 88, 76, 84, 
      92, 92, 0, 68, 76, 88, 57.5, 35.5, 16.5), Galium_saxatile = c(0.8, 
      24, 18.4, 4, 12, 24, 32, 24, 4, 4, 24, 8, 4, 44, 56, 48, 
      56, 4, 0, 36, 52, 15.5, 23.5, 11.5), Carex_pilulifera = c(31.2, 
      15.2, 25.6, 24, 0, 76, 0, 20, 92, 24, 0, 80, 4, 0, 12, 0, 
      0, 0, 20, 0, 0, 0, 0, 0), Festuca_ovina = c(0, 1.6, 0, 16, 
      0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 12, 0, 0, 44, 0, 20, 7.5, 
      27.5, 12), Juncus_squarrosus = c(68L, 56L, 0L, 44L, 0L, 0L, 
      0L, 0L, 0L, 96L, 0L, 0L, 40L, 0L, 0L, 0L, 0L, 0L, 12L, 0L, 
      0L, 0L, 0L, 0L), Hylocomium_splendens = c(21.6, 32.8, 0.8, 
      28, 0, 0, 68, 56, 0, 16, 24, 0, 8, 8, 0, 12, 12, 0, 0, 4, 
      0, 0, 0, 0), Rhytidiadelphus_loreus = c(70.4, 40, 32.8, 32, 
      48, 40, 56, 64, 76, 12, 52, 40, 0, 24, 0, 8, 24, 24, 0, 44, 
      12, 1, 3, 0)), class = &quot;data.frame&quot;, row.names = c(NA, -24L
  ))
  </code></pre>
  <p>This is the code I've used to run the NMDS and produce the spider graph plot in <code>ggplot2</code>.</p>
  <pre><code>#subset the dataframe on which to base the ordination: the species data (dataframe 1)
  data_1 &lt;- dataNMDS[,9:17]
  
  #Identify the columns that contains the descriptive/environmental data (dataframe 2)
  data_2 &lt;- dataNMDS[,2:8]
  
  
  #ordination by NMDS
  NMDS &lt;- metaMDS(data_1, distance = &quot;bray&quot;, k = 2)
  stressplot(NMDS)
  
  
  #########################
  #Data visualisation
  
  
  #Extract the axes scores
  datascores &lt;- as.data.frame(scores(NMDS))  #extract the site scores
  
  #Add/calculate spider diagram
  scores &lt;- cbind(as.data.frame(datascores), Habitat = data_2$Plot_gradient)
  centroids &lt;- aggregate(cbind(NMDS1, NMDS2) ~ Habitat, data = scores, FUN = mean)
  seg &lt;- merge(scores, setNames(centroids, c('Habitat','oNMDS1','oNMDS2')),
               by = 'Habitat', sort = FALSE)
  
  #NMDS spider graph plot
  ggplot(scores, aes(x = NMDS1, y = NMDS2, colour = Habitat)) +
    geom_segment(data = seg,
                 mapping = aes(xend = oNMDS1, yend = oNMDS2)) + # add spiders
    geom_point(data = centroids, size = 4) +                    # add centroids
    geom_point() +                                              
    coord_fixed()+                                              
    theme_bw()+
    theme(legend.position=&quot;right&quot;,legend.text=element_text(size=10),legend.direction='vertical')
  </code></pre>
  <p>This code produces this image below...</p>
  <p><a href="https://i.stack.imgur.com/H8ZOl.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/H8ZOl.png" alt="enter image description here" /></a></p>
  <p>However, I'd like to change the colour of the habitat types showing a gradient from Grassland to heath.</p>
  <p>I've tried things like...</p>
  <pre><code>scale_fill_manual(values = c(&quot;darkgreen&quot; , &quot;green4&quot; , &quot;green3&quot; , &quot;green2&quot;, &quot;orangered&quot; , &quot;orangered2&quot; , &quot;orangered3&quot; , &quot;orangered4&quot;),
                      breaks = c(&quot;Grassland&quot;,&quot;EGA&quot;,&quot;EGB&quot;,&quot;EGC&quot;,&quot;EHC&quot;,&quot;EHB&quot;,&quot;EHA&quot;,&quot;Heath&quot;))
  </code></pre>
  <p>But my limited R skill causes me to fail though, from what I can tell it seems to relate to the <code>aes (colour = Habitat)</code>.</p>
  <p>I've also tried to add the contours of the Shannon diversity which I've calculated with...
  <code>ordisurf(NMDS, diversity(data_1), add = T)</code>
  and</p>
  <pre><code>shannon &lt;- ordisurf(NMDS, diversity(data_1), add = T)
  shannoncurve &lt;- ordisurfgrid.long(shannon)
  + geom_contour(data=shannoncurve,aes(x=x, y=y, z=z))
  </code></pre>
  <p>Though I can't get that to work either.</p>
  <p>Sorry for the incredibly long-winded and probably simple question but any help would be extremely appreciated.</p>