library(lme4)
library(lmerTest)
library(effects)
library(Hmisc)
library(cAIC4)
library(car)
library(MuMIn)
library(MASS)
library(party)
library(gridExtra)
library(latticeExtra)
library(dismo)
library(rgdal)
library(raster)
library(adehabitatHR)
library(rgeos)
library(dplyr)
library(RColorBrewer)
library(viridis)

########################
##### load rasters #####

setwd(".../DATA")
# load AOI, extent and rasters from 1) raster processing.R

######################
##### load model #####

### use model Mi4c from 2) model building.R or run code:
{dat=read.csv("rsf_table.csv", sep=";")
  dat$pb = factor(dat$pb)
  dat$anid = factor(dat$anid)}

{ dat3 = subset(dat, select=-c(d, a, nbiot, sibiot, cvbiot, nveg, cvveg, dbuh, dmun, da)) # g,
  dat3$dw = dat3$dw/1000
  dat3$dpuo = dat3$dpuo/1000
  dat3$dd = dat3$dd/1000
  dat3$dg = dat3$dg/1000}

# Mi4c: Mi4b -hc:dd2 -hc:buh2 -r:mun2 -r:dd2 -r:buh12 -hc:mun2 -siveg2
Mi4c = glmer(pb ~ dpuo + I(dpuo^2) + dw + I(dw^2) + dg + I(dg^2) + hc + r + I(r^2) + 
               mun + I(mun^2) + buh + I(buh^2) + dd + I(dd^2) + 
               DG + siveg + 
               hc : (mun + buh + dd + DG + siveg) + 
               r  : (mun + dd + DG) + 
               (1|anid), 
             data = dat3, family=binomial)
summary(Mi4c)   # AIC=8265.7
summary(Mi4c)$coef


######################
##### projection #####

### raster distances to km
dw   = dw/1000
dpuo = dpuo/1000
dd   = dd/1000
dg   = dg/1000

tmp = summary(Mi4c)$coefficients[,1]
for(i in 1:length(tmp)){
  if(i==1) Mrast = noquote(paste(tmp[[i]])) else
    Mrast = noquote(paste(Mrast, "+", tmp[[i]], "*", sub('I', '', sub(':', '*', sub('I', '', names(tmp)[[i]])))))
}
logit = exp(formula(Mrast))

plot(logit / (1 + logit))
plot((logit / (1 + logit)) * abs(w - 1)) # prevent holes in raster because of NA

## plot projection result
hvk = shapefile("Hauptvorkommen.shp")
hvk = hvk[16:19,]
hvk = spTransform(hvk, CRS(proj4string(logit)))
png("../Mi4c_HSI_proj.png", 
    width=800, height=1200, units="px", pointsize=20)
cuts = seq(0, 1, 0.1)
pal = magma(length(cuts), direction=-1)
plot(logit / (1 + logit) * abs(w-1), axes=F, legend=F, breaks=cuts, col=pal)
plot(logit / (1 + logit) * abs(w-1), legend.only=TRUE, col=pal, breaks=cuts, 
     legend.width=1.5, legend.shrink=0.3,
     axis.args=list(cex.axis=.8), 
     legend.args=list(text='Habitat suitability', side=4, font=2, line=2.5, cex=.8))
plot(hvk, add=T)
dev.off()


#######################################
##### extract area of suitability #####

link.log = logit / (1 + logit) * abs(w-1)
tmp = raster("rast_hc.grd") * AOI
arlog = data.frame(val=link.log@data@values, hc=tmp@data@values)
arlog = na.omit(arlog)
arlog$cuts = 0
cuts = seq(0, 1, 0.25)
arlog$cuts[arlog$val>=cuts[1] & arlog$val<cuts[2]] = paste0("<", cuts[2])
arlog$cuts[arlog$val>=cuts[2] & arlog$val<cuts[3]] = paste0("<", cuts[3])
arlog$cuts[arlog$val>=cuts[3] & arlog$val<cuts[4]] = paste0("<", cuts[4])
arlog$cuts[arlog$val>=cuts[4] & arlog$val<cuts[5]] = paste0("<", cuts[5])
arlog$cuts[arlog$val>=cuts[5] & arlog$val<cuts[6]] = paste0("<", cuts[6])
arlog$cuts[arlog$val>=cuts[6] & arlog$val<cuts[7]] = paste0("<", cuts[7])
arlog$cuts[arlog$val>=cuts[7] & arlog$val<cuts[8]] = paste0("<", cuts[8])
arlog$cuts[arlog$val>=cuts[8] & arlog$val<cuts[9]] = paste0("<", cuts[9])
arlog$cuts[arlog$val>=cuts[9] & arlog$val<cuts[10]] = paste0("<", cuts[10])
arlog$cuts[arlog$val>=cuts[10] & arlog$val<cuts[11]] = paste0("<", cuts[11])

# Raster cell size: 10mx10m = 100mē -> number cells * 100 / 1.000.000 = Area [kmē]
by(arlog$cuts, arlog$cuts, length)*100/1000000                    # area by suitability [kmē]
round(by(arlog$cuts, arlog$cuts, length)*100/nrow(arlog), 2)      # area percentage by suitability [%]