setwd(".../DATA")
library(dismo)
library(rgdal)
library(raster)
library(adehabitatHR)
library(rgeos)
library(dplyr)
library(lubridate)

#################
## raster data ##
#################
# mask
{'AOI <- shapefile("habitat_nsg.shp")
AOI
e <- extent(AOI)
e
mask <- raster(e, crs=AOI, resolution=10, vals=1)
mask
plot(mask)
AOI <- rasterize(AOI, mask, field=1)
plot(AOI, add=T)
plot(AOI)
writeRaster(AOI,"AOI.grd",overwrite=T)'}

AOI <- raster("AOI.grd")
e <- extent(AOI)
mask <- raster(e, crs=AOI, resolution=10, vals=1)

# Shapefile to raster biotope types
{'
{layers <- shapefile("layer_creation.shp")
layers$DG_S_NAHRU[layers$nveg==0]=NA
layers$siveg[layers$nveg==0]=NA
layers$nveg[layers$nveg==0]=NA

names(layers)
head(layers)
cor.test(layers$DG_S_NAHRU,layers$QK_S_NAHRU)}

{w <- rasterize(layers, mask, field=layers$W)
buh <- rasterize(layers, mask, field=layers$BuH)
suf <- rasterize(layers, mask, field=layers$SuF)
mun <- rasterize(layers, mask, field=layers$MuN)
d <- rasterize(layers, mask, field=layers$D)
hc <- rasterize(layers, mask, field=layers$HC)
r <- rasterize(layers, mask, field=layers$R)
g <- rasterize(layers, mask, field=layers$G)
a <- rasterize(layers, mask, field=layers$A)
u <- rasterize(layers, mask, field=layers$U)
puo <- rasterize(layers, mask, field=layers$PuO)} # rasterize biotope types

{DG <- rasterize(layers, mask, field=layers$DG_S_NAHRU)
nveg <- rasterize(layers, mask, field=layers$nveg)
siveg <- rasterize(layers, mask, field=layers$siveg)} # rasterize further variables

{writeRaster(w,"rast_w.grd",overwrite=T)
writeRaster(buh,"rast_buh.grd",overwrite=T)
writeRaster(suf,"rast_suf.grd",overwrite=T)
writeRaster(mun,"rast_mun.grd",overwrite=T)
writeRaster(d,"rast_d.grd",overwrite=T)
writeRaster(hc,"rast_hc.grd",overwrite=T)
writeRaster(r,"rast_r.grd",overwrite=T)
writeRaster(g,"rast_g.grd",overwrite=T)
writeRaster(a,"rast_a.grd",overwrite=T)
writeRaster(u,"rast_u.grd",overwrite=T)
writeRaster(puo,"rast_puo.grd",overwrite=T)
writeRaster(HE,"rast_he.grd",overwrite=T)
writeRaster(nbiot,"rast_nbiot.grd",overwrite=T)
writeRaster(cvbiot,"rast_cvbiot.grd",overwrite=T)
writeRaster(sibiot,"rast_sibiot.grd",overwrite=T)
writeRaster(DG,"rast_dg.grd",overwrite=T)
writeRaster(QK,"rast_qk.grd",overwrite=T)
writeRaster(nveg,"rast_nveg.grd",overwrite=T)
writeRaster(siveg,"rast_siveg.grd",overwrite=T)
writeRaster(cvveg,"rast_cvveg.grd",overwrite=T)} #writeRaster
'}

##### load Rasters
{w=raster("rast_w.grd")
buh=raster("rast_buh.grd")
suf=raster("rast_suf.grd")
mun=raster("rast_mun.grd")
d=raster("rast_d.grd")
hc=raster("rast_hc.grd")
r=raster("rast_r.grd")
g=raster("rast_g.grd")
a=raster("rast_a.grd")
u=raster("rast_u.grd")
puo=raster("rast_puo.grd")
DG=raster("rast_dg.grd")
nveg=raster("rast_nveg.grd")
siveg=raster("rast_siveg.grd")}

# Rasterstack
temp=stack(w,buh,suf,mun,d,hc,r,g,a,u,puo)
temp
names(temp)=c("w","buh","suf","mun","d","hc","r","g","a","u","puo")
# define moving Window
rast=raster(ncol=21,nrow=21,xmn=0)
gf <- focalWeight(rast, 50, "circle")     # d=100
gf=ifelse(gf>0,1,gf)

##### moving Window Raster
w=focal(x=temp[["w"]],w=gf,fun=mean,pad=T,padValue=NA,na.rm=T)
buh=focal(x=temp[["buh"]],w=gf,fun=mean,pad=T,padValue=NA,na.rm=T)
suf=focal(x=temp[["suf"]],w=gf,fun=mean,pad=T,padValue=NA,na.rm=T)
mun=focal(x=temp[["mun"]],w=gf,fun=mean,pad=T,padValue=NA,na.rm=T)
d=focal(x=temp[["d"]],w=gf,fun=mean,pad=T,padValue=NA,na.rm=T)
hc=focal(x=temp[["hc"]],w=gf,fun=mean,pad=T,padValue=NA,na.rm=T)
r=focal(x=temp[["r"]],w=gf,fun=mean,pad=T,padValue=NA,na.rm=T)
g=focal(x=temp[["g"]],w=gf,fun=mean,pad=T,padValue=NA,na.rm=T)
a=focal(x=temp[["a"]],w=gf,fun=mean,pad=T,padValue=NA,na.rm=T)
u=focal(x=temp[["u"]],w=gf,fun=mean,pad=T,padValue=NA,na.rm=T)
puo=focal(x=temp[["puo"]],w=gf,fun=mean,pad=T,padValue=NA,na.rm=T)

##### distance Raster

dw=raster("dist_w.tif")
dw=dw*AOI
dbuh=raster("dist_buh.tif")
dbuh=dbuh*AOI
dsuf=raster("dist_suf.tif")
dsuf=dsuf*AOI
dmun=raster("dist_mun.tif")
dmun=dmun*AOI
dd=raster("dist_d.tif")
dd=dd*AOI
dg=raster("dist_g.tif")
dg=dg*AOI
da=raster("dist_a.tif")
da=da*AOI
dpuo=raster("dist_puogis.tif")
dpuo=dpuo*AOI

# Raster Stack
biotstack=stack(dw,buh,suf,mun,d,hc,r,g,a,u,dpuo)
names(biotstack)=c("dw","buh","suf","mun","d","hc","r","g","a","u","dpuo")
parstack=stack(DG,nveg,siveg)
names(parstack)=c("DG","nveg","siveg")
biotdist=stack(dbuh,dsuf,dmun,dd,dg,da)
names(biotdist)=c("dbuh","dsuf","dmun","dd","dg","da")

plot(biotstack)
plot(parstack)
plot(biotdist)

# write Raster
{'
  writeRaster(dw,"rast_dw.tif",overwrite=T)
  writeRaster(buh,"rast_buh.tif",overwrite=T)
  writeRaster(suf,"rast_suf.tif",overwrite=T)
  writeRaster(mun,"rast_mun.tif",overwrite=T)
  writeRaster(d,"rast_d.tif",overwrite=T)
  writeRaster(hc,"rast_hc.tif",overwrite=T)
  writeRaster(r,"rast_r.tif",overwrite=T)
  writeRaster(g,"rast_g.tif",overwrite=T)
  writeRaster(a,"rast_a.tif",overwrite=T)
  writeRaster(u,"rast_u.tif",overwrite=T)
  writeRaster(dpuo,"rast_dpuo.tif",overwrite=T)
  writeRaster(HE,"rast_HE.tif",overwrite=T)
  writeRaster(nbiot,"rast_nbiot.tif",overwrite=T)
  writeRaster(cvbiot,"rast_cvbiot.tif",overwrite=T)
  writeRaster(sibiot,"rast_sibiot.tif",overwrite=T)
  writeRaster(DG,"rast_DG.tif",overwrite=T)
  writeRaster(nveg,"rast_nveg.tif",overwrite=T)
  writeRaster(siveg,"rast_siveg.tif",overwrite=T)
  writeRaster(cvveg,"rast_cvveg.tif",overwrite=T)
  '}


############################################################################
######### Background Point Generation ######################################

##### load GPS-locations
loc=shapefile("telemetrie_nsg.shp")
loc$kw=lubridate::week(loc$study_ts)
loc$anid=substr(loc$ind_ident,1,4)

##### create Background Points - rdm
# MCPs
{i1101=subset(loc,ind_ident=='1101_1734')
i1201=subset(loc,ind_ident=='1201_1946')
i1202=subset(loc,ind_ident=='1202_1948')
i1204=subset(loc,ind_ident=='1204_2146')
i1205=subset(loc,ind_ident=='1205_1695')
i1206=subset(loc,ind_ident=='1206_1945_2605')
i1207=subset(loc,ind_ident=='1207_2145')}

# single-MCPs to Polygone
{mcp1poly=mcp(i1101, percent=99.9)
mcp2poly=mcp(i1201, percent=99.9)
mcp3poly=mcp(i1202, percent=99.9)
mcp4poly=mcp(i1204, percent=99.9)
mcp5poly=mcp(i1205, percent=99.9)
mcp6poly=mcp(i1206, percent=99.9)
mcp7poly=mcp(i1207, percent=99.9)}


# buffer MCP Polygons
bufdist=603 # Buffer-Radius
{mcp1polybuf=raster::buffer(mcp1poly,width=bufdist)
mcp2polybuf=raster::buffer(mcp2poly,width=bufdist)
mcp3polybuf=raster::buffer(mcp3poly,width=bufdist)
mcp4polybuf=raster::buffer(mcp4poly,width=bufdist)
mcp5polybuf=raster::buffer(mcp5poly,width=bufdist)
mcp6polybuf=raster::buffer(mcp6poly,width=bufdist)
mcp7polybuf=raster::buffer(mcp7poly,width=bufdist)}

# rasterize buffered MCPs
{mcp1=rasterize(mcp1polybuf,mask)
mcp2=rasterize(mcp2polybuf,mask)
mcp3=rasterize(mcp3polybuf,mask)
mcp4=rasterize(mcp4polybuf,mask)
mcp5=rasterize(mcp5polybuf,mask)
mcp6=rasterize(mcp6polybuf,mask)
mcp7=rasterize(mcp7polybuf,mask)}


################################
### create Backgroung Points ###

{nrdm=2000 # number of Background Points
rdm1=randomPoints(mcp1,nrdm)
rdm2=randomPoints(mcp2,nrdm)
rdm3=randomPoints(mcp3,nrdm)
rdm4=randomPoints(mcp4,nrdm)
rdm5=randomPoints(mcp5,nrdm)
rdm6=randomPoints(mcp6,nrdm)
rdm7=randomPoints(mcp7,nrdm)}

# bind Background points
rdm = rbind(rdm1, rdm2, rdm3, rdm4, rdm5, rdm6, rdm7)
rdm = cbind(rdm, sort(rep(c(1101, 1201, 1202, 1204, 1205, 1206, 1207), nrdm)))
rdm = data.frame(rdm); names(rdm)=c("x","y","anid")

# write.table(rdm, "background_points_2k.csv", dec=".", sep=";", row.names=F)

rdm = read.csv("background_points_2k.csv", sep=";")

sprdm = SpatialPointsDataFrame(coords=rdm[,1:2], rdm, proj4string=CRS(proj4string(loc)))
mcp.area(sprdm, unout="km2", plotit=F)
mcp(sprdm, percent=99.9)

############################################################################
######### Point Extraction #################################################

## Extraction Rasterstack
predictors = stack(biotstack, parstack, biotdist)

## Extract Background Raster Information to Points
presvals = raster::extract(predictors, loc)
absvals = raster::extract(predictors, rdm[,1:2])

pb = c(rep(1, nrow(presvals)), rep(0, nrow(absvals)))
anid = c(substr(loc$ind_ident, 1, 4), sort(rep(c(1101, 1201, 1202, 1204, 1205, 1206, 1207), nrdm)))

##### Create RSF-Data Table
dat = data.frame(cbind(pb, anid, rbind(presvals, absvals)))
dat$anid=factor(dat$anid)

write.table(dat, "rsf_table.csv", dec=".", sep=";", row.names=F)