## 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)
# Specify  directory for saved functions
FunctDirect <- "C:/Users/jamesltracy61/Documents/R/win-library/"
setwd(FunctDirect)
#
# Read in User Defined Functions
source("EvalStatVars_Summary_Function2.R")
source("SpatialFilter_Function.R")
source("Round2_Function.R")
#
InDirect <- "K:/Monarch"
setwd(InDirect)
Species <- "Asclepias_capricornu"
# Read in data or estimated roadkill/100 m from 2016 hotspot and nonhotspot roadkill data
PhenologyDat.df <- data.frame(read.csv(paste0(Species, "_PhenologyUnthinnedWeekYearMxtCoreSC2020.csv"), header = TRUE, sep=',', stringsAsFactors=FALSE))
head(PhenologyDat.df)
tail(PhenologyDat.df)
nrow(PhenologyDat.df)
#str(PhenologyDat.df)
# Only keep data after 2010
PhenologyDat.df1 <- PhenologyDat.df[which(PhenologyDat.df$Year>2010),]
head(PhenologyDat.df1)
# Keep only iNat data
PhenologyDat.df1 <- PhenologyDat.df1[which(PhenologyDat.df1$Source=="iNat"),]
# Thin data by 30 km within same week and year
YearList <- unique(PhenologyDat.df1$Year)
WeekList <- unique(PhenologyDat.df1$WeekNum)
PhenologyDatList <- list()
#
Count <- 0
for(Year in YearList) {
  #Year <- 2019
  PhenologyDat.df2 <-  PhenologyDat.df1[which(PhenologyDat.df1$Year==Year),]
  for(Week in WeekList) {
    #Week <- 37
    PhenologyDat.df3 <-  PhenologyDat.df2[which(PhenologyDat.df2$WeekNum==Week),]
    head(PhenologyDat.df3)
    nrow(PhenologyDat.df3)
    if(nrow(PhenologyDat.df3)< 3) {
      PhenologyDat.df4 <- PhenologyDat.df3
    } else {
      # Spatially thin to 30 km
      # Select only lat lon and id columns
      pres.df <- PhenologyDat.df3[,c("lat","lon", "ID")]
      head(pres.df)
      nrow(pres.df)
      # Omit any rows with NA values
      pres.df <- na.omit(pres.df)
      # Convert point data.frame to SpatialPointsDataFrame
      # First, specify xy coordinates
      xy <- pres.df[,c(2,1)]
      # Create spatial points data frame
      # Set projection
      CRS.WGS84 <- CRS("+init=epsg:4326")
      #
      pres.spdf <- SpatialPointsDataFrame(coords=xy, data=pres.df, proj4string=CRS.WGS84)
      #head(pres.spdf)
  #    ## Plot presence data with country background
  #    # Read ESRI shapefile into CountryShp below.
  #    AdminShp <- readOGR(".", paste0("admin00"))
  #    #
  #    # xlim sets the LONGITUDE limits and ylim the LATITUDE limits
  #    dev.new()
  #    plot(pres.spdf, axes=TRUE, col='red')
  #    plot(AdminShp,  add=TRUE, col='light yellow')
  #    plot(pres.spdf, add=TRUE, col='red')
  #    #
      ## Spatially thin presence points to 30 km
      presthin.spdf <- SpatialFilter(pres.spdf, dist=30, mapUnits=F)
  #    plot(presthin.spdf, add=TRUE, col='blue')
      ## Optional, use zm() to zoom into points on plot
      #library(zoom)
      #zm()
      # Convert spatial points data frame back to normal data frame
      presthin.df <- data.frame(presthin.spdf, stringsAsFactors=FALSE)
      head(presthin.df)
      nrow(presthin.df)
      # Keep only names of thinned transects
      presthin.df <- presthin.df[c("ID")]
      #str(presthin.df)
      # Convert presthin.df$Name from factor to character
      presthin.df$ID <- as.character(presthin.df$ID)
      #
      # Merge data with original roadkill by lat/long
      PhenologyDat.df4 <- merge(presthin.df, PhenologyDat.df3, by.x="ID", by.y="ID")
      head(PhenologyDat.df4)
      nrow(PhenologyDat.df4)
    }
    # Save data to list
    Count <- Count + 1
    PhenologyDatList[[Count]] <- PhenologyDat.df4
  }
}
##########################
PhenologyDatAll.df <- do.call(rbind, PhenologyDatList)
head(PhenologyDatAll.df)
tail(PhenologyDatAll.df)
nrow(PhenologyDatAll.df)
write.table(PhenologyDatAll.df, file=paste0(Species, "_PhenologyThinnedWeekYear30kmMxtCoreSC2020.csv"), sep=",", row.names=FALSE)
##
##################################################################################
## Randomly thin to 525 points
PhenologyDatAll.df2 <- PhenologyDatAll.df[sample(nrow(PhenologyDatAll.df), 525), ]
head(PhenologyDatAll.df2)
tail(PhenologyDatAll.df2)
nrow(PhenologyDatAll.df2)
write.table(PhenologyDatAll.df2, file=paste0(Species, "_PhenologyThinned525WeekYear30kmMxtCoreSC2020.csv"), sep=",", row.names=FALSE)


