setwd("~/XXXXX") # edit this to point to wherever you have the unzipped folder located

library(sf)
library(raster)
library(spatstat)
library(ald)
library(foreach)
library(doParallel)
library(stringi) #used to manipulate Sys.time() for appending to output file names

load("helper_objects_V5.2_,20240925_114930.067491")# has all functions and other objects needed for analysis

dev.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="BLMoutput/"

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

#only need to generate once, so put inside loop that won't execute when sourcing the script, change if need to run 
if(0){  #set to 1 if simulating new dens, 0 if using ones already created
  #Simulate a set of simulated dens for use in disturbance function
  library(foreach)
  library(doParallel)
  
    for (e.y in 2025){  
      print(Sys.time())
      writeLines(paste("starting den sim for emergence year",e.y))
      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=e.y)
        
        stopCluster(cl)
        print(Sys.time()-temp.time)
      }
      saveRDS(sim.dens,paste0(out.path,"sim.dens",e.y-1,"-",e.y,"_", insert.time(), ".rds"))
    }#for (e.y)  
  
sim.dens2425 =sim.dens
}#end if(0)


############################
############################
#read in development data
######################################
#read coastline (used in plotting)
ak.coast  = st_read("AK coast/ak_coast63k_aa154n83_Line_Ncoast.shp")
ak.coast = st_zm(ak.coast) #drop 3rd dimension as we don't need it (won't do anything if already 2d)
ak.coast = st_transform(ak.coast,akproj)


dev.layers = st_layers("denning_analysis.gdb") #gives info on features within the GDB
# Available layers:
#          layer_name                 geometry_type features fields              crs_name
# 1            routes 3D Measured Multi Line String        8      4 NAD83 / Alaska Albers
# 2             lakes              3D Multi Polygon       11      5 NAD83 / Alaska Albers
# 3 combined_features              3D Multi Polygon       14      5 NAD83 / Alaska Albers
# 4             wells                      3D Point        2      3 WGS 84 + EGM96 height
# 5 routes_1mile_buff                 Multi Polygon        1      2 NAD83 / Alaska Albers
# 6 wells__1mile_buff                 Multi Polygon        2      7 WGS 84 + EGM96 height


dev.layers.names = dev.layers$name

dev.layers.names[3] #only need combined_features

#use the original layer names and assign the feature info to that name
combined_features = import.shapefile(shape.dir ="", file.name = "denning_analysis.gdb", layer.name = dev.layers.names[3], start.stop.names = c("Start_25","End_25")  )  
       
#check if any have start with no end and vice versa  2FLIR
check.start.stop(combined_features,start.stop.names)  

boundingbox(combined_features[[1]]) #get bounding box, perhaps to use with plots. 
#window: rectangle = [37551.51, 157959.45] x [2227042.7, 2311078.3] units

bbox.x = c(37551.51, 157959.45)
bbox.y = c(2227042.7, 2311078.3)

#pull and compact dev for each emergence year

BLM.2425.dev.2flir = combine.compact.feature(combined_features,col.names = start.stop.names) 


##########################
##########################
#generate FLIR data
##########################

#in ALL years                        #2 FLIR: 12/1-12-25; 12/15-1/10

FLIR.2x.2425 = do.call(union.owin,BLM.2425.dev.2flir$unique.dev.list.buff) #combine all objects doing 2 flir

#plot them if you like
plot(FLIR.2x.2425, main="FLIR 2425")


FLIR.owin.list.2425 = list(FLIR.2x.2425, FLIR.2x.2425)

FLIR.2x.dates.2425 = data.frame(start = c(strptime('2024-12-01','%Y-%m-%d'),strptime('2024-12-25','%Y-%m-%d')),
                                end   = c(strptime('2024-12-15','%Y-%m-%d'),strptime('2025-01-10','%Y-%m-%d')))

FLIR.dates.2425 = FLIR.2x.dates.2425

#create complete compacted dev objects with buffers
BLM.2425.dev    = BLM.2425.dev.2flir

#add FLIR info to dev.objects

BLM.2425.dev    = append(BLM.2425.dev, list(FLIR.owin.list=FLIR.owin.list.2425, FLIR.dates=FLIR.dates.2425, study.area.owin = FLIR.2x.2425 ) ) 

#save all yearly dev objects
saveRDS(BLM.2425.dev, paste0(out.path,"BLM.2425.dev_", insert.time(),  ".rds"))

#################

#plot all the FLIR objects with dates of flir  


pdf(paste0(out.path,"Dev.FLIR.plots_2425_",insert.time(), ".pdf"),width = 10,height = 7.5)
tmp=BLM.2425.dev
plot(tmp[[1]],xlim=bbox.x,ylim=bbox.y,col="black",
     main=paste("BLM 2024-2025","development footprint (unbuffered)"))
plot(ak.coast[,1],col="red",lwd=2,add=T)
for(i in 1:length(tmp$FLIR.owin.list)){
  plot(tmp$FLIR.owin.list[[i]],col=i+1,xlim=bbox.x,ylim=bbox.y,
       main=paste("2024-2025 FLIR date",tmp$FLIR.dates[i,1],"to",tmp$FLIR.dates[i,2])); plot(ak.coast[,1],col="black",lwd=2,add=T)}
dev.off()


##########################
##########################
#we now have all the input data available
# can run disturbance sim
library(foreach)
library(doParallel)

# #if you want to read dev object back:
# BLM.2425.dev = readRDS(paste0(out.path,list.files(path=out.path,pattern = "BLM\\.2425\\.dev.*\\.rds$")))
list2env(BLM.2425.dev,globalenv())  #put all elements of BLM.2425.dev list into global environment
names(unique.dev.dates) = c("start","end") #need to not have the years in the column names

#if you need to read in existing sim.dens
# sim.dens = readRDS(paste0(out.path,"sim.dens2024-2025_20240807_112853.783789.rds"))

#run den sim

  cl <- makePSOCKcluster(15) 
  
  {
    temp.time=Sys.time()
    registerDoParallel(cl)
    x=vector()
    x=foreach(i=1:length(sim.dens),.packages=c("terra","st","spatstat"),.errorhandling='pass',.verbose=F,.multicombine=TRUE) %dopar%
      disturbance.func.parallel.FLIRvar(dens=sim.dens[[i]],
                                        dev.obj          =unique.dev.owin,
                                        dev.obj.list     =unique.dev.list,
                                        dev.obj.list.buff=unique.dev.list.buff, 
                                        dev.obj.dates    =unique.dev.dates,
                                        FLIR.regions     =FLIR.owin.list,
                                        FLIR.dates       =FLIR.dates)
    stopCluster(cl)
    print(Sys.time()-temp.time)
    
    sim.dist=list()
    sim.dist[[1]]=lapply(x,FUN=function(xx) xx[[1]])
    sim.dist[[2]]=do.call(rbind,lapply(x,FUN = function(xx) xx[[2]]))
    sim.dist[[3]]=lapply(x,FUN = function(xx) xx[[3]]) 
    
  }
  saveRDS(sim.dist,paste0(out.path,"sim.dist2425_", insert.time(),".rds"))
  assign("sim.dist2425",sim.dist) #rename a copy to have in the workspace
  env <- foreach:::.foreachGlobals; rm(list=ls(name=env), pos=env); rm(env) #do some cleanup a la https://stackoverflow.com/questions/25097729/un-register-a-doparallel-cluster


#######################
#######################
#post processing and summary

#read in existing saved disturbance objects, presumes only one copy for each emergence year
# uncomment if needed, should be in workspace if above code run
# sim.dist2425    = readRDS(paste0(out.path,list.files(path=out.path,pattern = "sim.dist24.*\\.rds$"       )))


sink(paste0(out.path,"BLM Summary Stat output_", insert.time(),".txt"),split=T)#write the output to a file split sends to new sink and current output stream

writeLines("BLM den disturbance output for denning year 2024-2025")
writeLines(paste0("using den disturbance model ",disturbance.model.version))
tmp = summary.stats(sim.dist2425)
  print(tmp, row.names = FALSE)
  writeLines("\n\n\n")

sink()
