### Initialization
cat("Initializing..\n")
### Removes everthing from environment
rm(list=ls())

#Set working directory to somewhere on your local machine
working_dir <- "~/MakeRasters"
setwd(working_dir)
setwd("C:/Users/dick smith/Desktop/2018/Directory/Raster Layers/Real/BIO1")
getwd()
install.packages("rgdal")
# Set user params
# MODIFY TO YOUR RESOLUTION (e.g. and the units of your raster - km? lats and longs?
res=0.008333333

# Change if you are running on a server
numCPUs=1

# Set values - the parameters that describe the resistance relationships (alpha and gamma). Modify these values according to your hypotheses and after looking at how the transformations affect your dataset (by plotting the rasters).

#Test values - makes 3 rasters: 100,0.01; 100,1; 100,100)
alphaVals=c(0,5,10,100,1000)
gammaVals=c(0.1,0.2,0.5,1,2,5,10)

#All combos - makes 5x7 rasters  = 35
#alphaVals=c(0,5,10,100,1000)
#gammaVals=c(0.01,0.1,0.5,1,5,10,100)

# Load in deps
#library(snow) - this is for parallel computing if you have a server
library(sp)
library(rgdal)
library(raster)
library(DT)
library(dbplyr)

#outputs data to current temp directory
rasterOptions(tmpdir="./")

##### Highlight and run from here onwards to the end in one chunk (after modifying raster data and names etc.)

##### Define functions - the main thing
get.multires.tc.lm <- function(curr_args) {
  
  # init
  library(sp)
  library(rgdal)
  library(raster)
  rasterOptions(tmpdir="~/MakeRasters")
  # set alpha and gamma
  curr_alpha=curr_args[["alpha"]]
  curr_gamma=curr_args[["gamma"]]
  
  i=curr_args[["num"]]
  #compile layers
  currOutput=Output
  Output[Pos.tc]
  
  #Consider here if you want to reverse values of the raster, otherwise apply equation directly:
  #this is where you deal with zero and negative values too. Add code as required.Check NA flag in file and convert values if needed to -9999
  
  currOutput[Pos.tc] <- 1 + curr_alpha * (((TC[Pos.tc])-1)/6.6) ^ curr_gamma
  
  #to format NA values
  currOutput[NA.tc] <- 1 + curr_alpha * (((TC[NA.tc])-1)/6.6) ^ curr_gamma
  
  
  # resampleRaster
  resRast <- resample(currOutput, tempRast)
  
  ###Write raster to file. 
  writeRaster(resRast,paste("BIO",i,".asc",sep=""),overwrite=TRUE, NAflag=-9999)
  return(1)
}

#### Preliminary processing
cat("Preliminary processing..\n")

## Load in data - modify path to your formatted raster data .asc
TC <- raster("invert.asc")
Combs.TC <- as.matrix(expand.grid(alpha=alphaVals, gamma=gammaVals))

combosList=list()
for (currRow in 1:nrow(Combs.TC)) {
  combosList[[currRow]] = c(num=currRow,alpha=Combs.TC[[currRow,"alpha"]], gamma=Combs.TC[[currRow,"gamma"]])
}


###Convert to Dataframe and write combinations file of parameters applied to each raster - this allows you to identify what parameters correspond with which generated raster file based on the file name ‘multi_x’.
write.table(data.frame(name=paste("BIO",1:nrow(Combs.TC),".asc", sep=""),Combs.TC),"./BIO_resistances.csv",row.names=FALSE,sep=",")

# watch directory above - may need to create a ‘raster’ folder

#Extract values 
#create output raster
Output <- TC

# Get positive TC and NAs in TC layer 
Pos.tc <- Which(!is.na(TC), cells=TRUE)
NA.tc <- Which (is.na(TC), cells = TRUE)

# make template raster based on your data layer
tempRast = raster(TC)

###### Main Processing
cat("Main processing..\n")

# generate rasters
result=lapply(combosList, FUN =get.multires.tc.lm)

curr_alpha=combosList[1]["alpha"]
curr_gamma=combosList[["gamma"]]

#stop

cat("DONE!\n")
warnings()

