library(rgdal) #some of these packages are no longer in general distribution.  Search CRAN archives for most recent versions
library(rgeos)
library(raster)
library(spatstat)
library(ald)
library(stringi) #used to manipulate Sys.time() for appending to output file names

load("...[edit path for your location of this file]/helper_objects19oct2023.RData")

setwd("~/...[edit path]") #common working directory, set as neeeded
dev.path="...[edit path]"       # this is the common directory for where development objects are stored (or their sub directories) can also hard code below if needed
out.path="OutputFiles/" #This is a folder (created as necessary) to contain output files to help keep files more organized
if (! dir.exists(out.path) ) dir.create(out.path) #create an output folder if doesn't exist


filename.info = "2024BLM"
summ.filename = paste0(out.path,"run_output_"      ,filename.info,insert.time(), ".txt")
 pdf.filename = paste0(out.path,"run_output_plots_",filename.info,insert.time(), ".pdf")
sink(summ.filename,split=T)#write the output to a file split sends to new sink and current output stream
pdf(pdf.filename,width = 11,height = 8.5)

 for (scenario in 0){ #loop for each scenario (0 is ALL)
writeLines(paste("STARTING SCENARIO",ifelse(scenario==0,"ALL",scenario)))
  
writeLines("starting function")
print(Sys.time())

redoDev=T    #T forces reimport of dev info from shape files (used if dev gis files have been modified)
redoDist=T   #forces redo of disturbance sim (used if dev info has been changed and model need rerunning despite existence of previous run)

###########
#create or import dev info
###########
writeLines("starting dev.info")
print(Sys.time())


tmp.files.dev=list.files(path=out.path,pattern = "dev.info.*\\.rds$")
if(length(tmp.files.dev)>0 & redoDev != T){ #dev.infoXXXXXXXX_XXXXXX.rds exists, so read in
  dev.file.used <<- paste0(out.path,tmp.files.dev[length(tmp.files.dev)])
  writeLines(paste("reading in dev.file",dev.file.used))
  dev.info = readRDS(dev.file.used) #always choose latest one by default since latest date/time
  list2env(dev.info,globalenv())  #put all elements of list into global environment
}else { 
  #else create the info.  can also be run manually if dev info has been updated. 
  {
    ##########
    #Roads
    ##########
    roads=import.shapefile(shape.dir=dev.path,
                           #the following file is slighltly modified to be able to work with rgdal:
                           #had to split it through the "holes" due to the error in reading file
                           file.name='Routes_IHA_Analysis_6_20_2023_buffer.shp', #added the desired mile buffer, will add extra mile internal to disturbance function as well
                                                                                 #also had to split it through the "holes" due to the error in reading file 
                           study.area.range=c(-442958,545262,2151448,2433230),#x1,x2,y1,y2
                           type='polygon',
                           proj="+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")
    
    
    roads_dates=data.frame(start23=rep(strptime("2024-02-01",'%Y-%m-%d'),length(roads[[2]])),end24=rep(strptime("2024-05-15",'%Y-%m-%d'),length(roads[[2]])))
      
   which2023 = which(!is.na(roads_dates$start23)) #which features actually operational in 2023-4
   
    roads = combine.feature(roads,which2023 )
    roads_dates = roads_dates[which2023,]
    names(roads_dates)=c("start","end")

    
    

  ##########
  #Lakes  #this is reading two separate sets of lakes, with different dates
  ##########
  lakes_Tulageak=import.shapefile(shape.dir=dev.path,
                         file.name='Water_Source_Lakes_Tulageak.shp',
                         study.area.range=c(-442958,545262,2151448,2433230),#x1,x2,y1,y2
                         type='polygon',
                         proj="+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")
  
  lakes_Tulageak_dates=data.frame(start23=rep(strptime("2024-02-01",'%Y-%m-%d'),length(lakes_Tulageak[[2]])),end24=rep(strptime("2024-04-30",'%Y-%m-%d'),length(lakes_Tulageak[[2]])))
  
  
  which2023 = which(!is.na(lakes_Tulageak_dates$start23)) #which features actually operational in 2023-4
  
  lakes_Tulageak = combine.feature(lakes_Tulageak,which2023 )
  lakes_Tulageak_dates = lakes_Tulageak_dates[which2023,]
  names(lakes_Tulageak_dates)=c("start","end")
  
  
  lakes_ESimpson=import.shapefile(shape.dir=dev.path,
                                  file.name='Water_Source_Lakes_ESimpson.shp',
                                  study.area.range=c(-442958,545262,2151448,2433230),#x1,x2,y1,y2
                                  type='polygon',
                                  proj="+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")
  
  lakes_ESimpson_dates=data.frame(start23=rep(strptime("2024-03-01",'%Y-%m-%d'),length(lakes_ESimpson[[2]])),end24=rep(strptime("2024-03-31",'%Y-%m-%d'),length(lakes_ESimpson[[2]])))
  
  
  which2023 = which(!is.na(lakes_ESimpson_dates$start23)) #which features actually operational in 2023-4
  
  lakes_ESimpson = combine.feature(lakes_ESimpson,which2023 )
  lakes_ESimpson_dates = lakes_ESimpson_dates[which2023,]
  names(lakes_ESimpson_dates)=c("start","end")
    
  
    #########################################
    #legacy wells-  camps are within 1 mile radius of wells, but need to buffer camp too.
    #creating mile buffer around well and then including as dev object will subsequently buffer a second mile so camps can be anywhere within a mile of the well point
    #
    
    wells =  spTransform(readOGR("~/MMM/projects/2024 BLM IHA legacy well program/2024 BLM legacy well GIS/Legacy_Wells.shp"),akproj)
    
    obj.list=list()
    for(i in 1:length(wells)){
      #tmp.xy=wells@coords[i,]
      tmp.owin=owin(xrange=c(-442958,545262),yrange=c(2151448,2433230))
      #tmp.point=ppp(tmp.xy[1],tmp.xy[2],window=tmp.owin)
      tmp.point=ppp(wells@coords[i,1],wells@coords[i,2],window=tmp.owin)
      tmp.owin=as.owin(bufftess(tmp.point,c(0,1610)))#1610 buffer added to point features to convert them to polygons
      obj.list[[i]]=tmp.owin 
    }
    
    #Combine into a object 
    obj=obj.list[[1]]
    if(length(obj.list)>1){
      for(i in 2:length(obj.list)){
        obj=union.owin(obj,obj.list[[i]])
      }
    }
  wells=list(obj,obj.list)
  wells_dates=data.frame(start23=rep(strptime("2024-02-01",'%Y-%m-%d'),length(wells[[2]])),end24=rep(strptime("2024-05-15",'%Y-%m-%d'),length(wells[[2]])))
  names(wells_dates)=c("start","end")
  
  }
  
    
  
    #########################
    #FLIR
    #########################
    
    
  flir1=union.owin(roads[[1]],lakes_Tulageak[[1]],lakes_ESimpson[[1]],wells[[1]])
  
    
    flir.owin.list=list()
    flir.owin.list[[1]]=flir1
    flir.owin.list[[2]]=flir1
      
    flir.dates=data.frame(start=c(strptime('2023-12-01','%Y-%m-%d'),strptime('2023-12-25','%Y-%m-%d')),
                          end  =c(strptime('2023-12-15','%Y-%m-%d'),strptime('2024-01-10','%Y-%m-%d')))
    #using same dates as AOGA, will likely be right after AOGA surveys but slightly earlier is more conservative

  
  
  ###########################################
  #Combine all owins
  
    #dev.owin=union.owin(roads[[1]],lakes[[1]],wells[[1]])
    dev.owin=union.owin(roads[[1]],lakes_Tulageak[[1]],lakes_ESimpson[[1]],wells[[1]])
    
  {
    #Develop a list of owin objects

    dev.list=append(roads[[2]],lakes_Tulageak[[2]])
    dev.list=append(dev.list,lakes_ESimpson[[2]])
    dev.list=append(dev.list,wells[[2]])
    
    
    #Develop date/time data.frame for owin objects
    dev.dates=rbind(roads_dates,lakes_Tulageak_dates,lakes_ESimpson_dates,wells_dates)
    
    print(paste("dev list length",length(dev.list)))
    print(paste("dev date length",dim(dev.dates)[1]))
    
    tmp = compact.dev(dev.dates,dev.list) #create compacted dev list and dates (and buffered version too) unique.dev.dates  unique.dev.list   unique.dev.list.buff
    list2env(tmp,globalenv())  #put all elements of returned list into global environment
    
    #get extent for the project area
    study.area.owin.old=study.area.owin  #imported study area is huge region, we can reduce it a bunch keeping only for reference later
    study.area.owin = expand.owin(dev.owin, distance=disturbance.distance)  #making the study area a buffered version of combined dev objects
    
    #exporting all dev and FLIR info 
    dev.info=list(unique.dev.dates    =unique.dev.dates,
                  unique.dev.list     =unique.dev.list,
                  unique.dev.list.buff=unique.dev.list.buff,
                  dev.owin            =dev.owin,
                  dev.list            =dev.list,
                  dev.dates           =dev.dates,
                  flir.owin.list      =flir.owin.list,
                  flir.dates          =flir.dates,
                  study.area.owin     =study.area.owin)
    # saveRDS(dev.info,"ACoutput/dev.info.rds") 
    dev.file.used <<- paste0(out.path,"dev.info", insert.time(),  ".rds")
    saveRDS(dev.info,dev.file.used )
  }
}#end import dev chunk
 
plot(dev.owin,main=paste("Scenario",ifelse(scenario==0,"ALL",scenario)))

######################################
######################################
#Simulate dens

writeLines("starting den sim")
print(Sys.time())


#Simulate a set of simulated dens for use in disturbance function
library(foreach)
library(doParallel)

tmp.files.dens=list.files(path=out.path,pattern = "sim.dens2024_.*\\.rds$") #file name format should be sim.densYYYY_[datecreatedinfo].rds  where YYYYis emergence year

if (length(tmp.files.dens)>0){     #read it back if already created
  den.file.used <<- paste0(out.path,tmp.files.dens[length(tmp.files.dens)])
  writeLines(paste("reading in sim.dens",den.file.used))
  
  sim.dens = readRDS(den.file.used)
}else {                                       #otherwise do den sim
  print("no existing sim.dens file, running den sim")
  print(Sys.time())
  
  n.sim =10000
  {
    temp.time=Sys.time()
    cl <- makePSOCKcluster(detectCores()-1)#  15)
    registerDoParallel(cl)
    sim.dens=vector()
    sim.dens=foreach(i=1:n.sim,.packages=c("raster","maptools","spatstat","ald"),.errorhandling='pass',.verbose=F,.multicombine=TRUE) %dopar%
      den.func.parallel(emerg.year=2024)
    
    stopCluster(cl)
    print(Sys.time()-temp.time)
  }
  den.file.used <<- paste0(out.path,"sim.dens2024_", insert.time(), ".rds")
  saveRDS(sim.dens,den.file.used)
}


#disturbance analysis

writeLines("starting disturbance section")
print(Sys.time())


tmp.files.dist=list.files(path=out.path,pattern = "sim.dist2024_.*\\.rds$") #file name format should be sim.distYYYY_[datecreatedinfo].rds  where YYYYis emergence year
if(length(tmp.files.dist)>0 & redoDist !=T){ #dev.infoXXXXXXXX_XXXXXX.rds exists, so read in
  dist.file.used <<-paste0(out.path,tmp.files.dist[length(tmp.files.dist)])
  writeLines(paste("reading in sim.dist file",dist.file.used))
  sim.dist = readRDS(dist.file.used) #always choose latest one by default since latest date/time
}else {
  #Run den disturbance model
  print("no existing sim.dist file or forced rerun, running disturbance model sim")
  
  library(foreach)
  library(doParallel)
  
  cl <- makePSOCKcluster(15)
  temp.time=Sys.time()
  
  registerDoParallel(cl)
  x=vector()
  x=foreach(i=1:length(sim.dens),.packages=c("raster","maptools","spatstat"),.errorhandling='pass',.verbose=F,.multicombine=TRUE) %dopar% 
    disturbance.func.parallel.FLIRvar(dens=sim.dens[[i]],
                                      activity.region=study.area.owin,
                                      dev.obj=dev.owin,
                                      dev.obj.list=unique.dev.list,
                                      dev.obj.list.buff=unique.dev.list.buff, ####adding in buffered set of owins for each unique date
                                      dev.obj.dates=unique.dev.dates,
                                      FLIR=T,#Include FLIR T/F
                                      FLIR.regions=flir.owin.list,
                                      FLIR.dates=flir.dates,
                                      min.days.between.flir=1,
                                      FLIR.detec.prob=data.frame(x=flir3),
                                      litter.surv.coef=surv.param,
                                      disturbance.distance=1610)
  stopCluster(cl)
  print(Sys.time()-temp.time)
  
  
  
  den.results=list()
  for(i in 1:length(x)){
    den.results[[i]]=x[[i]][[1]][[1]]
    if(i==1) results.df=x[[i]][[2]] else results.df=rbind(results.df,x[[i]][[2]])
    if(i==1) flir.dates.df=x[[i]][[3]] else flir.dates.df=rbind(flir.dates.df,x[[i]][[3]])
  }
  
  sim.dist=list()
  sim.dist[[1]]=den.results
  sim.dist[[2]]=results.df
  sim.dist[[3]]=flir.dates.df
  
  dist.file.used <<- paste0(out.path,"sim.dist2024_", insert.time(),".rds")
  saveRDS(sim.dist,dist.file.used)
  # saveRDS(sim.dist,"ACoutput/simdist2024AOGAalone.rds") 
}#end run dist function


# sim.dist = readRDS("ACoutput/simdist2024AOGAalone.rds")#read it back 

summ.stats=summary.stats(sim.dist,LevANonSerious.names=c('nonSeriousLevA.ncubs'),LevALethal.names=c('Lethal.early.ncubs','Lethal.late.ncubs','Lethal.post.ncubs','SeriousLevA.ncubs'))
writeLines(filename.info)
writeLines(paste("development info file:\n",dev.file.used))
writeLines(paste("den sim file:\n",den.file.used))
writeLines(paste("disturbance sim file:\n",dist.file.used))
print(summ.stats)

writeLines("ending function")
print(Sys.time())

writeLines(paste("ENDING SCENARIO",ifelse(scenario==0,"ALL",scenario)))
}#for i
sink()
dev.off()

