## This program spatio-temporally thins iNaturalist milkweed data from specified 
## dates to within 30 km for unique week/year combinations for phenological analysis
library(maptools)
library(rgdal)
library(sp)
library(raster)
library(geosphere)
library(spThin)
#
FunctDirect <- "C:/Users/jamesltracy61/Documents/R/win-library/"
setwd(FunctDirect)
#
# Read in User Defined Functions
source("EvalStatVars_Summary_Function2.R")
source("Round2_Function.R")
#
InDirect <- "K:/Monarch"
setwd(InDirect)
Species <- "Danaus_plexippus"
PopGridName <- paste0("sc2010pop1k")
#
PhenologyDat.df <- data.frame(read.csv(paste0(Species, "_ImmaturePhenologyUnthinnedWeekYearSCStatesRandom.csv"), header = TRUE, sep=',', stringsAsFactors=FALSE))
colnames(PhenologyDat.df)
head(PhenologyDat.df)
tail(PhenologyDat.df)
nrow(PhenologyDat.df)
unique(PhenologyDat.df$Phenological.Stage)
#str(PhenologyDat.df)
# Only keep data after 2010
PhenologyDat.df1 <- PhenologyDat.df[which(PhenologyDat.df$Year>2010 & PhenologyDat.df$Year<2021),]
nrow(PhenologyDat.df1)
# Drop exotic Apocynaceae
Drop <- c("Gomphocarpus", "curassavica")
PhenologyDat.df1 <- PhenologyDat.df1[!grepl(paste(Drop, collapse= "|"), PhenologyDat.df1$Host), ]
nrow(PhenologyDat.df1)
head(PhenologyDat.df1)
unique(PhenologyDat.df1$Host)
# Drop blanks for Host
PhenologyDat.df1 <- PhenologyDat.df1[!(PhenologyDat.df1$Host=="."), ]
nrow(PhenologyDat.df1)
#
# Delete records with comment of "Ornamental"
PhenologyDat.df1 <- PhenologyDat.df1[!grepl("Ornamental", PhenologyDat.df1$Comments), , drop = FALSE] 
head(PhenologyDat.df1)
nrow(PhenologyDat.df1)
#
# Delete records from David Berman's late season study
PhenologyDat.df1 <- PhenologyDat.df1[!grepl("Berman", PhenologyDat.df1$Comments), , drop = FALSE] 
head(PhenologyDat.df1)
nrow(PhenologyDat.df1)
#
# Randomize order of data
set.seed(42)
rows <- sample(nrow(PhenologyDat.df1))
PhenologyDat.df1 <- PhenologyDat.df1[rows, ]
#
##################################################################################
### Thin data for total year
##################################################################################
##
# Thin data by 10 km 
PhenologyDatList <- list()
SpeciesList <- unique(PhenologyDat.df1$Host)
length(SpeciesList)
#
Count <- 0
#
for(SpeciesName in SpeciesList) {
  #SpeciesName <- "A. variegata"
  PhenologyDat.df2 <- PhenologyDat.df1[which(PhenologyDat.df1$Host==SpeciesName),]
  nrow(PhenologyDat.df2)
  # Randomize order of rows
  set.seed(42)
  rows <- sample(nrow(PhenologyDat.df2))
  PhenologyDat.df3 <- PhenologyDat.df2[rows, ]
  # Remove duplicate values for Lat/Lon
  PhenologyDat.df4 <- PhenologyDat.df3[!duplicated(PhenologyDat.df3[c("lat", "lon")]),]
  nrow(PhenologyDat.df4)
  if(nrow(PhenologyDat.df4)< 3) {
    if(nrow(PhenologyDat.df4)==2) {
      # calculate distance between the two points
      PointDist <- distm(c(PhenologyDat.df4$lon[1], PhenologyDat.df4$lat[1]), c(PhenologyDat.df4$lon[2], PhenologyDat.df4$lat[2]), fun = distHaversine)
      # Delete second record if closer than 10 km (10,000 m)
      if(PointDist < 10000) {
         PhenologyDat.df5 <- PhenologyDat.df4[1,]
      } else {
         PhenologyDat.df5 <- PhenologyDat.df4
      }
    } else {
      PhenologyDat.df5 <- PhenologyDat.df4
    }
  } else {
    # Spatially thin to 10 km
    # Select only lat lon and id columns
    pres.df <- PhenologyDat.df4[,c("lat","lon", "ID")]
    head(pres.df)
    nrow(pres.df)
    # Omit any rows with NA values
    pres.df <- na.omit(pres.df)
    pres.df$Species <- Species
    ## Spatially thin presence points to 10 km
    presthin.df <- as.data.frame(thin(pres.df, lat.col="lat", long.col="lon", spec.col = "Species", thin.par=10, reps=1, write.file=FALSE, locs.thinned.list.return = TRUE))
    colnames(presthin.df) <- c("lon", "lat")
    # Merge data with original roadkill by lat/long
    PhenologyDat.df5 <- merge(presthin.df, PhenologyDat.df4, by=c("lon", "lat"))
    head(PhenologyDat.df5)
    nrow(PhenologyDat.df5)
  }
  # Save data to list
  Count <- Count + 1
  PhenologyDatList[[Count]] <- PhenologyDat.df5
}  
##########################
PhenologyDatAll.df <- do.call(rbind, PhenologyDatList)
colnames(PhenologyDatAll.df)
# Re-order columns
PhenologyDatAll.df <- PhenologyDatAll.df[, c(3,4,1,2,5:ncol(PhenologyDatAll.df))]
head(PhenologyDatAll.df)
tail(PhenologyDatAll.df)
nrow(PhenologyDatAll.df)
write.table(PhenologyDatAll.df, file=paste0(Species, "_PhenologyThinned10kmSC2011_20_AllNativeMilkweeds.txt"), sep=",", row.names=FALSE)
##
###########################
## Obtain frequency of plants by species
HostFreq.df <- as.data.frame(table(PhenologyDatAll.df$Host), stringsAsFactors=FALSE)
colnames(HostFreq.df)[1] <- c("Host")
#
write.table(HostFreq.df, file=paste0(Species, "_HostFrequencyforThinned10kmSC2011_20_AllNativeMilkweeds.csv"), sep=",", row.names=FALSE)
##################################################################################
##################################################################################
### Thin data for Spring Vs Fall
##################################################################################
##
WeekNumMinList <- c(11, 32)
WeekNumMaxList <- c(28, 47)
SeasonList <- c("Spring", "Fall")
# Thin data by 10 km 
for(i in 1:2) {
  #i=1
  WeekNumMin <- WeekNumMinList[i]
  WeekNumMax <- WeekNumMaxList[i]
  Season <- SeasonList[i]    
  PhenologyDat.df1a <- PhenologyDat.df1[ which(PhenologyDat.df1$WeekNum>=WeekNumMin & PhenologyDat.df1$WeekNum<=WeekNumMax),]
  nrow(PhenologyDat.df1a)
  PhenologyDat.df1a$WeekNum
  Count <- 0
  PhenologyDatList <- list()
  SpeciesList <- unique(PhenologyDat.df1a$Host)
  length(SpeciesList)
  for(SpeciesName in SpeciesList) {
    #SpeciesName <- "A. variegata"
    PhenologyDat.df2 <- PhenologyDat.df1a[which(PhenologyDat.df1a$Host==SpeciesName),]
    nrow(PhenologyDat.df2)
    # Randomize order of rows
    set.seed(42)
    rows <- sample(nrow(PhenologyDat.df2))
    PhenologyDat.df3 <- PhenologyDat.df2[rows, ]
    # Remove duplicate values for Lat/Lon
    PhenologyDat.df4 <- PhenologyDat.df3[!duplicated(PhenologyDat.df3[c("lat", "lon")]),]
    nrow(PhenologyDat.df4)
    if(nrow(PhenologyDat.df4)< 3) {
      if(nrow(PhenologyDat.df4)==2) {
        # calculate distance between the two points
        PointDist <- distm(c(PhenologyDat.df4$lon[1], PhenologyDat.df4$lat[1]), c(PhenologyDat.df4$lon[2], PhenologyDat.df4$lat[2]), fun = distHaversine)
        # Delete second record if closer than 10 km (10,000 m)
        if(PointDist < 10000) {
           PhenologyDat.df5 <- PhenologyDat.df4[1,]
        } else {
           PhenologyDat.df5 <- PhenologyDat.df4
        }
      } else {
        PhenologyDat.df5 <- PhenologyDat.df4
      }
    } else {
      # Spatially thin to 10 km
      # Select only lat lon and id columns
      pres.df <- PhenologyDat.df4[,c("lat","lon", "ID")]
      head(pres.df)
      nrow(pres.df)
      # Omit any rows with NA values
      pres.df <- na.omit(pres.df)
      pres.df$Species <- Species
      ## Spatially thin presence points to 10 km
      presthin.df <- as.data.frame(thin(pres.df, lat.col="lat", long.col="lon", spec.col = "Species", thin.par=10, reps=1, write.file=FALSE, locs.thinned.list.return = TRUE))
      colnames(presthin.df) <- c("lon", "lat")
      # Merge data with original roadkill by lat/long
      PhenologyDat.df5 <- merge(presthin.df, PhenologyDat.df4, by=c("lon", "lat"))
      head(PhenologyDat.df5)
      nrow(PhenologyDat.df5)
    }
    # Save data to list
    Count <- Count + 1
    PhenologyDatList[[Count]] <- PhenologyDat.df5
  }   
  PhenologyDatAll.df <- do.call(rbind, PhenologyDatList)
  colnames(PhenologyDatAll.df)
  # Re-order columns
  PhenologyDatAll.df <- PhenologyDatAll.df[, c(3,4,1,2,5:ncol(PhenologyDatAll.df))]
  head(PhenologyDatAll.df)
  tail(PhenologyDatAll.df)
  nrow(PhenologyDatAll.df)
  write.table(PhenologyDatAll.df, file=paste0(Species, "_Phenology", Season, "Thinned10kmSC2011_20_AllNativeMilkweeds.csv"), sep=",", row.names=FALSE)
  ## Obtain frequency of plants by species
  HostFreq.df <- as.data.frame(table(PhenologyDatAll.df$Host), stringsAsFactors=FALSE)
  colnames(HostFreq.df)[1] <- c("Host")
  #
  write.table(HostFreq.df, file=paste0(Species, "_HostFrequencyfor", Season, "Thinned10kmSC2011_20_AllNativeMilkweeds.csv"), sep=",", row.names=FALSE)
}


