#################################
rm(list = ls())
require("raster")
require("rgdal")
setwd("P:/Hihi/Manuscript_data_objects/Hihi map")

load("final_model.rda")
load("projected_raster.rda")


Var<- read.csv("Variables_original.csv")


DfW <-Var$DfW
sd.DfW <- sd(DfW)
mean.DfW <- mean(DfW)
DfW <-  (DfW - mean.DfW) / sd.DfW
DfW <- as.matrix(DfW)


DfR <-Var$DfR
sd.DfR <- sd(DfR)
mean.DfR <- mean(DfR)
DfR <-  (DfR - mean.DfR) / sd.DfR
DfR <- as.matrix(DfR)

TPI <- Var$TPI
sd.TPI <- sd(TPI)
mean.TPI <- mean(TPI)
TPI <-(TPI - mean.TPI) / sd.TPI
TPI <- as.matrix(TPI)

Aspect<- as.matrix(Var$Aspect)

#create empty raster
blankr <- raster(ncol=3330, nrow=2102, xmn=2630229.89, xmx=2632569.23, ymn=6192974.59, ymx=6195149.12, crs="+proj=nzmg +lat_0=-41 +lon_0=173 +x_0=2510000 +y_0=6023150 +ellps=intl +towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993 +units=m +no_defs")
res(blankr) <- 15
blankr[]<-0


#function to rasterize .shp


shp2raster <- function(shp, mask.raster, label, value, transform = FALSE, proj.from = NA,
                       proj.to = NA, map = TRUE) {
  require(raster, rgdal)
  
  # use transform==TRUE if the polygon is not in the same coordinate system as
  # the output raster, setting proj.from & proj.to to the appropriate
  # projections
  if (transform == TRUE) {
    proj4string(shp) <- proj.from
    shp <- spTransform(shp, proj.to)
  }
  
  # convert the shapefile to a raster based on a standardised background
  # raster
  r <- rasterize(shp, mask.raster)
  # set the cells associated with the shapfile to the specified value
  r[!is.na(r)] <- value
  # merge the new raster with the mask raster and export to the working
  # directory as a tif file
  r <- mask(merge(r, mask.raster), mask.raster, filename = label, format = "GTiff",
            overwrite = T)
  
  # plot map of new raster
  if (map == TRUE) {
    plot(r, main = label, axes = F, box = F)
  }
  
  names(r) <- label
  return(r)
}


#get data on watercourses
watercourses <- readOGR(dsn = ".", layer = "Watercourses")


watercoursesraster <- shp2raster(shp = watercourses,
                                 mask.raster = blankr, label = "Distance from Watercourse", transform = FALSE, value = 1)


NAvalue(watercoursesraster) <- 0


DfWraster<-distance(watercoursesraster)
#plot(DfWraster)



#TPI
TopPI <- raster("15m_TPI_clipped50.tif")
TPIraster = projectRaster(TopPI, DfWraster)
#plot(TPIraster)

#Aspect
Aspectmap<- raster("Aspect.tif")
Aspectraster = projectRaster(Aspectmap, DfWraster)
#plot(Aspectraster)

#Normalise
DfW.s <- (DfWraster-mean.DfW)/sd.DfW
TPI.s <- (TPIraster-mean.TPI)/sd.TPI

TPIRoto<-as.vector(TPI.s)
DfWRoto<-as.vector(DfW.s)
AR<-as.vector(Aspectraster)


TPIRoto[ is.na(TPIRoto) ] <- 0


#break Aspect in to 4 categories
ARnew <- cut(AR, breaks = seq(0, 360, by = 90), labels=c("NtoE", "EtoS", "StoW", "WtoN"))
ARnew[ is.na(ARnew)] <- "NtoE"

#create a dataframe of covariates
newscovs <- data.frame(DfW=DfWRoto,
                       TPI=TPIRoto,
                       Aspect=ARnew
)


DATEnew <- paste0("T", rep(1:4, each = 16))

Datemat<-matrix(DATEnew, nrow=22620, ncol=64, byrow=T)

newyscovs <- list(date=Datemat)

#project <- function(mod, SiteCovs, YearlySiteCovs){
  nssn <- mod@data@numPrimary
  print(paste("DOING SEASON", 1, "OF", nssn))
  psi <- matrix(nrow=nrow(SiteCovs), ncol=nssn)
  psi[,1] <- predict(mod, "psi", SiteCovs)$Predicted
  for(ssn in 2:nssn){
    print(paste("DOING SEASON", ssn, "OF", nssn))
    ysc <- data.frame(lapply(YearlySiteCovs, function(x) x[,ssn-1]))
    for(i in 1:ncol(ysc)) if(!is.numeric(ysc[,i])) {
      ysc[,i] <- factor(ysc[,i], levels=levels(mod@data@yearlySiteCovs[,names(ysc)[i]]))
    nd <- cbind(SiteCovs, ysc)
    col <- predict(mod, "col", nd)$Predicted
    ext <- predict(mod, "ext", nd)$Predicted
    psi[,ssn] <- psi[,ssn-1]*(1-ext) + (1-psi[,ssn-1])*col
  }
  psi
}


#res <- project(finalmodel, newscovs, newyscovs)

d2mvector<-res[,64]

d2mmatrix<-as.matrix(res, ncol=3330, nrow=2102)


finaldata<- raster(ncol=3330, nrow=2102, 
                   xmn=2630229.89, xmx=2632569.23, ymn=6192974.59, ymx=6195149.12, 
                   crs="+proj=nzmg +lat_0=-41 +lon_0=173 +x_0=2510000 +y_0=6023150 +ellps=intl +towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993 +units=m +no_defs")


res(finaldata) <- 15
finaldata[]<-c(d2mvector)
op <- par(mar=rep(0, 4))
plot(finaldata, asp=NA)
par(op)

#install.packages("prettymapr")
library(prettymapr)



writeRaster(finaldata, "finaldata", format = "GTiff", overwrite=TRUE)


lake<- readOGR(dsn = ".", layer = "Lake2")
reserve<- readOGR(dsn = ".", layer = "Reserve")

lakehole <- mask(finaldata, lake, inverse = TRUE)
plot(lakehole)



#install.packages("colorspace")
library(colorspace)

#pal<-choose_palette()
#colourscheme<-pal(24)

load("colourscheme.rda")


HihiPref <- mask(lakehole, reserve, inverse = FALSE, axes=FALSE)
op <- par(mar=rep(1, 4))
plot(HihiPref, axes=FALSE, box=FALSE, asp=1, col=rev(colourscheme))
title("19th May 2017", adj=1, line=-16.5)
addnortharrow(pos = "topleft", padin = c(0.01, 0.001), scale = 0.75,
              lwd = 1, border = "black", cols = c("white", "black"),
              text.col = "black")
addscalebar(plotunit = NULL, plotepsg = 2000, widthhint = 0.25,
            unitcategory = "metric", htin = 0.1, padin = c(0.01, 0.01),
            style = "bar", bar.cols = c("black", "white"), lwd = 1,
            linecol = "black", tick.cex = 0.7, labelpadin = 0.08, label.cex = 0.8,
            label.col = "black", pos = "bottomleft")
par(op)


#writeRaster(HihiPref, "Hihi preference map", format = "GTiff", overwrite=TRUE)