#'---
#' title: "Macrophyte Niches in MN"
#' author: "Mike Verhoeven, Wes Glisson, Dan Larkin"
#' output: 
#'    html_document:
#'       toc: true
#'       theme: default
#'       toc_depth: 2
#'       toc_float:
#'           collapsed: false
#'---

# header ------------------------------------------------------------------

#' # Preface
#' 
#' DRUM submission version (uses cleaned up dataset and filepaths, adds footer)
#' 
#' This code will pull in dataset and analyze the niche occupancy of aquatic plants in three
#' dimensions-- depth, light availability, phenology--and will produce statistics
#' tests, and visualizations used in the accompanying manuscript.
#' 
#' 
#' ### Load libraries:

# +warning=FALSE, message=FALSE 
  library(knitr)
  library(ggplot2)  
  library(ezknitr)
  library(data.table)
  library(TPD)
  library(gridExtra)
  library(ggnewscale)
  library(GGally)
  library(ngram)
  library(sf)
  library(maps)
  library(rgdal)
  library(ggsn)
  library(moments)




# load in data ------------------------------------------------------------

#' ### Refresh data folder:
  # run this line to update your local copy of project data folder files
  # source("scripts/data_folder_pull.R")

#' ### Load in Data:
  #plant surveys
  plants <- fread(file = "plants_DRUM.csv", drop = 1)
  
# explore dataset ---------------------------------------------------------

  str(plants)
  
  
#' # Data Exploration and Prep
#' 
#'   Ensure complete on all vars we'll be using. 
  
  #whats in the dataset?
  plants[ , .N , POINT_ID ] #353,148 sampled locations
  plants[, .N ,SURVEY_ID.x ] # Surveys
  plants[, .N ,SURVEY_ID.y ] # Surveys with covariates
  
  plants[ , .N , .(DOWLKNUM, YEAR.SURVEY) ] #Lake-years
  plants[ , .N, DOWLKNUM]# Lakes
  plants[ , .N, POINT_ID]# Survey points
  plants[!is.na(TAXON) , .N, OBS_ID]
  surveypoints <- plants[ , length(unique(POINT_ID)), SURVEY_ID.x][ , V1,] #survey points per taxon
  hist(surveypoints, breaks = 250, xlim = c(0,500))
  mean(surveypoints)
  sd(surveypoints)
  skewness(surveypoints)
  kurtosis(surveypoints)
  
  observations <- plants[ , length(unique(OBS_ID)), POINT_ID][ , V1,] #survey points per taxon
  hist(observations)
  mean(observations)
  sd(observations)
  skewness(observations)
  kurtosis(observations)
  
  
  
  #survey point per survey
  
  
  #retain only locs for which we have all vars to estimate on
    plants[is.na(DEPTH_FT) == F &
           is.na(gdd) == F &
           is.na(proplight) == F, .N, POINT_ID] # of sampled locs, 292,824 have all data
  
  #percent coverage by sample location for complete cases
  nrow(plants[is.na(DEPTH_FT) == F &
                is.na(gdd) == F &
                is.na(proplight) == F, .N, POINT_ID])/nrow(plants[ , .N , POINT_ID ])
  
  #percent coverage by species observationssample location for 
  plants[is.na(DEPTH_FT) == F &
                is.na(gdd) == F &
                is.na(proplight) == F, .N,]/plants[ , .N ,]

#' ### Include only the complete cases
  #retain only complete cases
  plants <- plants[is.na(DEPTH_FT) == F &
                     is.na(gdd) == F &
                     is.na(proplight) == F, ,]
  
  plants[ , .N, POINT_ID]# Survey points
  plants[!is.na(TAXON) , .N, OBS_ID]# plant obs

#' ### WX distance from sample lakes
  # grab each lake's distance (take mean w/in lake)
  WXlk_dist <- plants[!is.na(distance) , mean(distance) , DOWLKNUM] # distance of lakes from WX
  # calc distance stats across lakes
  WXlk_dist[ ,c(meandist=mean(V1), maxdist = max(V1), mindist = min(V1), sddist = sd(V1)),  ,]
  

#' ## Survey Effort
#' 
#' Check out the data. We want to have a look at the overall distribution of
#' survey effort, and all of the niche parameters we hope to estimate on 
#'  
  
  # locations sampled (done by including POINT_ID in "by")
  # summary(plants)
  locssamp <-  plants[ , .(Secchi = mean(Secchi_m.mean) , DEPTH_FT = mean(DEPTH_FT), proplight = mean(proplight), gdd = mean(gdd), year = mean(YEAR.SURVEY)) , .(DOWLKNUM, DATASOURCE, SURVEY_DATE, POINT_ID, SURVEY_ID.x) ]
  
  hist(locssamp$proplight)
  hist(locssamp$DEPTH_FT, breaks = seq(0,1000,1), xlim = c(0,30))
  hist(locssamp$gdd, breaks = seq(0,3500,100))
  
  #sample effort by lakes, years, lake year, and surveys 
  locssamp[ , .N, .(DOWLKNUM, year) ]
  locssamp[ , .N, .(DOWLKNUM) ]
  locssamp[ , .N, .(year) ]
  locssamp[ , .N, SURVEY_ID.x]
  

#' ## PDF of sampling:
#' What does the probability distribution fn of our sampling look like?
#'  
  #kernel density plots
  kde_depths <- density(locssamp$DEPTH_FT)
  plot(kde_depths, main = "TPD depth")
  plot(kde_depths, main = "TPD depth", xlim = c(0,30))
  
  
  kde_light <- density(locssamp$proplight)
  plot(kde_light, main = "TPD light")
  
  kde_gdd <- density(locssamp$gdd)
  plot(kde_gdd, main = "TPD gdd")

#' ## Niche Axes Orthogonality
#' What about covariation in our niche axes:

  # plot covariance and correllation coefs (r) from the niche axes (uncomment for visulaization thats computationally expensive)
  ggpairs(locssamp[ ,.(proplight, DEPTH_FT, gdd) ,], lower = list(continuous = "smooth"), ggplot2::aes(fill = NULL))
  cor.test(locssamp[ ,c(proplight) ,], locssamp[ ,c(DEPTH_FT) ,])
  cor.test(locssamp[ ,c(gdd) ,], locssamp[ ,c(proplight) ,])
  cor.test(locssamp[ ,c(gdd) ,], locssamp[ ,c(DEPTH_FT) ,])
  
#' Let's take a better look at those light-depth relationships:
  ggplot(locssamp, aes(DEPTH_FT,proplight))+
           geom_point(alpha = .01)+
           xlim(c(0,23))

  hist(locssamp$Secchi)

  summary(locssamp$Secchi)
  
#' We can see from this that we do have some non-orthogonality in the light-
#' depth relationship (Correlation Coef 0.416). It doesn't seem too high to me, 
#' but by the pearson correllation we can reject the null that there's no 
#' correlation. In additon, we can see that the light-depth relationship is well
#' represented considering that light exponentially decays with depth.Thus, the 
#' variation we see w/ light at depth is due to variation in water clarity-- 
#' which we consider good because it provides us with a range of light vals for 
#' any given depth. I don't see guidance in the lit about bounds on 
#' reasonable covariance for this application, so we are going to run with it.    
#' 
#' ## Plant Occurrences
#' 
#' We want to keep only the species observations for niche models. Clean up the
#' Taxon == NA. These are cases where no veg was found. 
#'  

  plants[is.na(TAXON) == T , .N ,] # no veg found cases
  plants[is.na(TAXON) == F , .N ,] # plant observations
  
  plantocc <- plants[is.na(TAXON) == F , ,]

#' Are the plants occurring at reasonable values of light avialability? Ususally
#' consider submerged plants to be adapted to 2-10% of incoming light.  
  plantocc[proplight < .10 & proplight > 0.02, .N , ]/plantocc[,.N,]
  plantocc[proplight > 0.02, .N , ]/plantocc[,.N,]
  hist(plantocc$proplight)
  
#' ## Depth outliers  
#' Need to trim some of the depth outliers. Let's have a look at what those 
#' outliers look like. 

  #explore outliers
  plantocc[, summary(DEPTH_FT) , ]
  hist(plantocc$DEPTH_FT, breaks = 200)
  hist(plantocc$DEPTH_FT, xlim = c(0,25), breaks = 200)
  boxplot(plantocc$DEPTH_FT, ylim = c(0,75))
  sd(plantocc$DEPTH_FT)
  
#' Lets call everything > 5 Standard Deviations form the mean an outlier for depth (Muthukrishnan et al 2018)
  outlier_cut <- mean(plantocc$DEPTH_FT)+5*sd(plantocc$DEPTH_FT) #conveniently this corresponds to 6.998m
  
  #who gets cut?
  plantocc[DEPTH_FT > outlier_cut |
             DEPTH_FT==0 , .N,]
    plantocc[DEPTH_FT > outlier_cut |
             DEPTH_FT==0, c(1:15) ]
  plantocc[ , .N,]
  #proportion
  plantocc[DEPTH_FT > outlier_cut |
             DEPTH_FT==0 , .N,]/plantocc[ , .N,]
  
  
  plantocc <- plantocc[!DEPTH_FT > outlier_cut , ,]
  plantocc <- plantocc[!DEPTH_FT==0]
  
# occurrence data by species ------------------------------------------------------

#' ## Descriptive Stats by Species
#' 
#' Calculate H Index Values for these species  

#' Nsurveys, Noccurrance, H-index by species
  hind <- plants[ , .N, .(SURVEY_DATE, DATASOURCE, DOWLKNUM, TAXON) ] #n obs per species per survey
  hind$TAXON <- as.factor(hind$TAXON)
  setorder(hind, TAXON, -N)
  hind[, index := 1:.N, by=TAXON]
  hind[N > index , hindex:= .N, by = TAXON]
  hind[ ,  nsurv := .N , by = TAXON]
  hind[is.na(hindex) == T , hindex:= 0,]
  hind[ , nocc:= sum(N) , by = .(TAXON)]
  hind2 <- hind[ , .(hindex=max(hindex)) , by = .(TAXON,nsurv,nocc)]
  setorder(hind2, -nocc)
  
  #label taxonomic level
  hind2[ , TAXON := as.character(TAXON), ]
  # strips brackets from non-species taxa
  hind2[, taxacl := gsub("\\s*\\([^\\)]+\\)","",TAXON),]
  hind2$taxawords <- sapply(hind2$taxacl, wordcount)
  hind2[taxawords>1, species := T,]
  hind2[is.na(species), species:= F,]
  hind2[ , .N , species]
  hind2[ nocc>10000, .N , species]
  
  # check out results:
  hind2[1:75]
  
  ggplot(hind,aes(hindex))+
    geom_density()
  ggplot(hind2, aes(hindex))+
    geom_histogram()
  plot(hind2$nsurv, hind2$hindex)
  plot(hind2$nocc, hind2$hindex)
  plot(hind2$nocc, hind2$nsurv)
  
#' ## Threshold for inclusion
#'  Then we simply choose a cutoff value and press forward. We will use 1000
#'  occurrences, and only species level IDs 

# Set cutoff values for inclusion -----------------------------------------


  # drop at some threshold
  
   ptrm <- plantocc[ TAXON %in% hind2[nocc > 1000 &
                                        species==T, TAXON, ], , ]
  
  ptrm[,.N , TAXON]# by taxon, n obs
  ptrm[,.N ,]# total n obs
  
# TPD ----------------------------------------------------------------

#' # Modeling Niches
#'  
#' Now that we have trimmed our dataset to 1) only locs where plants occurred,
#' 2) excluding observations that are depth outliers, 3) only taxa with more
#' than 1000 occurrences, we can move on to modeling the species' niches.
#' 
#' ## Total niche model

  # calculating probability densities for multiple species
  TPDs_plants <- TPDs(species = ptrm$TAXON, traits = ptrm[,.(scale(DEPTH_FT),scale(proplight),scale(gdd))] )

#' Checkout evaluation grid. Each axis was divided into 50 intervals, creating a
#' 3-dimensional matrix (50^3 = 125k cells) with a crossed set of those 50 part 
#' axes. Each species' TPD ought to sum to one. Lets have a look at those 
#' things.
  head(TPDs_plants$data$evaluation_grid) #first chunk of matrix only varies "V1"
  nrow(TPDs_plants$data$evaluation_grid) #size looks right
  names(TPDs_plants$TPDs) #Estimated for all plants
  sum(TPDs_plants$TPDs$'Zanichellia palustris') 
  
  # Do any not sum to one?
  sapply(TPDs_plants$TPDs, sum)[sapply(TPDs_plants$TPDs, sum) != 1 ] #NOPE


# full TPD dissimilarity table --------------------------------------------

#' ## Niche Dissimilarity
  #dissimilarity across in 3 dimensions
  dissim_TPDs_plants <- dissim(TPDs_plants)
  inv_dissim <- setDT(as.data.frame(dissim_TPDs_plants$populations$dissimilarity[ , c('Myriophyllum spicatum', 'Potamogeton crispus')]), keep.rownames = "Taxon")
  inv_dissim[ , clpdis := `Potamogeton crispus`-`Myriophyllum spicatum` ,]
  
  # relative contribution of abundance in shared space vs. spaces only occupied by that species
  inv_pshared <- setDT(as.data.frame(dissim_TPDs_plants$populations$P_shared[ , c('Myriophyllum spicatum', 'Potamogeton crispus')]), keep.rownames = "Taxon")
  inv_pnonshared <- setDT(as.data.frame(dissim_TPDs_plants$populations$P_non_shared[ , c('Myriophyllum spicatum', 'Potamogeton crispus')]), keep.rownames = "Taxon")
  
  invsumm <- cbind(cbind(inv_dissim,inv_pnonshared[,2:3]),inv_pshared[, 2:3])
  setnames(invsumm, c("Taxon", "MspicatumDissim", "PcrispusDissim", "PClessMS", "MspicatumNonShared", "PcrispusNonShared", "MspicatumShared", "PcrispusShared"))
  
#' Test for proportional contrib of non-shared in CLP v EWM:
  # nonshared prop of niche differentiation in CLP exceeds EWM (# of native species comparisons)
  summary(0>invsumm[ , MspicatumNonShared-PcrispusNonShared , ])
  # proportional contribution of shared exceeds nonshared in each?
  summary(0>invsumm[ , PcrispusShared-PcrispusNonShared , ]) # CLP
  summary(0>invsumm[ , MspicatumShared-MspicatumNonShared , ]) # EWM
  
  # Changing shared and non-shared to absolute dissims (i.e., sum to total dissim rather than relativized to sum to 1)
  invsumm[, MspicatumNonShared := MspicatumNonShared*MspicatumDissim]
  invsumm[, MspicatumShared := MspicatumShared*MspicatumDissim]
  invsumm[, PcrispusNonShared := PcrispusNonShared*PcrispusDissim]
  invsumm[, PcrispusShared := PcrispusShared*PcrispusDissim]

  # Changing column orders and dropping "PClessMS"
  setcolorder(invsumm, c("Taxon", "PcrispusDissim", "PcrispusNonShared", "PcrispusShared", 
                         "MspicatumDissim", "MspicatumNonShared", "MspicatumShared"))
  invsumm[,PClessMS:=NULL]
  
  # Adding, and ordering by, mean dissimilarity of each species to all other species
  sppMeanDissim <- rowSums(dissim_TPDs_plants$populations$dissimilarity)/(dim(dissim_TPDs_plants$populations$dissimilarity)[2]-1)
  invsumm[, sppMeanDissim := sppMeanDissim]
  setcolorder(invsumm, c("Taxon", "sppMeanDissim"))
  setorder(invsumm, sppMeanDissim)
  invsumm[is.na(invsumm)] <- 0
  
  # Export Table 1:
  write.csv(invsumm, file = "invsumm.csv", row.names = FALSE)

  # Creating heatmap-style figure from table
  invsumm.m <- melt(invsumm, id.vars = c("Taxon"))  
  invsumm.m$Taxon <- factor(invsumm.m$Taxon, levels = invsumm$Taxon)
  levels(invsumm.m$variable) <- c(
    "All species\nTotal",
    "P. crispus\nTotal",
    "P. crispus\nNon-shared",
    "P. crispus\nShared",
    "M. spicatum\nTotal",
    "M. spicatum\nNon-shared",
    "M. spicatum\nShared")
  
#' # Dissimilarity Table

  dissim_table <- ggplot(data = invsumm.m, aes(x=variable, y=Taxon, fill=value)) + 
    geom_tile() +
    scale_fill_gradient(low = "white", high = "steelblue", name = "Dissimilarity") +
    geom_text(aes(label = sprintf("%0.2f", invsumm.m$value)), size=3) +
    labs(x="Comparison") +
    scale_x_discrete(expand = c(0, 0), position = "top") +
    theme(axis.ticks = element_blank())+
    theme(axis.text.y = element_text(face= "italic"))+
    theme(legend.position = "none")


  # Write to file:
  tiff(file = "Fig_DissimTable.tif", width=8.5, height=10, units="in", res = 600, compression = "lzw")
  dissim_table
  dev.off()
  
# TPD niche axes one at a time --------------------------------------------
  
#' # Decompose by niche axis
#' We will want to 
#' 1) TPD each axis individually
#' 2) Calculate dissimilarity for each axis
#' 3) Decompose each axis' dissimilarity into shared and non-shared niche space
  
  traits_plants_d1 <- ptrm[, c("gdd")]
  sp_plants <- ptrm$TAXON
  TPDs_plants_d1_a1 <- TPDs(species = sp_plants, traits = traits_plants_d1, alpha= 0.95, n_divisions = 100)
  # plotTPD(TPD = TPDs_plants_d1_a1, nRowCol = c(10,5))
  
  traits_plants_d2 <- ptrm[, c(proplight)]
  sp_plants <- ptrm$TAXON
  TPDs_plants_d2_a2 <- TPDs(species = sp_plants, traits = traits_plants_d2, alpha=0.95, n_divisions = 100)
  # plotTPD(TPD = TPDs_plants_d2_a2, nRowCol = c(10,5) ) 
  ptrm[TAXON=="Potamogeton crispus", summary(proplight)]
  
  
  traits_plants_d3 <- ptrm[, c("DEPTH_FT")]
  sp_plants <- ptrm$TAXON
  TPDs_plants_d3_a3 <- TPDs(species = sp_plants, traits = traits_plants_d3, alpha=0.95, n_divisions = 27)
  # plotTPD(TPD = TPDs_plants_d3_a3, nRowCol = c(10,5))   


# calc dissimilarity and decompose by each axis ------------------------------
#' ## Single niche dissimilarity
  #dissimilarity across for each dimension
  dissim_TPDs_gdd <- dissim(TPDs_plants_d1_a1)
  dissim_TPDs_light <- dissim(TPDs_plants_d2_a2)
  dissim_TPDs_depth <- dissim(TPDs_plants_d3_a3)


# dimensional contribution to dissim --------------------------------------------------

#'  We can make a plot that shows the relative contrib from each niche dim:
#'  
#'  
  # Grab mean dissimilarity of each species to all other species for total and each niche axis
  dimdissim <- data.table(Total = rowSums(dissim_TPDs_plants$populations$dissimilarity)/(dim(dissim_TPDs_plants$populations$dissimilarity)[2]-1))
  dimdissim[ , Taxon:= rownames(dissim_TPDs_plants$populations$dissimilarity)  , ]
  dimdissim[ , GDD := rowSums(dissim_TPDs_gdd$populations$dissimilarity)/(dim(dissim_TPDs_gdd$populations$dissimilarity)[2]-1)]
  # dimdissim[ , gddTaxon:= rownames(dissim_TPDs_gdd$populations$dissimilarity)  , ]
  # dimdissim[!Taxon == gddTaxon, , ]
  dimdissim[ , Light := rowSums(dissim_TPDs_light$populations$dissimilarity)/(dim(dissim_TPDs_light$populations$dissimilarity)[2]-1), ]
  # dimdissim[ , lightTaxon:= rownames(dissim_TPDs_light$populations$dissimilarity)  , ]
  # dimdissim[!Taxon == lightTaxon, , ]
  dimdissim[ , Depth :=  rowSums(dissim_TPDs_depth$populations$dissimilarity)/(dim(dissim_TPDs_depth$populations$dissimilarity)[2]-1), ]
  # dimdissim[ , depthTaxon:= rownames(dissim_TPDs_depth$populations$dissimilarity)  , ]
  # dimdissim[!Taxon == depthTaxon, , ]
  
  #order by total dissim   
  setcolorder(dimdissim, c("Taxon", "GDD", "Light","Depth", "Total"))
  setorder(dimdissim, Total)
  
  # Creating heatmap-style figure from table
  dimdissim.m <- melt(dimdissim, id.vars = c("Taxon"))  
  dimdissim.m$Taxon <- factor(dimdissim.m$Taxon, levels = dimdissim$Taxon)

#' # Dimensional Dissimilarity
  
  full_dimfig <- ggplot(data = dimdissim.m, aes(x=variable, y=Taxon, fill=value)) + 
    geom_tile() +
    scale_fill_gradient(low = "white", high = "steelblue", name = "Dissimilarity") +
    geom_text(aes(label = sprintf("%0.2f", dimdissim.m$value)), size=3) +
    labs(x="Dimensional Dissimilarity") +
    scale_x_discrete(expand = c(0, 0), position = "top") +
    theme(axis.ticks = element_blank())+
    theme(legend.position = "none")+
    theme(axis.text.y = element_text(face=c(ifelse(dimdissim$Taxon == "Myriophyllum spicatum" | 
                                                     dimdissim$Taxon == "Potamogeton crispus", "bold.italic", "italic"))))
  
  ifelse(dimdissim$Taxon == "Myriophyllum spicatum" | 
           dimdissim$Taxon == "Potamogeton crispus", "bold", "plain")
  
  simp_dim <- dimdissim.m[Taxon == "Myriophyllum spicatum" | 
                         Taxon == "Potamogeton crispus", , ]
  
  #Simplified heatmap for only the invaders
  simp_dimfig <- ggplot(data = simp_dim, aes(x=variable, y=Taxon, fill=value)) + 
    geom_tile() +
    scale_fill_gradient(low = "white", high = "steelblue", name = "Dissimilarity") +
    geom_text(aes(label = sprintf("%0.2f", simp_dim$value)), size=3) +
    labs(x="Dimensional Dissimilarity") +
    scale_x_discrete(expand = c(0, 0), position = "top") +
    theme(axis.ticks = element_blank())+
    theme(legend.position = "none")+
    theme(axis.text.y = element_text(face= "italic"))
  
  
  # # Write to file:
  tiff(file = "Fig_DimensFull.tif", width=8.5, height=10, units="in", res = 600, compression = "lzw")
  full_dimfig
  dev.off()

  tiff(file = "Fig_DimensSimp.tif", width=8.5, height=2, units="in", res = 600, compression = "lzw")
    simp_dimfig
  dev.off()


# plot all three axis of interest -----------------------------------------

#' # Niche visualizations
  
  #custom plotting all species in a single plot GDD
  gddTPD <- data.table(matrix(unlist(TPDs_plants_d1_a1$TPDs), ncol = length(TPDs_plants_d1_a1$TPDs), byrow = F))
  setnames(gddTPD, names(TPDs_plants_d1_a1$TPDs))
  
  gddTPD[, gddinterval := as.vector(TPDs_plants_d1_a1$data$evaluation_grid) , ]
  
  gddTPD <- melt(gddTPD, id.vars = "gddinterval", variable.name = "taxon" )
  
  levels(gddTPD$taxon)
  gddTPD$taxon <- factor(gddTPD$taxon,levels(gddTPD$taxon)[c(1:8,10:30,32:34,9,31)])
  
  
  gdd_inv <- ggplot(gddTPD, aes(gddinterval, value))+
    geom_line(aes(color = taxon), size = 1)+
    scale_color_manual( values = c(ifelse(levels(gddTPD$taxon) %in% c("Potamogeton crispus") ==T, "black" , ifelse(levels(gddTPD$taxon) %in% c("Myriophyllum spicatum") ==T, "red" , "gray")))) +
    theme(legend.position = "none") +
    xlab("Growing Degrees Days")+
    xlim(c(0,3000))+
    ylab("Density")
  
  #custom plotting all species in a single plot Proplight
  plTPD <- data.table(matrix(unlist(TPDs_plants_d2_a2$TPDs), ncol = length(TPDs_plants_d2_a2$TPDs), byrow = F))
  setnames(plTPD, names(TPDs_plants_d2_a2$TPDs))
  plTPD[, plinterval := as.vector(TPDs_plants_d2_a2$data$evaluation_grid) , ]
  
  plTPD <- melt(plTPD, id.vars = "plinterval", variable.name = "taxon" )
  #bound probability viz to force null at 0 > x > 1
  plTPD[plinterval < 0 | plinterval > 1, value := 0]
  
  plTPD[taxon=="Potamogeton crispus"]
  levels(plTPD$taxon)
  plTPD$taxon <- factor(plTPD$taxon,levels(plTPD$taxon)[c(1:8,10:30,32:34,9,31)])
  
  #add a color selector for legend
  plTPD[ , colorInv := as.factor(ifelse(taxon == "Potamogeton crispus", "black" , ifelse(taxon=="Myriophyllum spicatum", "red" , "gray"))) , ]
  plTPD$colorInv <- factor(plTPD$colorInv, levels = c("black", "red", "gray"))
  
  
  pl_inv <- ggplot(plTPD, aes(group = taxon, x = plinterval, y = value, color = colorInv))+
    geom_line( size = 1)+
    scale_color_manual(name = '', values=c('black','red', 'gray'),
                       labels = c( expression(italic("Potamogeton crispus")), 
                                   expression(italic("Myriophyllum spicatum")),
                                   "Native species")) +
    theme(legend.position = c(0.75,0.5),
          legend.key = element_blank(),
          legend.background = element_blank(),
          legend.text=element_text(size=12.5),
          legend.text.align = 0)+
    xlab("Proportion of Surface Light Remaining at Sediments")+
    xlim(c(-0.01,1.01))+
    # coord_cartesian(ylim = c(0,0.20))+
    ylab("Density")+
    scale_y_continuous(breaks = c(0.00, 0.10,  0.20), labels = c("0.00", "0.10", "0.20"))
  
  #custom plotting all species in a single plot Depth
  depTPD <- data.table(matrix(unlist(TPDs_plants_d3_a3$TPDs), ncol = length(TPDs_plants_d3_a3$TPDs), byrow = F))
  setnames(depTPD, names(TPDs_plants_d2_a2$TPDs))
  
  depTPD[, depinterval := as.vector(TPDs_plants_d3_a3$data$evaluation_grid) , ]
  
  depTPD <- melt(depTPD, id.vars = "depinterval", variable.name = "taxon" )
  depTPD[depinterval < 0 , value := 0]
  
  
  levels(depTPD$taxon)
  depTPD$taxon <- factor(depTPD$taxon,levels(depTPD$taxon)[c(1:8,10:30,32:34,9,31)])
  
  
  dep_inv <- ggplot(depTPD, aes(depinterval*0.3048, value))+
    geom_line(aes(color = taxon), size = 1)+
    scale_color_manual( values = c(ifelse(levels(depTPD$taxon) %in% c("Potamogeton crispus") ==T, "black" , ifelse(levels(depTPD$taxon) %in% c("Myriophyllum spicatum") ==T, "red" , "gray")))) +
    theme(legend.position = "none")+
    xlab("Depth (m)")+
    xlim( c(-0.05,7))+
    ylab("Density")+
    scale_y_continuous(breaks = c(0.00, 0.20,  0.40), labels = c("0.00", "0.20", "0.40"))

# Make panel by niche dimension plot -----------------------------------------------------------
  
# ## One dimensional vizualizations
  
  # # Figure with 3 axes broken out with overlapping species niches
  fig_3axes <- grid.arrange(gdd_inv, pl_inv, dep_inv)
  tiff("fig_3axes.tiff", res = 600, width = 9, height = 12, units = "in", compression = "lzw")
  plot(fig_3axes) # Make plot
  dev.off()
  

# Model TPD 2 dimensions, 3 ways for all species  ------------------------------------------
  
  #gdd v dep
  TPDs_gdddep <- TPDs(species = ptrm[, TAXON ], traits = ptrm[, .(gdd,DEPTH_FT*0.3048) ] )
  
  #gdd v proplight
  TPDs_gddlight <- TPDs(species = ptrm[, TAXON ], traits = ptrm[, .(gdd,proplight) ] )
  
  #dep v proplight
  TPDs_deplight <- TPDs(species = ptrm[, TAXON ], traits = ptrm[, .(DEPTH_FT*0.3048, proplight) ] )
  
  
# Library of niches --------------------------------------------------

#' Great. That covers plotting the sampling distribution in the background. Now
#' scale this to allow calc and paneling by crossed niche and species.
#' 
#'  ## Sampling "Niche" visualization
#'  This is a way to condier how the sampling was biased across our dataset

  #TPD for full cross of gdd, depth, proplight for sample effort dist
  locssamp[ , surveyed := "Samplous distributionous",]
  TPDs_effort_gd <- TPDs(species = locssamp[, surveyed , ], traits = locssamp[, .(gdd, DEPTH_FT*0.3048) ]  )
  TPDs_effort_gl <- TPDs(species = locssamp[, surveyed , ], traits = locssamp[, .(gdd, proplight) ] )
  TPDs_effort_ld <- TPDs(species = locssamp[, surveyed , ], traits = locssamp[, .(proplight, DEPTH_FT*0.3048) ] )
  
  # effort gdd by dep
  e_gd <- data.table(gdd = TPDs_effort_gd$data$evaluation_grid$gdd,
                               dep = TPDs_effort_gd$data$evaluation_grid$V2 ,
                               tpd = TPDs_effort_gd$TPDs$`Samplous distributionous`)
  e_gd[tpd == 0 , tpd := NA , ]
  
  gd_samps <- ggplot(e_gd, aes(x = gdd, y = dep))+
    geom_tile(data = e_gd[is.na(tpd) == F, ], aes(fill = tpd ))
     # scale_fill_gradient(low = "gray", high = "gray")+
     # new_scale_fill()+
     # geom_tile(data = gd[is.na(tpdp) == F, ], aes(fill = tpdp ))
  
  # effort light by dep
  e_gp <- data.table(gdd = TPDs_effort_gl$data$evaluation_grid$gdd,
                     pl = TPDs_effort_gl$data$evaluation_grid$proplight ,
                     tpd = TPDs_effort_gl$TPDs$`Samplous distributionous`)
  e_gp[tpd == 0 , tpd := NA , ]
  
  gl_samps <- ggplot(e_gp, aes(x = gdd, y = pl))+
    geom_tile(data = e_gp[is.na(tpd) == F, ], aes(fill = tpd ))

  # effort light by dep
  e_pd <- data.table(dep = TPDs_effort_ld$data$evaluation_grid$V2,
                     pl = TPDs_effort_ld$data$evaluation_grid$proplight ,
                     tpd = TPDs_effort_ld$TPDs$`Samplous distributionous`)
  e_pd[tpd == 0 , tpd := NA , ]
  
  dl_samps <- ggplot(e_pd, aes(x = pl, y = dep))+
    geom_tile(data = e_pd[is.na(tpd) == F, ], aes(fill = tpd ))
  
  fig_sampleeffort <- grid.arrange(gd_samps, gl_samps, dl_samps, ncol = 3)
  
  # write to file
  # tiff("fig_sampledist.tiff", res = 600, width = 12, height = 4, units = "in", compression = "lzw")
  # plot(fig_sampleeffort) # Make plot
  # dev.off()


#' Okay, so we have a way to plot what we want, now we have to build a table
#' from which we can plot all of the species data. This need to be by crossed set
#' 
#' ## Library of Niches Visualizations
  

# library of niches gdd x depth -------------------------------------------
  
  # gdd x depth
  #effort data
  e_gd <- data.table(gdd = TPDs_effort_gd$data$evaluation_grid$gdd,
                     dep = TPDs_effort_gd$data$evaluation_grid$V2 ,
                     tpd = TPDs_effort_gd$TPDs$`Samplous distributionous`)
  #make species TPDs into data.table
  setDT(TPDs_gdddep$TPDs)
  speciesdat <- data.table(gdd = TPDs_gdddep$data$evaluation_grid$gdd,
                           dep = TPDs_gdddep$data$evaluation_grid$V2,
                           TPDs_gdddep$TPDs
                           )
  #make long format
  speciesdat.m <- melt(speciesdat, id.vars = c('gdd', 'dep' ))
  #plot:
  lib_GDDDepth <- ggplot(e_gd, aes(x = gdd, y = dep))+
    geom_tile(data = e_gd[tpd > 0 , ], aes(fill = tpd ))+
    scale_fill_gradient(low = "gray", high = "gray")+
    theme(legend.position = NULL)+
    new_scale_fill()+
    geom_tile(data = speciesdat.m[value > 0, ], aes(fill = value, group = variable ))+
    facet_wrap(~ variable, ncol = 5 )+
    theme(legend.position = NULL)+
    xlab("Growing degree days")+
    ylab("Depth (m)")+
    theme(strip.text = element_text(face = "italic"))
  
  #plot only the ususal suspects
  inv_GDDDepth <- ggplot(e_gd, aes(x = gdd, y = dep))+
    geom_tile(data = e_gd[tpd > 0 , ], aes(fill = tpd ))+
    scale_fill_gradient(low = "gray", high = "gray")+
    theme(legend.position = NULL)+
    new_scale_fill()+
    geom_tile(data = speciesdat.m[value > 0 &
                                    variable %in% c("Potamogeton crispus", "Myriophyllum spicatum"), ], aes(fill = value))+
    facet_wrap(~ variable, ncol = 2 )+
    # scale_fill_gradient(low = "lightblue" , high = "steelblue", name = "Dissimilarity") +
    theme(legend.position = "none")+
    xlab("Growing degree days")+
    ylab("Depth (m)")+
    theme(strip.text = element_text(face = "italic"))

# library of niches gdd x proplight -------------------------------------------
  
  # gdd x light
  #effort data
  e_gl <- data.table(gdd = TPDs_effort_gl$data$evaluation_grid$gdd,
                     light = TPDs_effort_gl$data$evaluation_grid$proplight ,
                     tpd = TPDs_effort_gl$TPDs$`Samplous distributionous`)
  #make species TPDs into data.table
  setDT(TPDs_gddlight$TPDs)
  speciesdat <- data.table(gdd = TPDs_gddlight$data$evaluation_grid$gdd,
                           light = TPDs_gddlight$data$evaluation_grid$proplight,
                           TPDs_gddlight$TPDs
  )
  #make long format
  speciesdat.m <- melt(speciesdat, id.vars = c('gdd', 'light' ))
  #plot:
  lib_GDDLight <- ggplot(e_gl, aes(x = gdd, y = light))+
    geom_tile(data = e_gl[tpd > 0 , ], aes(fill = tpd ))+
    scale_fill_gradient(low = "gray", high = "gray")+
    theme(legend.position = 'none')+
    new_scale_fill()+
    geom_tile(data = speciesdat.m[value > 0, ], aes(fill = value, group = variable ))+
    facet_wrap(~ variable, ncol = 5 )+
    theme(legend.position = NULL)+
    xlab("Growing degree days")+
    ylab("Proportion of surface irradiance\nremaining at substrate")+
    theme(strip.text = element_text(face = "italic"))
  
  #plot only the ususal suspects
  inv_GDDLight <- ggplot(e_gl, aes(x = gdd, y = light))+
    geom_tile(data = e_gl[tpd > 0 , ], aes(fill = tpd ))+
    scale_fill_gradient(low = "gray", high = "gray")+
    theme(legend.position = NULL)+
    new_scale_fill()+
    geom_tile(data = speciesdat.m[value > 0 &
                                    variable %in% c("Potamogeton crispus", "Myriophyllum spicatum"), ], aes(fill = value))+
    facet_wrap(~ variable, ncol = 2 )+
    # scale_fill_gradient(low = "steelblue", high = "red", name = "Dissimilarity") +
    theme(legend.position = "none")+
    xlab("Growing degree days")+
    ylab("Proportion of surface irradiance\nremaining at substrate")+
    theme(strip.text = element_text(face = "italic"))


# library of niches depth x light -----------------------------------------

  # depth x light
  #effort data
  e_dl <- data.table(dep = TPDs_effort_ld$data$evaluation_grid$V2,
                     light = TPDs_effort_ld$data$evaluation_grid$proplight ,
                     tpd = TPDs_effort_ld$TPDs$`Samplous distributionous`)
  #make species TPDs into data.table
  setDT(TPDs_deplight$TPDs)
  speciesdat <- data.table(dep = TPDs_deplight$data$evaluation_grid$V1,
                           light = TPDs_deplight$data$evaluation_grid$proplight,
                           TPDs_deplight$TPDs
  )
  #make long format
  speciesdat.m <- melt(speciesdat, id.vars = c('dep', 'light' ))
  #plot:
  lib_DepthLight <- ggplot(e_dl, aes(x = dep, y = light))+
    geom_tile(data = e_dl[tpd > 0 , ], aes(fill = tpd ))+
    scale_fill_gradient(low = "gray", high = "gray")+
    theme(legend.position = 'none')+
    new_scale_fill()+
    geom_tile(data = speciesdat.m[value > 0, ], aes(fill = value, group = variable ))+
    facet_wrap(~ variable, ncol = 5 )+
    theme(legend.position = NULL)+
    ylab("Proportion of surface irradiance\nremaining at substrate")+
    xlab("Depth (m)")+
    theme(strip.text = element_text(face = "italic"))
  
  #plot only the ususal suspects
  inv_DepthLight <- ggplot(e_dl, aes(x = dep, y = light))+
    geom_tile(data = e_dl[tpd > 0 , ], aes(fill = tpd ))+
    scale_fill_gradient(low = "gray", high = "gray")+
    theme(legend.position = NULL)+
    new_scale_fill()+
    geom_tile(data = speciesdat.m[value > 0 &
                                    variable %in% c("Potamogeton crispus", "Myriophyllum spicatum"), ], aes(fill = value))+
    facet_wrap(~ variable, ncol = 2 )+
    # scale_fill_gradient(low = "steelblue", high = "red", name = "Dissimilarity") +
    theme(legend.position = "none")+
    ylab("Proportion of surface irradiance\nremaining at substrate")+
    xlab("Depth (m)")+
    theme(strip.text = element_text(face = "italic"))
  
  

  

# invader niche plots -----------------------------------------------------

  fig_invaderNiche <- grid.arrange(inv_GDDDepth, inv_GDDLight, inv_DepthLight)  
  # Write to file:
  tiff(file = "Fig_InvaderNiches.tif", width=8.5, height=11, units="in", res = 600, compression = "lzw")
  grid.arrange(inv_GDDDepth, inv_GDDLight, inv_DepthLight)
  dev.off()


# library of niches dimensional -------------------------------------------------------

  # lib_GDDDepth
  # lib_GDDLight
  # lib_DepthLight
  
  # Write to file:
  tiff(file = "Fig_LibraryGDDDepth.tif", width=10, height=13, units="in", res = 600, compression = "lzw")
  lib_GDDDepth
  dev.off()
  tiff(file = "Fig_LibraryGDDLight.tif", width=10, height=13, units="in", res = 600, compression = "lzw")
  lib_GDDLight
  dev.off()
  tiff(file = "Fig_LibraryDepthLight.tif", width=10, height=13, units="in", res = 600, compression = "lzw")
  lib_DepthLight
  dev.off()

# stat diff in niche dimensional diffs between clp & ewm -------------------------------


#' # Statistical tests referenced in text  
  
  mean_traits <- plantocc[ ,.(gdd_hat = mean(gdd), proplight_hat = mean(proplight) , depthmeter_hat = mean(DEPTH_FT*0.3048) ) , by = TAXON ]  
  
  community_means <- plantocc[!TAXON %in% c('Myriophyllum spicatum', 'Potamogeton crispus'),.(gdd_hat = mean(gdd), proplight_hat = mean(proplight) , depthmeter_hat = mean(DEPTH_FT*0.3048) ) ]
  
  invader_means <- plantocc[TAXON %in% c('Myriophyllum spicatum', 'Potamogeton crispus'),.(gdd_hat = mean(gdd), proplight_hat = mean(proplight) , depthmeter_hat = mean(DEPTH_FT*0.3048) ), TAXON ]
  
  inv_occ <- plantocc[ TAXON %in% c('Myriophyllum spicatum', 'Potamogeton crispus'),.(TAXON,gdd, proplight, depm=DEPTH_FT*0.3048 )]
  
  inv_occ[ , gdd_cs := scale(gdd) ,] 
  inv_occ[ , proplight_cs := scale(proplight) ,] 
  inv_occ[ , depm_cs := scale(depm) ,] 
  
  #' comparing TPDs by axis:
  
  t.test(inv_occ[TAXON == 'Myriophyllum spicatum', gdd ], inv_occ[TAXON == 'Potamogeton crispus', gdd ])
  t.test(inv_occ[TAXON == 'Myriophyllum spicatum', proplight ], inv_occ[TAXON == 'Potamogeton crispus', proplight ])
  t.test(inv_occ[TAXON == 'Myriophyllum spicatum', depm ], inv_occ[TAXON == 'Potamogeton crispus', depm ])
  #with z-transformed vars
  t.test(inv_occ[TAXON == 'Myriophyllum spicatum', gdd_cs ], inv_occ[TAXON == 'Potamogeton crispus', gdd_cs ])
  t.test(inv_occ[TAXON == 'Myriophyllum spicatum', proplight_cs ], inv_occ[TAXON == 'Potamogeton crispus', proplight_cs ])
  t.test(inv_occ[TAXON == 'Myriophyllum spicatum', depm_cs ], inv_occ[TAXON == 'Potamogeton crispus', depm_cs ])
  
  #Kolmogorov–Smirnov test
  ks.test(inv_occ[TAXON == 'Myriophyllum spicatum', gdd ], inv_occ[TAXON == 'Potamogeton crispus', gdd ])
  ks.test(inv_occ[TAXON == 'Myriophyllum spicatum', proplight ], inv_occ[TAXON == 'Potamogeton crispus', proplight ])
  ks.test(inv_occ[TAXON == 'Myriophyllum spicatum', depm ], inv_occ[TAXON == 'Potamogeton crispus', depm ])
  #with z-transformed vars
  ks.test(inv_occ[TAXON == 'Myriophyllum spicatum', gdd_cs ], inv_occ[TAXON == 'Potamogeton crispus', gdd_cs ])
  ks.test(inv_occ[TAXON == 'Myriophyllum spicatum', proplight_cs ], inv_occ[TAXON == 'Potamogeton crispus', proplight_cs ])
  ks.test(inv_occ[TAXON == 'Myriophyllum spicatum', depm_cs ], inv_occ[TAXON == 'Potamogeton crispus', depm_cs ])
  
#' Now testing for who has more cases of non-shared exceeding  
  

# map of data -------------------------------------------------------------
#' # Map the data
#' Here's all the data we started from:
  #plant surveys
  plants_map <- fread(file = "plants_DRUM.csv", drop = 1)
  
  locs <- plants_map[is.na(TAXON) == F , .(.N, unique(center_utm), unique(center_u_1)), by = DOWLKNUM]
  weatherlocs <- plants_map[is.na(TAXON) == F , .(.N, longitude = unique(Long_WX), latitude = unique(Lat_WX)), by = STATION.NAME]
  
  #drop added plants layer
  rm(plants_map)
  summary(locs)
  #drop waterbodies w/o any location data (these are NWI numbered wetlands)
  locs <- na.omit(locs, cols = c("V2","V3"))
  locsutm <- SpatialPointsDataFrame(locs[,.(V2,V3),], locs, proj4string=CRS("+proj=utm +zone=15N +datum=WGS84"))  
  locsgeo <- spTransform(locsutm, CRS("+proj=longlat +datum=WGS84"))
  str(locsgeo)
  locs_sf <- st_as_sf(locsgeo)
  summary(locs_sf$N)
  locs_sf$bin <- cut(locs_sf$N, breaks = c(0,50, 100, 250, 500, 1000, 5000, 10000, 15000), 
                     labels = c("50", "100", "250", "500", "1000", "5000", "10000", "12000") )
  #drop WX station NAs:
  weatherlocs <- na.omit(weatherlocs, cols = c("longitude","latitude"))
  wxlocsutm <- SpatialPointsDataFrame(weatherlocs[,.(longitude, latitude),], weatherlocs, proj4string=CRS("+proj=longlat +zone=15N +datum=WGS84"))  
  wxlocsgeo <- spTransform(wxlocsutm, CRS("+proj=longlat +datum=WGS84"))
  str(wxlocsgeo)
  wxlocs_sf <- st_as_sf(wxlocsgeo)
  
  
  
  #other data for fig
  usa <- map_data("usa")
  canada <- map_data("world", region = "canada")
  states <- map_data("state")
   mn_df <- subset(states, region == c("minnesota"))

  
   # MN DNR Ecological Classification System retrievd from: https://gisdata.mn.gov/dataset/geos-ecological-class-system
   mn_er <- st_read("ecs_provinces_of_mn_v99a.shp")
  
  #TNC ecoregions from resilient sites analysis retrievd from: https://gisdata.mn.gov/dataset/env-resilient-sites-tnc
  tnc_er<- st_read("Ecoregions.shp")
  tnc_er <- tnc_er[c(2,3,4,8) , ]
  
  #main
    study_map <- ggplot(locs_sf)+
      geom_sf(data = tnc_er, aes(fill = ECO_NAME), alpha = 0.2 )+
    #geom_sf(data = mn_er, aes(fill = PROVNAME), alpha = 0.2 )+
    scale_fill_viridis_d()+
    geom_sf( size = 3, shape = 16, colour = "black",  alpha = 0.3)+
    geom_sf(data = wxlocs_sf, shape = 10 , colour = "red", alpha = 1.0, size = 5)+
    geom_polygon(data = states, aes(x=long, y = lat, group = group), color = "black", alpha = .0, lwd = .75)+
    coord_sf(xlim = c(-98, -87.5),  ylim = c(43, 49.5))+
    theme(text = element_text(size=15))+
    ylab("Latitude")+
    xlab("Longitude")+
    theme(legend.position = c(.755,.29),
          legend.background = element_rect(fill = "white"))+
    guides(fill=guide_legend(title="Ecoregion"))
  
  #inset
  inset <- ggplotGrob(
    ggplot() +
      geom_polygon(data = usa,
                   aes(x=long, y = lat, group = group),
                   fill = "gray90", color = "gray50", size = 0.3) +
      # geom_polygon(data = canada, aes(x=long, y = lat, group = group),
                   # fill = "gray90", color = "gray50", size = 0.3) +
      coord_fixed(xlim = c(-110, -55),  ylim = c(25, 60), ratio = 1.2) +
      geom_polygon(data = mn_df, aes(x = long, y = lat), color = "black", size = 0.3,
                   fill = "gray40") +
      geom_polygon(data = states, aes(x=long, y = lat, group = group), color = "black", alpha = .0, lwd = .75)+
      theme_bw() +
      theme(axis.title = element_blank(), axis.ticks = element_blank(), axis.text = element_blank())+
      coord_map("polyconic")
  )    
  
  #combine
  g3 <- study_map +
    north(mn_df, anchor = c(x = -88.15, y = 44)) +
    annotation_custom(grob = inset, xmin = -91.5, xmax = -87.75,
                     ymin = 42.7)+
    scalebar(mn_df, dist = 50, dist_unit = "km",
             transform = TRUE, model = "WGS84", st.size = 3)
  
  # figure construction -----------------------------------------------------
  #Map
  tiff("Map.tiff", res = 600, width = 9, height = 12, units = "in", compression = "lzw", type = "cairo")
  plot(g3) # Make plot
  dev.off()
  
  # footer ------------------------------------------------------------------
  #' ## Document footer 
  #' 
  #' Session Information:
  #+ sessionInfo
  sessionInfo()