# ------ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ------ #
#     All the information which is not automatically computed, should be added here.           #     
#               This is done with the function KORAtool::KORAMonitoringBericht                 #
# ------ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ------ #

  #Trap night calendar information:

  effective.TN = params$effective.TN 
  potentially.TN= params$potentially.TN 
  #Title:
  Title= params$Title 
  Language= params$Language 
  #Authors: ex."Luc Le Grand, Florin Kunz, Fridolin Zimmermann"
  Authors= params$Authors
  #KORAphoto data csv name:
  KORAphoto= params$KORAphoto
  #Lynx master data csv name:
  Lynxmaster= params$Lynxmaster
  #Photo Predator csv name:
  PhotoPredator= params$PhotoPredator
  #Lynx Observation CSV name:
  LynxObs=params$LynxObs
  #Compartment: 
  #1:"Nordostschweiz"
  #2:"Misox (Mesolcina)"
  #3:"Tessin"
  #4:"Surselva"
  #5:"Zentralschweiz Ost"
  #6:"Rhone-Nord"
  #7:"Unterwallisüd"
  #8:"Simme-Saane"
  #9:"Engadin"
  #10:"Jura Nord"
  #11:"Jura Süd"
  #12:"Mittelbünden"
  #13:"Zentralschweiz West"
  #14:"Zentralschweiz Mitte"
  #15:"Berner Oberland Ost"
  #16:"Oberwallis")
  Compartment= params$Compartment
  #Sutdy Area (Reference Area):
  # 1: "Western central Switzerland"
  # 2: "North-eastern Switzerland"  
  # 3: "North-western Alps"
  # 4: "Mid-central Switzerland"
  # 5: "Southern Jura"
  # 6: "Central Jura"
  # 7: "Northern Jura" 
  # 8: "Northern Valais"  
  Refarea = params$Refarea
  #Suitable habitat (km2):
  favorable= params$favorable
  #Period length selected (3 other 5 days)
  Period.selected= params$Period.selected
  #JuvasMother
  JuvasMother= params$JuvasMother
  #Species.other to be ploted at the end of the report
  Species.other = params$Species.other
  # Which best Model (1 to 8 or 9 if ext.model provided)
  which.best.model.selected = params$which.best.model.selected
  #External Model results to be used in the report
  External.Model= params$External.Model
devtools::install_github("CrocutaLupus/KORAtool")
library(KORAtool)


library(data.table)
library(knitr)#tables
library(pander)
library(flextable)



# ------ Projection : ####
#Projection to be used for the map (CH1903 / LV03):
#Warnings OK
suppressWarnings(CRS<- sp::CRS("+init=epsg:21781"))
##### Import all the raw data needed ####

#KORA photo data:
Data<-data.table::fread(KORAphoto)
Data$TIME<-as.POSIXct(paste(Data$exposure_date,Data$exposure_time ,sep=" "), format= '%Y-%m-%d %H:%M:%S')

#Lynx master data table:
L.master<-data.table::fread(Lynxmaster)
L.master$mother<-as.character(L.master$mother)
L.master$yearOfBirth<-as.numeric(L.master$yearOfBirth)

#Text table for thanking part and method's foot note:
Text.table<-read.csv("Det_Mon_Report_Data_per_Reference_Area.csv")


##### Add sex and Juv info to Data: ####

#Select Juveniles's year of birth (previous year):
juv.year<-data.table::year(max(Data$TIME))-1

# Add info from Master Lynx Table: (L.master)

for(i in 1:length(L.master$name)){

  Data[Data$id_individual==L.master[i,name],"Sex"]<-L.master[i,sex]
  Data[Data$id_individual==L.master[i,name],"Juv"]<-L.master[i,yearOfBirth]

}

Data$Juv<-as.character(Data$Juv)

for(i in 1:nrow(Data)){
    if(Data[i ,Juv]==juv.year & !is.na(Data[i ,Juv])){
    Data[i ,"Juv"]<-"Juv"
  }
}

Data[Data$Sex=="","Sex"]<-NA
Data[Data$Juv!="Juv","Juv"]<-NA

#check if sex info is provided for every lynx
#unique(Data[is.na(Data$Sex) & Data$animal_species=="Lynx lynx",id_individual])
U.length<-length(Data[Data$id_individual=="U" &
       Data$animal_species=="Lynx lynx",file_name])

tot.length<-length(Data[Data$animal_species=="Lynx lynx",file_name])
#Table all sites:
Table.Sites<-data.table::data.table(SiteName=sort(unique(Data[,site_name])),
                                    Npictures=0,
                                    Npicutreslynx=0)


for(i in 1:length(Table.Sites$SiteName)){

  Table.Sites[i,"Npictures"]<-length(Data[Data$site_name==Table.Sites[i,SiteName],file_name])
  Table.Sites[i,"Npicutreslynx"]<-length(Data[Data$site_name==Table.Sites[i,SiteName] &
                                                Data$animal_species=="Lynx lynx",file_name])
}


Table.Sites$RunningNumber<-1:length(Table.Sites$SiteName)

Table.Sites<-flextable(Table.Sites)
Table.Sites <- set_caption(Table.Sites, caption = " Summary Sites Session. Please control in KORA photo that each site has at least one picture of something.")
Table.Sites<-font(Table.Sites,fontname = "Calibri", part = "all")
Table.Sites<-fontsize(Table.Sites,size = 12, part = "body")
Table.Sites <- bold(Table.Sites, bold = TRUE, part = "header")
Table.Sites<- autofit(Table.Sites) 
Table.Sites<-height_all(Table.Sites, height = .26)
Table.Sites<- hrule(Table.Sites, rule = "exact")
#Note: Individual considered as Juv are not plotted

KORAtool::KORAmapsex(KORA.Photo.Output=Data, 
  species="Lynx lynx",
  Buffer.map=0.03,
  Zoom.map=10,
  IDremove="U",
  Start= min(Data$TIME), #For the main map, we want all the lynx data plotted
  Stop=max(Data$TIME), #For the main map, we want all the lynx data plotted
  Pattern="TxxxxT",#pattern not resent but used not to block the function
  Buffer.label=2000, 
  Name.Map="Map_Report",
  Refarea=Refarea,
  Red.point.ID="U",
  Sexremove=c(""))
Start<-as.POSIXct(paste(Data[1,session_study_start_date],
                          Data[1,session_study_start_time],sep=" "),
                    format= '%Y-%m-%d %H:%M:%S')

Stop<-as.POSIXct(paste(Data[1,session_study_end_date],
                         Data[1,session_study_end_time],sep=" "),
                   format= '%Y-%m-%d %H:%M:%S')

#Is there any lynx that was seen before or after the 60 days? 

Seen.external<-unique(Data[Data$animal_species== "Lynx lynx",id_individual])[which(!(unique(Data[Data$animal_species== "Lynx lynx",id_individual]) %in%
unique(Data[Data$TIME>Start & Data$TIME<Stop & Data$animal_species== "Lynx lynx",id_individual])))]

#To check for a specific individual:
#x<-"R292"
#Before
#x %in% unique(Data[Data$TIME<Start & Data$animal_species== "Lynx lynx",id_individual])
#During
#x %in% unique(Data[Data$TIME>Start & Data$TIME<Stop & Data$animal_species== "Lynx lynx",id_individual])
#After
#x %in% unique(Data[Data$TIME>Start & Data$animal_species== "Lynx lynx",id_individual])

#Subset Data so that only information during the 60 days is used
DataFULL<-Data
Data<-Data[Data$TIME>Start & Data$TIME< Stop,]
#Number of unique sites (based on coordinates (uses all data also before/after 60 days)
Sites<-length(unique(paste(DataFULL$x,DataFULL$y,sep="")))

#Count Lynx
Lynx.Table<-data.table::data.table(ID=unique(Data[Data$animal_species=="Lynx lynx" &
                                   Data$id_individual!="U",id_individual]))


for(i in 1:length(Lynx.Table$ID)){
  Lynx.Table[i,"YOB"]<-L.master[L.master$name==Lynx.Table$ID[i],yearOfBirth]
  Lynx.Table[i,"Sex"]<-L.master[L.master$name==Lynx.Table$ID[i],sex]
  Lynx.Table[i,"mother"]<-L.master[L.master$name==Lynx.Table$ID[i],mother]
}

#Table with only adults
Lynx.Table.Adult<-Lynx.Table[Lynx.Table$YOB!=juv.year | is.na(Lynx.Table$YOB),]

#Count minimum adults
N.ind.lynx<-length(Lynx.Table.Adult$ID)

N.juv<-length(Lynx.Table[Lynx.Table$YOB==juv.year ,ID])



#Event:

## create column with ID of captured Individual, juvenile catches count as catch of the mother
Data$ID.mummy<-Data$id_individual

for(i in 1:length(Lynx.Table[Lynx.Table$YOB==juv.year,ID])){
  if(Lynx.Table[Lynx.Table$ID==Lynx.Table[Lynx.Table$YOB==juv.year,ID][i],mother]!=""){
    Data[Data$id_individual==Lynx.Table[Lynx.Table$YOB==juv.year,ID][i],"ID.mummy"]<-Lynx.Table[Lynx.Table$ID==Lynx.Table[Lynx.Table$YOB==juv.year,ID][i],mother]
  } #mother not known
}


#Select if Juveniles are counted as mother or not
if(JuvasMother==TRUE){
  Data$ID<-Data$ID.mummy
}
if(JuvasMother==FALSE){
  Data$ID<-Data$id_individual
}

# Create Event Table
Event.data<-Data[!grepl(paste(paste(Lynx.Table[Lynx.Table$YOB==juv.year,ID],collapse="|"),sep="|"),
                                         Data$ID) & Data$animal_species=="Lynx lynx",]
Event.data<-Event.data[Event.data$ID!="U",]

Event<-KORAtool::EvCountR(Table_Object=Event.data,Event_length=30)
Event.report<-Event

Event.report<-merge(x = Event.report, y = Lynx.Table, by="ID", all.x = TRUE)
Event.report$YOB<-as.character(Event.report$YOB)

Event.report[Event.report$Sex=="unknown","Sex"]<-NA
Event.report[Event.report$Sex=="male","Sex"]<-"M"
Event.report[Event.report$Sex=="female","Sex"]<-"F"

# --- Add first seen

#Import Raw Data (CSV from KORA Data): 
  LynxMaster.table<-suppressWarnings(data.table::fread(Lynxmaster))
  LynxObs.table<-data.table::fread(LynxObs,select=c("lynx_ID","image_uid","date","x","y"))
  PhotoPredator.table<-data.table::fread(PhotoPredator,select=c("id_individual","exposure_date","session_name","x","y"))
#Select interesting Columns:  
  First.seen.table<-LynxMaster.table[,c("name","coatPattern","yearOfBirth","deathDate","mother")]

#Fill in info First/Last seen: 
 for(i in 1:length(First.seen.table$name)){
   First.seen.table[i,"FirstSeen"]<-suppressWarnings(min(min(as.POSIXct(LynxObs.table[LynxObs.table$lynx_ID==as.character(First.seen.table[i,name]),date],format='%Y-%m-%d',tz="")),
                                   min(as.POSIXct(PhotoPredator.table[PhotoPredator.table$id_individual==as.character(First.seen.table[i,name]),exposure_date],format='%Y-%m-%d'))))
   First.seen.table[i,"LastSeen"]<- suppressWarnings(max(max(as.POSIXct(LynxObs.table[LynxObs.table$lynx_ID==as.character(First.seen.table[i,name]),date],format='%Y-%m-%d',tz="")),
                                   max(as.POSIXct(PhotoPredator.table[PhotoPredator.table$id_individual==as.character(First.seen.table[i,name]),exposure_date],format='%Y-%m-%d'))))
}

Event.report<-merge(Event.report, First.seen.table[,c("name","FirstSeen")],all.x = TRUE,by.x = "ID",by.y = "name")

#Check if some lynx were discovered during the session:

for(i in 1:nrow(Event.report)){
  Event.report[i,"FirstTimeSeen"]<-Event.report[i,FirstSeen]>as.POSIXct(paste(unique(Data$session_start_date),unique(Data$session_start_time)),format="%Y-%m-%d %H:%M:%S")
}

New.lynx<-Event.report[Event.report$FirstTimeSeen==TRUE,ID]

Event.report$FirstTimeSeen<-NULL

#First seen only show year in final table
Event.report$FirstSeen<-as.character(year(Event.report$FirstSeen))

# --- Add juveniles info to each known mother
for(i in 1:length(unique(Lynx.Table[Lynx.Table$YOB==juv.year,mother]))){
  Event.report[Event.report$ID==unique(Lynx.Table[Lynx.Table$YOB==juv.year,mother])[i],"Juv"]<-paste(
    Lynx.Table[Lynx.Table$mother==unique(Lynx.Table[Lynx.Table$YOB==juv.year,mother])[i] & !is.na(Lynx.Table$mother) & Lynx.Table$YOB==juv.year,ID],collapse=", ")
}



#Add canton seen:
for(i in 1:length(Event.report$ID)){

  if(JuvasMother==TRUE){
   Event.report[i,"Canton"]<-paste(unique(Data[Data$ID.mummy==Event.report[i,ID],loc_canton]),collapse = ", ")
  }  

   if(JuvasMother==FALSE){
   Event.report[i,"Canton"]<-paste(unique(Data[Data$id_individual==Event.report[i,ID],loc_canton]),collapse = ", ")
  }  

}


#Rename Table
Event.report<-Event.report[,c("ID","Event","N pictures","N site","YOB","FirstSeen","Sex","mother","Juv","Canton")]


#Remove some columns (discussion with Florin and Fridolin 2021-09-01)
Event.report<-Event.report[,-c("N pictures","N site","YOB")]

#Rename columns depending of language:
if(Language=="FR"){names(Event.report)<-c("ID","Événement(s)","Connu depuis","Sexe","Mère","Juvénile(s)","Canton(s)")}
if(Language=="DE"){names(Event.report)<-c("ID","Ereignisse","Bekannt seit","Geschlecht","Mutter","Juvenil(e)","Kanton(e)")}

# Change Female ("F") to w
if(Language=="DE"){
  Event.report[Event.report$Geschlecht=="F","Geschlecht"]<-"W"
}

#Total as last:
#bindtotal<-Event.report[Event.report$ID=="Total",]
#Event.report<-Event.report[Event.report$ID!="Total",]
#Event.report<-rbind(Event.report,bindtotal)

#Remove Total from Table
Event.report<-Event.report[Event.report$ID!="Total",]

#New lynx in bold (get line number for later)
New.lynx<-Event.report[,which(Event.report$ID %in% New.lynx)]

# Data table into flextable to plot:



Event.report<-flextable(Event.report)

if(Language=="FR"){
 Event.report <- set_caption(Event.report, caption = "Lynx indépenants détectés dans l’aire de référence pendant les 60 nuits d’échantillonnage. Un événement comprend toutes les images d'un individu à un site donné séparées de moins de 30 minutes. Connu depuis\ua0: année de la première observation. Mère\ua0: mère des lynx indépendants quand elle est connue. Juvénile(s)\ua0: juvénile(s) des lynx femelles indépendantes quand elles sont connues. Canton(s)\ua0: canton(s) dans le(s)quels les lynx ont été observés lors de la session. Les lynx qui ont été détectés pour la première fois lors de la session sont mis en évidence en **gras**.",autonum = officer::run_autonum(
  seq_id = "table",
  pre_label = "Tab. ",
  post_label = ". ",
  bkm = NULL,
  bkm_all = FALSE,
  prop = officer::fp_text(bold= TRUE, font.size = 10, font.family= "Calibri") ,
  start_at = ""
)) 
}

if(Language=="DE"){
 Event.report <- set_caption(Event.report, caption = "selbständige Luchse, die in den 60 Nächten des deterministischen Durchgangs fotografiert wurden. Ein Ereignis umfasst alle Bilder eines Luchsindividuums am selben Standort, welche weniger als 30 Minuten auseinander liegen. Bekannt seit: Jahr des ersten Nachweises. Mutter: Mutter von selbständigen Luchsen, sofern bekannt. Juvenil(e): Jungtiere von selbständigen Luchsen, sofern bekannt. Kanton(e): Kantone, in welchen sie während des Durchgangs fotografiert wurden. **Fett** hervorgehoben: Luchse, die während des Durchgangs zum ersten Mal nachgewiesen wurden.",autonum = officer::run_autonum(
  seq_id = "table",
  pre_label = "Tab. ",
  post_label = ". ",
  bkm = NULL,
  bkm_all = FALSE,
  prop = officer::fp_text(bold= TRUE, font.size = 10, font.family= "Calibri") ,
  start_at = ""
)) 
}

# --- Add footnote:

#Texte juvéniles sans mère:

if(Language== "FR"){
  txt.juv<-"Les juvénile(s) sans mère connue suivants n'apparaissent pas dans le tableau\ua0: "

  #Some Juv. have mothers which are not known
  if(length(Lynx.Table[Lynx.Table$YOB==juv.year & Lynx.Table$mother=="",ID])>0){
    for(i in 1:length(Lynx.Table[Lynx.Table$YOB==juv.year & Lynx.Table$mother=="",ID])){

      txt.juv<- paste(txt.juv,
                      Lynx.Table[Lynx.Table$YOB==juv.year & Lynx.Table$mother=="",ID][i]," (Nb photo = ",
                      length(Data[Data$id_individual==Lynx.Table[Lynx.Table$YOB==juv.year & Lynx.Table$mother=="",ID][i],file_name]),
                      ") ",sep="")

    } 
  #All mothers are known:
  }else{

    #Not all mothers have been seen
    if(!all(Lynx.Table[Lynx.Table$YOB==juv.year,mother] %in% Lynx.Table$ID)){

      mom.present<-Lynx.Table[Lynx.Table$YOB==juv.year,mother] %in% Lynx.Table$ID


      if(length(mom.present[mom.present == FALSE])==1){
        txt.juv<-paste("Le juvénile ", Lynx.Table[Lynx.Table$YOB==juv.year,ID][!mom.present], " a été observé sans sa mère ", Lynx.Table[Lynx.Table$YOB==juv.year,mother][!mom.present],".",sep="")
      }else{
        txt.juv<-paste("Juvéniles observés sans mère\ua0: ",paste(Lynx.Table[Lynx.Table$YOB==juv.year,ID],collapse=", "))
      }
    #All mothers have been seen  
    }else{txt.juv<-"Tous les juvéniles photographiés ont pu être associé à leur mère. Ils figurent donc tous dans le tableau."}

  } 
}


if(Language== "DE"){
  txt.juv<-"Folgende(r) Jungluchse(r) ohne bekannte Mutter erscheint/erscheinen nicht in der Tabelle: "

  #Some Juv. have mothers which are not known
  if(length(Lynx.Table[Lynx.Table$YOB==juv.year & Lynx.Table$mother=="",ID])>0){
    for(i in 1:length(Lynx.Table[Lynx.Table$YOB==juv.year & Lynx.Table$mother=="",ID])){

      txt.juv<- paste(txt.juv,
                      Lynx.Table[Lynx.Table$YOB==juv.year & Lynx.Table$mother=="",ID][i]," (Anzahl Bilder = ",
                      length(Data[Data$id_individual==Lynx.Table[Lynx.Table$YOB==juv.year & Lynx.Table$mother=="",ID][i],file_name]),
                      ") ",sep="")

    } 
    #All mothers are known:
  }else{

    #Not all mothers have been seen
    if(!all(Lynx.Table[Lynx.Table$YOB==juv.year,mother] %in% Lynx.Table$ID)){

      mom.present<-Lynx.Table[Lynx.Table$YOB==juv.year,mother] %in% Lynx.Table$ID


      if(length(mom.present[mom.present == FALSE])==1){
        txt.juv<-paste("Der juvenile ", Lynx.Table[Lynx.Table$YOB==juv.year,ID][!mom.present], " wurde ohne seine Mutter (",Lynx.Table[Lynx.Table$YOB==juv.year,mother][!mom.present],") beobachtet.",sep="")
      }else{
        txt.juv<-paste("Jungtiere, die ohne Mutter beobachtet wurden: ",paste(Lynx.Table[Lynx.Table$YOB==juv.year,ID],collapse=", "))
      }
      #All mothers have been seen  
    }else{txt.juv<-"Von allen fotografierten Jungtieren konnte die Mutter eruiert werden. Somit erscheinen alle Jungtiere in der Tabelle."}

  } 
}

#Text JuvasMother:

if(Language=="FR"){
 if(JuvasMother==TRUE){

txt.juvasmother<-"A cause de leur faible taux de détection et leur fort taux de disparition (mortalité et dispersion) les individus juvéniles photographiés au cours du monitoring sont identifiés et comptabilisés comme une capture de leur mère dans le calendrier de capture." 

}else{

txt.juvasmother<-"A cause de leur faible taux de détection et leur fort taux de disparition (mortalité et dispersion) les individus juvéniles photographiés au cours du monitoring sont identifiés, mais pas comptabilisés comme individus dans le calendrier de capture."

} 
}

if(Language=="DE"){
 if(JuvasMother==TRUE){

txt.juvasmother<-"Aufgrund ihrer geringen Erfassbarkeit und hohen Verschwinderate (Mortalität und Dispersal) werden Jungtiere, die während des Durchgangs fotografiert wurden, identifiziert und im Fangkalender als Fang ihrer Mutter eingetragen." 

}else{

txt.juvasmother<-"Aufgrund ihrer geringen Erfassbarkeit und hohen Verschwinderate (Mortalität und Dispersal) werden Jungtiere, die während des Durchgangs fotografiert wurden, identifiziert aber nicht als Individuen im Fangkalender berücksichtigt."

} 
}


#Create footnote:
Event.report<- footnote( Event.report, i = 1,j=c(1,6),
            value = as_paragraph(
                      c(txt.juv,
                        txt.juvasmother)),
            ref_symbols = c("1","2"),
            part = "header")

Event.report <- bold(Event.report, bold = TRUE, part = "footer")

#New lynx in bold (get line number for later)
Event.report <- bold(Event.report, bold = TRUE, part = "body",i=New.lynx)
Event.report <- bold(Event.report, bold = TRUE, part = "header")

Event.report<-font(Event.report,fontname = "Calibri", part = "all")
Event.report<-fontsize(Event.report,size = 12, part = "body")

Event.report<- autofit(Event.report) 
Event.report<-height_all(Event.report, height = .26)
Event.report<- hrule(Event.report, rule = "exact")
Event.report<-fontsize(Event.report, size = 8, part = "footer")
# KORAmapopp using LynxObs (Opp. mon observations)


Plot.Opp.Mon<-KORAtool::KORAmapopp(LynxObs= LynxObs,
                              Start=Start,
                              Stop=Stop,
                              Compartment=Compartment,
                              Refarea=Refarea,
                              IDremove=c("U","lx","XXXX"),
                              Buffer.polygon=2000)

ggplot2::ggsave("Plot_Opp_Mon.jpeg",plot=Plot.Opp.Mon,
                    units = "cm",
                    width = 16.5,
                    height = 20)

# -- Plot function Other species (used in appendix)

Plot_KORA_Report<-function(Species){
  #Import shapefile compartment
Komp <- sf::read_sf(dsn = "MAP_Data/Luchs_Komp_21_07_2015_new.shp")

# -- Study Area #
study_area<-sp::bbox(as.matrix(Komp$geometry[Compartment][[1]]))
study_area.origin<-sp::bbox(as.matrix(Komp$geometry[Compartment][[1]]))
#Add x% around the BBox to have some extra map area

  study_area[1,1]<-study_area[1,1]-round(study_area[1,1]*0.02) #x min
  study_area[2,1]<-study_area[2,1]-round(study_area[2,1]*6*0.02) #y min
  study_area[1,2]<-study_area[1,2]+round(study_area[1,2]*0.02) #x max
  study_area[2,2]<-study_area[2,2]+round(study_area[2,2]*0.02) #y max

study_area<-rgeos::readWKT(paste("POLYGON((",
                                   study_area[1,1]," ",study_area[2,1],",",
                                   study_area[1,1]," ",study_area[2,2],",",
                                   study_area[1,2]," ",study_area[2,2],",",
                                   study_area[1,2]," ",study_area[2,1],",",
                                   study_area[1,1]," ",study_area[2,1],"))",sep=""),p4s= CRS)



# --- Import shapefile Study area Polygon
suppressWarnings(Rcompartment <- raster::shapefile("MAP_Data/Reference_areas_all_new2022.shp"))
Rcompartment<-Rcompartment[Refarea,]

#Import shapefile compartment
suppressWarnings(Komp <- raster::shapefile("MAP_Data/Luchs_Komp_21_07_2015_new.shp"))
Komp<-Komp[Compartment,]

# ------ Open Source Map

#study area BBox in lat long
study_area_lat_long <- sp::spTransform(study_area, sp::CRS("+init=epsg:4326"))

  LAT1 = study_area_lat_long@bbox[2,1] ; LAT2 = study_area_lat_long@bbox[2,2]
  LON1 = study_area_lat_long@bbox[1,1] ; LON2 = study_area_lat_long@bbox[1,2]

#Get the map
#warnings OK
suppressWarnings(map <- OpenStreetMap::openmap(c(LAT2,LON1), c(LAT1,LON2), zoom = NULL, #can be replaced by NULL
                                                 type = c("stamen-terrain")[1],
                                                 mergeTiles = TRUE))

#Correct projection
#warnings OK
suppressWarnings(map <- OpenStreetMap::openproj(map, projection = "+init=epsg:21781"))

# ------ Build map 

# --- Open Source Map in ggplot:
map<-OpenStreetMap::autoplot.OpenStreetMap(map)+
ggplot2::theme(axis.title=ggplot2::element_blank(),
                   axis.text=ggplot2::element_blank(),
                   axis.ticks=ggplot2::element_blank(),
                   panel.border = ggplot2::element_rect(colour = "white", fill=NA, size=20))


# --- Add Scale:

  # distance on x axes:
  dist.scale<-plyr::round_any(round(((study_area@bbox[1,2]-study_area@bbox[1,1])/1000)/4), 10, f = ceiling)

  # scale thickness:
  s.thick<-(0.01*(study_area@bbox[2,2]-study_area@bbox[2,1]))

  # scale black rectangle
  xleft<-study_area@bbox[1,1]+((study_area.origin[1,1]-study_area@bbox[1,1])/3)#left
  xright<-xleft+dist.scale*1000#right
  ybottom<-study_area.origin[2,1]-4*s.thick#bottom
  ytop<-study_area.origin[2,1]-3*s.thick#top

  map<-map+
    ggplot2::geom_rect(mapping=ggplot2::aes(xmin=xleft, xmax=xright, ymin=ybottom, ymax=ytop),
                       fill=c("black"),
                       inherit.aes = FALSE)+
    ggplot2::geom_text(x=xright+(2/5)*dist.scale*1000, y=ytop, label=paste(dist.scale,"Km",sep=" "), cex=10, color ="black")


# --- Add other layers:

map<-map+
#KOMP
ggplot2::geom_polygon(data = Komp@polygons[[1]], ggplot2::aes(x=long, y=lat,group = group),
                               colour="black", fill="white", alpha=0.4, lwd=1.1)+
#Refarea
ggplot2::geom_polygon(data = Rcompartment@polygons[[1]], ggplot2::aes(x=long, y=lat,group = group),
                               colour="darkblue", fill=NA, alpha=1, lwd=1.1)+
#Points Sites
ggplot2::geom_point(data=DataFULL,ggplot2::aes(x,y), col="white", pch=19,cex=3.5)+
      ggplot2::geom_point(data=DataFULL,ggplot2::aes(x,y),col="black", pch=1,cex=3.5)+
#Points Sites positive for species
ggplot2::geom_point(data=DataFULL[DataFULL$animal_species==Species,],ggplot2::aes(x,y), col="black", pch=19,cex=1.5)



ggplot2::ggsave(paste(Species,".jpeg",sep=""),plot=map,
                width = 16.5,
                    height = 17,
                units = "cm")

}
Plot_KORA_Report_2_SP<-function(Species1,Species2){
  #Import shapefile compartment
  Komp <- sf::read_sf(dsn = "MAP_Data/Luchs_Komp_21_07_2015_new.shp")

  # -- Study Area #
  study_area<-sp::bbox(as.matrix(Komp$geometry[Compartment][[1]]))
  study_area.origin<-sp::bbox(as.matrix(Komp$geometry[Compartment][[1]]))
  #Add x% around the BBox to have some extra map area

  study_area[1,1]<-study_area[1,1]-round(study_area[1,1]*0.02) #x min
  study_area[2,1]<-study_area[2,1]-round(study_area[2,1]*6*0.02) #y min
  study_area[1,2]<-study_area[1,2]+round(study_area[1,2]*0.02) #x max
  study_area[2,2]<-study_area[2,2]+round(study_area[2,2]*0.02) #y max

  study_area<-rgeos::readWKT(paste("POLYGON((",
                                   study_area[1,1]," ",study_area[2,1],",",
                                   study_area[1,1]," ",study_area[2,2],",",
                                   study_area[1,2]," ",study_area[2,2],",",
                                   study_area[1,2]," ",study_area[2,1],",",
                                   study_area[1,1]," ",study_area[2,1],"))",sep=""),p4s= CRS)



  # --- Import shapefile Study area Polygon
  suppressWarnings(Rcompartment <- raster::shapefile("MAP_Data/Reference_areas_alln.shp"))
  Rcompartment<-Rcompartment[Refarea,]

  #Import shapefile compartment
  suppressWarnings(Komp <- raster::shapefile("MAP_Data/Luchs_Komp_21_07_2015_new.shp"))
  Komp<-Komp[Compartment,]

  # ------ Open Source Map

  #study area BBox in lat long
  study_area_lat_long <- sp::spTransform(study_area, sp::CRS("+init=epsg:4326"))

  LAT1 = study_area_lat_long@bbox[2,1] ; LAT2 = study_area_lat_long@bbox[2,2]
  LON1 = study_area_lat_long@bbox[1,1] ; LON2 = study_area_lat_long@bbox[1,2]

  #Get the map
  #warnings OK
  suppressWarnings(map <- OpenStreetMap::openmap(c(LAT2,LON1), c(LAT1,LON2), zoom = NULL, #can be replaced by NULL
                                                 type = c("stamen-terrain")[1],
                                                 mergeTiles = TRUE))

  #Correct projection
  #warnings OK
  suppressWarnings(map <- OpenStreetMap::openproj(map, projection = "+init=epsg:21781"))

  # ------ Build map 

  # --- Open Source Map in ggplot:
  map<-OpenStreetMap::autoplot.OpenStreetMap(map)+
    ggplot2::theme(axis.title=ggplot2::element_blank(),
                   axis.text=ggplot2::element_blank(),
                   axis.ticks=ggplot2::element_blank(),
                   panel.border = ggplot2::element_rect(colour = "white", fill=NA, size=30),
                   legend.position="none")


  # --- Add Scale:

  # distance on x axes:
  dist.scale<-plyr::round_any(round(((study_area@bbox[1,2]-study_area@bbox[1,1])/1000)/4), 10, f = ceiling)

  # scale thickness:
  s.thick<-(0.01*(study_area@bbox[2,2]-study_area@bbox[2,1]))

  # scale black rectangle
  xleft<-study_area.origin[1,1]#left
  xright<-xleft+dist.scale*1000#right
  ybottom<-study_area.origin[2,1]-4*s.thick#bottom
  ytop<-study_area.origin[2,1]-3*s.thick#top

  map<-map+
    ggplot2::geom_rect(mapping=ggplot2::aes(xmin=xleft, xmax=xright, ymin=ybottom, ymax=ytop),
                       fill=c("black"),
                       inherit.aes = FALSE)+
    ggplot2::geom_text(x=xright+(2/5)*dist.scale*1000, y=ytop, label=paste(dist.scale,"Km",sep=" "), cex=10, color ="black")


  # --- Add other layers:

  map<-map+
    #KOMP
    ggplot2::geom_polygon(data = Komp@polygons[[1]], ggplot2::aes(x=long, y=lat,group = group),
                          colour="black", fill="white", alpha=0.4, lwd=1.1)+
    #Refarea
    ggplot2::geom_polygon(data = Rcompartment@polygons[[1]], ggplot2::aes(x=long, y=lat,group = group),
                          colour="darkblue", fill=NA, alpha=1, lwd=1.1)+
    #Points Sites
    ggplot2::geom_point(data=DataFULL,ggplot2::aes(x,y), col="white", pch=19,cex=6)+
    ggplot2::geom_point(data=DataFULL,ggplot2::aes(x,y),col="black", pch=1,cex=6)


  # --- Add Pie charts

  #Create table with all sites
  Pie.table<-data.table(Site=unique(paste(DataFULL$x,DataFULL$y,sep="/")))
  Pie.table$x<-substr(Pie.table$Site,1,6)                   
  Pie.table$y<-substr(Pie.table$Site,8,13)
  Pie.table$x<-as.numeric(Pie.table$x)
  Pie.table$y<-as.numeric(Pie.table$y)

  #Fill in Event Info per site
  suppressWarnings(for(i in 1:nrow(Pie.table)){

  #Lynx events (ID of lynx is not used to prevent over estimation in comparaison of wolf observations)    
  Loop<-DataFULL[DataFULL$animal_species==Species1 &
                 DataFULL$x==as.numeric(Pie.table[i,x]) &
                 DataFULL$y==as.numeric(Pie.table[i,y]),c("id_individual","TIME","x","y")]
  Loop$id_individual<-""
  names(Loop)[1]<-"ID"
  Loop<-KORAtool::EvCountR(Loop,30)
  Pie.table[i,"Event.Lynx"]<- Loop[Loop$ID=="Total",Event]

  #Wolf events (does not take wolf/Species2 numbers into account)
  Loop<-DataFULL[DataFULL$animal_species==Species2 &
                   DataFULL$x==as.numeric(Pie.table[i,x]) &
                   DataFULL$y==as.numeric(Pie.table[i,y]),c("id_individual","TIME","x","y")]
  names(Loop)[1]<-"ID"
  Loop<-KORAtool::EvCountR(Loop,30)
  Pie.table[i,"Event.Wolf"]<- Loop[Loop$ID=="Total",Event]

  })

  #Add to the map
  map<-map+
  scatterpie::geom_scatterpie(ggplot2::aes(x=x, y=y), data=Pie.table, cols=c("Event.Lynx","Event.Wolf"),pie_scale = 1.1)

  print(map)

}

#jpeg(filename = "Plot_Lynx_Wolf_Surselva.jpeg", width = 16.5, height = 15,
#     units = "cm", quality = 100, bg = "white", res=1000)

#Plot_KORA_Report_2_SP("Lynx lynx","Canis lupus")

#dev.off()
# Create table occasion:
Occ.table<-data.table(ID=unique(Lynx.Table.Adult$ID))


# To carry the analyses of the lynx monitoring, we test periods of 3 days and periods of 5 days

for(Period in c(3,5)){

# ---- Create Empty Occurrence Matrix  
occasions<-seq(from = Start,to = Stop, by = paste(Period,"days", sep=" "))
Occ.table<-data.table::data.table(ID = unique(Lynx.Table.Adult$ID))
if(Period==3){
  Occ.table[,c(as.character(occasions[-21])):= 0]
}
if(Period==5){
  Occ.table[,c(as.character(occasions[-13])):= 0]
}

# ---- Fill in Matrix 

#observation of juvenile counts as an observation of the mother

if(JuvasMother==TRUE){
  for(i in 1:length(Occ.table$ID)){
  for(j in 2:ncol(Occ.table)){
    Data_loop<- Data[Data$TIME > as.POSIXct(names(Occ.table)[j])  &
                         Data$TIME < as.POSIXct(names(Occ.table)[j]) + Period*86400 & #select days occasion
                         Data$ID.mummy == Occ.table[i,ID],c("TIME","file_name")]
    if(length(Data_loop$file_name)!=0){
      Occ.table[i,j]<-1
    }else{Occ.table[i,j]<-0}
  }
}
}

#observation of juvenile does not count as an observation of the mother

if(JuvasMother==FALSE){
  for(i in 1:length(Occ.table$ID)){
  for(j in 2:ncol(Occ.table)){
    Data_loop<- Data[Data$TIME > as.POSIXct(names(Occ.table)[j])  &
                         Data$TIME < as.POSIXct(names(Occ.table)[j]) + Period*86400 & #select days occasion
                         Data$id_individual == Occ.table[i,ID],c("TIME","file_name")]
    if(length(Data_loop$file_name)!=0){
      Occ.table[i,j]<-1
    }else{Occ.table[i,j]<-0}
  }
}
}


# ---- Check if an individual has no single catch within the study period (60nights)
for(i in 1:length(Occ.table$ID)){
  if(sum(Occ.table[1,2:ncol(Occ.table)])==0){
    warning("warning: ID has not been seen during the 60 days.")
  }
}


# ---- Occ.table for Plot cumulative (2 table, one per period length)

if(Period==3){
Occ.table.3<-Occ.table
Occ.table.plot.3<-dcast(melt(Occ.table, id.vars = "ID"), variable ~ ID)
names(Occ.table.plot.3)[1]<-"TIME"
#number of captures per occasion
Occ.table.plot.3$SUM.occasion<-rowSums(Occ.table.plot.3[,-1])
#cumulative number of captures per occasion
Occ.table.plot.3$SUM.occasion.cumulative<-cumsum(Occ.table.plot.3$SUM.occasion)

#Cumulative number of ID lynx

Occ.table.plot.3$SUM.ID.cumulative<-sum(Occ.table.plot.3[1,2:(length(Occ.table.plot.3)-2)])
Sum.ID<-as.numeric(Occ.table.plot.3[1,2:(length(Occ.table.plot.3)-3)])

for(i in 2:20){
 Sum.ID<-Sum.ID+as.numeric(Occ.table.plot.3[i,2:(length(Occ.table.plot.3)-3)])

 Sum.ID<-replace(Sum.ID, Sum.ID==2, 1)

 Occ.table.plot.3[i,"SUM.ID.cumulative"]<-sum(Sum.ID)

}

}

if(Period==5){
Occ.table.5<-Occ.table
Occ.table.plot.5<-dcast(melt(Occ.table, id.vars = "ID"), variable ~ ID)
names(Occ.table.plot.5)[1]<-"TIME"
#number of captures per occasion
Occ.table.plot.5$SUM.occasion<-rowSums(Occ.table.plot.5[,-1])
#cumulative number of captures per occasion
Occ.table.plot.5$SUM.occasion.cumulative<-cumsum(Occ.table.plot.5$SUM.occasion)

#Cumulative number of ID lynx

Occ.table.plot.5$SUM.ID.cumulative<-sum(Occ.table.plot.5[1,2:(length(Occ.table.plot.5)-2)])
Sum.ID<-as.numeric(Occ.table.plot.5[1,2:(length(Occ.table.plot.5)-3)])

for(i in 2:12){
 Sum.ID<-Sum.ID+as.numeric(Occ.table.plot.5[i,2:(length(Occ.table.plot.5)-3)])

 Sum.ID<-replace(Sum.ID, Sum.ID==2, 1)

 Occ.table.plot.5[i,"SUM.ID.cumulative"]<-sum(Sum.ID)

}

}

}
# ---- Closed Test:

for(Period in c(3,5)){

#Create table Sites info + read trap
Lynx_trapXY<-unique(Data[,c("site_name","x","y")])
names(Lynx_trapXY)[1]<-"trapID"
traps.lynx <- secr::read.traps(data = Lynx_trapXY)

#Create table observation info
Lynx_captXY<-data.table(Session=1,
                        ID=rep(Lynx.Table.Adult$ID,each=(60/Period)*length(Lynx_trapXY$trapID)),
                        Occasion=rep(rep(1:(60/Period),each=length(Lynx_trapXY$trapID)),length(Lynx.Table.Adult$ID)),
                        X=rep(Lynx_trapXY$x),
                        X=rep(Lynx_trapXY$y),
                        TrapID=rep(Lynx_trapXY$trapID),
                        Seen="NO")

if(Period==3){
  occasions<-seq(from = Start,to = Stop, by = paste(Period,"days", sep=" "))
  occasions<-occasions[-21]

  #Fill in information if lynx seen during period at given site
  for(i in 1:length(Lynx_captXY$Session)){

  if(JuvasMother==TRUE){
    Data_loop<- Data[Data$TIME > as.POSIXct(occasions[Lynx_captXY[i,Occasion]])  &
                         Data$TIME < as.POSIXct(occasions[Lynx_captXY[i,Occasion]]) + Period*86400 & #select days occasion
                         Data$ID.mummy == Lynx_captXY[i,ID] &
                         Data$site_name== Lynx_captXY[i,TrapID],c("TIME","file_name")]
  }

  if(JuvasMother==FALSE){
    Data_loop<- Data[Data$TIME > as.POSIXct(occasions[Lynx_captXY[i,Occasion]])  &
                         Data$TIME < as.POSIXct(occasions[Lynx_captXY[i,Occasion]]) + Period*86400 & #select days occasion
                         Data$id_individual == Lynx_captXY[i,ID] &
                         Data$site_name== Lynx_captXY[i,TrapID],c("TIME","file_name")]
  }

  if(length(Data_loop$file_name)>0){
    Lynx_captXY[i,"Seen"]<-"YES"
  }
  }

}
if(Period==5){

  #Occasions starting data/time
  occasions<-seq(from = Start,to = Stop, by = paste(Period,"days", sep=" "))
  occasions<-occasions[-13]

  #Fill in information if lynx seen during period at given site
  for(i in 1:length(Lynx_captXY$Session)){

  if(JuvasMother==TRUE){
  Data_loop<- Data[Data$TIME > as.POSIXct(occasions[Lynx_captXY[i,Occasion]])  &
                         Data$TIME < as.POSIXct(occasions[Lynx_captXY[i,Occasion]]) + Period*86400 & #select days occasion
                         Data$ID.mummy == Lynx_captXY[i,ID] &
                         Data$site_name== Lynx_captXY[i,TrapID],c("TIME","file_name")]
  }

  if(JuvasMother==FALSE){
  Data_loop<- Data[Data$TIME > as.POSIXct(occasions[Lynx_captXY[i,Occasion]])  &
                         Data$TIME < as.POSIXct(occasions[Lynx_captXY[i,Occasion]]) + Period*86400 & #select days occasion
                         Data$id_individual == Lynx_captXY[i,ID] &
                         Data$site_name== Lynx_captXY[i,TrapID],c("TIME","file_name")]
  }  


  if(length(Data_loop$file_name)>0){
    Lynx_captXY[i,"Seen"]<-"YES"
  }
  }


}



Lynx_captXY<-Lynx_captXY[Lynx_captXY$Seen=="YES",-c("TrapID","Seen")]

# --- Closure Test

if(Period==3){
  # ---- Create capthist 
  suppressWarnings(Lynx.CHxy.3  <- secr::make.capthist (Lynx_captXY, traps.lynx, fmt = "XY"))
  #warnings OK it produces the same results as in the program

  # ---- Closure test
  suppressWarnings(closure.result.3<-secr::closure.test(Lynx.CHxy.3,SB=TRUE))
  #warnings if pop too small

  # ---- Output for Mark program

  #Empty table
  inp.Mark<-data.table::data.table(ID=Lynx.Table.Adult$ID,ch="empty")
  #needs data.frame
  Occ.table.plot.3.1<-as.data.frame(Occ.table.plot.3)
  for(i in 1:length(Lynx.Table.Adult$ID)){
    inp.Mark[i,"ch"]<-paste( Occ.table.plot.3.1[,inp.Mark$ID[i]] ,collapse = "")
  }
  #process (RMark format)
  inp.Mark=RMark::process.data(inp.Mark)
  #output
  RMark::export.chdata(data=inp.Mark, filename=paste(ReportName,"period_3",sep="_"), covariates = NULL, replace = TRUE)

}

if(Period==5){
  # ---- Create capthist  
  suppressWarnings(Lynx.CHxy.5  <- secr::make.capthist (Lynx_captXY, traps.lynx, fmt = "XY"))
  #warnings OK it produces the same results as in the program

  # ---- Closure test
  suppressWarnings(closure.result.5<-secr::closure.test(Lynx.CHxy.5,SB=TRUE))
  #warnings if pop too small - do it any way like in program

  # ---- Output for Mark program

  #Empty table
  inp.Mark<-data.table::data.table(ID=Lynx.Table.Adult$ID,ch="empty")
  #needs data.frame
  Occ.table.plot.5.1<-as.data.frame(Occ.table.plot.5)
  for(i in 1:length(Lynx.Table.Adult$ID)){
    inp.Mark[i,"ch"]<-paste( Occ.table.plot.5.1[,inp.Mark$ID[i]] ,collapse = "")
  }
  #process (RMark format)
  inp.Mark=RMark::process.data(inp.Mark)
  #output
  RMark::export.chdata(data=inp.Mark, filename=paste(ReportName,"period_5",sep="_"), covariates = NULL, replace = TRUE)
}

}
# ----------------------------- !!!!!----------------------------------------- #
# It is assumed that MARK.EXE is located in directory "C:/Program Files/Mark". #
# https://cran.r-project.org/web/packages/RMark/RMark.pdf                      #
# ----------------------------- !!!!!----------------------------------------- #



#Capture History Period=3:
CaptHist.table.3<-data.table(ID=Occ.table.3$ID)
for( i in 1:length(Occ.table.3$ID)){
CaptHist.table.3[i,"ch"]<-paste(as.numeric(Occ.table.3[i,2:length(Occ.table.3)]),collapse = "")
}

#Capture History Period=5:
CaptHist.table.5<-data.table(ID=Occ.table.5$ID)
for( i in 1:length(Occ.table.5$ID)){
CaptHist.table.5[i,"ch"]<-paste(as.numeric(Occ.table.5[i,2:length(Occ.table.5)]),collapse = "")
}

#R MArk Analzses function
lynx.Mark.models=function(CaptHist.table)
{
  #
  # Define parameter models
  #
  pdotshared=list(formula=~1,share=TRUE)
  ptimeshared=list(formula=~time,share=TRUE)
  pmixture=list(formula=~mixture)
  #
  # Run assortment of models
  #
  #
  # Capture Closed models
  #
  # M0: constant p=c
  Lynx.m0=RMark::mark(CaptHist.table,model="Closed",
                    model.parameters=list(p=pdotshared))

  # Mt: Variation accross occasions. time varying p=c
  Lynx.mt=RMark::mark(CaptHist.table,model="Closed",
                    model.parameters=list(p=ptimeshared))
  # Mh: Variation among individuals
  Lynx.Mh=RMark::mark(CaptHist.table,model="HetClosed",
                     model.parameters=list(p=pmixture))

  # Return model table and list of models
  #
  return(RMark::collect.models() )
}


# fit models in mark by calling function created above
# warnings ok

#suppressWarnings(Lynx.Mark.results.3<-lynx.Mark.models(CaptHist.table.3))
#suppressWarnings(Lynx.Mark.results.5<-lynx.Mark.models(CaptHist.table.5))


#Select best model:
# Create Mark Table:
Mark.Table<-data.table(Period=rep(c(3,5),each=4),
                       Model=c("M0","Mb","Mt","Mh"),
                       AIC=0,
                       N=0,
                       Min=0,
                       Max=0)


# Add AICc from RMARK to MArk.Table 
#(to get same value as in Mark programm, unadjusted for Mt and Mh, ex. Lynx.Mark.results.3$Lynx.Mh$results$AICc.unadjusted)
#Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="M0","AICc"]<-round(Lynx.Mark.results.3$Lynx.m0$results$AICc,2)
#Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="Mt","AICc"]<-round(Lynx.Mark.results.3$Lynx.mt$results$AICc,2)
#Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="Mh","AICc"]<-round(Lynx.Mark.results.3$Lynx.Mh$results$AICc,2)

#Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="M0","AICc"]<-round(Lynx.Mark.results.5$Lynx.m0$results$AICc,2)
#Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="Mt","AICc"]<-round(Lynx.Mark.results.5$Lynx.mt$results$AICc,2)
#Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="Mh","AICc"]<-round(Lynx.Mark.results.5$Lynx.Mh$results$AICc,2)

# Estimation N lynx
closedN.3<-secr::closedN(Lynx.CHxy.3, estimator = c("null","zippin","darroch","betabinomial"),
                       level = 0.95, maxN = 1e+07, dmax = 10 )

closedN.5<-secr::closedN(Lynx.CHxy.5, estimator = c("null","zippin","darroch","betabinomial"),
                       level = 0.95, maxN = 1e+07, dmax = 10 )



# Add info in Mark.Table Period = 3
Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="M0","AIC"]<-round(closedN.3$AIC[1],2)
Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="Mb","AIC"]<-round(closedN.3$AIC[2],2)
Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="Mt","AIC"]<-round(closedN.3$AIC[3],2)
Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="Mh","AIC"]<-round(closedN.3$AIC[4],2)

Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="M0","N"]<-round(closedN.3$Nhat[1],0)
Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="Mb","N"]<-round(closedN.3$Nhat[2],0)
Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="Mt","N"]<-round(closedN.3$Nhat[3],0)
Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="Mh","N"]<-round(closedN.3$Nhat[4],0)

Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="M0","SE"]<-round(closedN.3$seNhat[1],2)
Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="Mb","SE"]<-round(closedN.3$seNhat[2],2)
Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="Mt","SE"]<-round(closedN.3$seNhat[3],2)
Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="Mh","SE"]<-round(closedN.3$seNhat[4],2)

Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="M0","Min"]<-round(closedN.3$lclNhat[1],0)
Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="Mb","Min"]<-round(closedN.3$lclNhat[2],0)
Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="Mt","Min"]<-round(closedN.3$lclNhat[3],0)
Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="Mh","Min"]<-round(closedN.3$lclNhat[4],0)

Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="M0","Max"]<-round(closedN.3$uclNhat[1],0)
Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="Mb","Max"]<-round(closedN.3$uclNhat[2],0)
Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="Mt","Max"]<-round(closedN.3$uclNhat[3],0)
Mark.Table[Mark.Table$Period==3 & Mark.Table$Model=="Mh","Max"]<-round(closedN.3$uclNhat[4],0)

# Add info in Mark.Table Period = 5
Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="M0","AIC"]<-round(closedN.5$AIC[1],0)
Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="Mb","AIC"]<-round(closedN.5$AIC[2],0)
Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="Mt","AIC"]<-round(closedN.5$AIC[3],0)
Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="Mh","AIC"]<-round(closedN.5$AIC[4],0)

Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="M0","N"]<-round(closedN.5$Nhat[1],0)
Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="Mb","N"]<-round(closedN.5$Nhat[2],0)
Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="Mt","N"]<-round(closedN.5$Nhat[3],0)
Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="Mh","N"]<-round(closedN.5$Nhat[4],0)

Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="M0","SE"]<-round(closedN.5$seNhat[1],2)
Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="Mb","SE"]<-round(closedN.5$seNhat[2],2)
Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="Mt","SE"]<-round(closedN.5$seNhat[3],2)
Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="Mh","SE"]<-round(closedN.5$seNhat[4],2)

Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="M0","Min"]<-round(closedN.5$lclNhat[1],0)
Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="Mb","Min"]<-round(closedN.5$lclNhat[2],0)
Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="Mt","Min"]<-round(closedN.5$lclNhat[3],0)
Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="Mh","Min"]<-round(closedN.5$lclNhat[4],0)

Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="M0","Max"]<-round(closedN.5$uclNhat[1],0)
Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="Mb","Max"]<-round(closedN.5$uclNhat[2],0)
Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="Mt","Max"]<-round(closedN.5$uclNhat[3],0)
Mark.Table[Mark.Table$Period==5 & Mark.Table$Model=="Mh","Max"]<-round(closedN.5$uclNhat[4],0)

# IC in one column:
Mark.Table$IC<-paste(Mark.Table$Min,Mark.Table$Max,sep="-")
Mark.Table$Min<-NULL
Mark.Table$Max<-NULL

#Order by AIC and period
Mark.Table.3<-Mark.Table[Mark.Table$Period==3,]
Mark.Table.3<-Mark.Table.3[order(Mark.Table.3$AIC),]

Mark.Table.5<-Mark.Table[Mark.Table$Period==5,]
Mark.Table.5<-Mark.Table.5[order(Mark.Table.5$AIC),]

Mark.Table<-rbind(Mark.Table.3,Mark.Table.5)

#Add An external model result to Mark.Table
if(length(External.Model)>1){
  Mark.Table.exct<-data.table(Period=External.Model[1],
                       Model=External.Model[2],
                       AIC="Ext. Model",
                       N=External.Model[3],
                       SE=External.Model[4],
                       IC=External.Model[5])


  Mark.Table<-rbind(Mark.Table,Mark.Table.exct)

  Mark.Table$N<-as.numeric(Mark.Table$N)
  Mark.Table$SE<-as.numeric(Mark.Table$SE)
}


#Make a copy to have shorter code in inline code (for abstract):
Mark.Table.inline<-Mark.Table

#Mark.Table as a felxtable
Mark.Table<-flextable(Mark.Table)
Mark.Table <- set_caption(Mark.Table, caption = " Selection of the best model and estimation of the number of lynx.")
Mark.Table<-font(Mark.Table,fontname = "Calibri", part = "all")
Mark.Table<-fontsize(Mark.Table,size = 12, part = "body")


#In Bold the Best model (lower AIC in period category)
which.best.model.3<-which(Mark.Table$body$dataset$AIC== min(Mark.Table$body$dataset[Mark.Table$body$dataset$Period==3,"AIC"]))
which.best.model.5<-which(Mark.Table$body$dataset$AIC== min(Mark.Table$body$dataset[Mark.Table$body$dataset$Period==5,"AIC"]))

Mark.Table <- bold(Mark.Table, bold = TRUE, i=c(which.best.model.3,which.best.model.5))


#Change name to fit language
Mark.Table <- set_header_labels(Mark.Table,
                        values = list(Period = "Period",
                                      Model = "Model",
                                      AICc = "AIC",
                                      N = "N lynx",
                                      IC = "IC 95%"))

#Create footnote:

if(Language=="FR"){
 Mark.Table<- footnote( Mark.Table, i =1 ,j=2,part = "header",
            value = as_paragraph("M0: P. are constant
Mb: P. vary by behavioural response to capture
Mt: P. vary with time
Mh: P. vary by individual animal"),ref_symbols = c(""))

}

#Autofit (last to make sur it works)
Mark.Table<-autofit(Mark.Table)
Mark.Table <- bold(Mark.Table, bold = TRUE, part = "header")
Mark.Table<-height_all(Mark.Table, height = .26)
Mark.Table<- hrule(Mark.Table, rule = "exact")
Mark.Table<- hline(Mark.Table, i=4)
Mark.Table<-fontsize(Mark.Table, size = 8, part = "footer")

if(length(External.Model)>1){
  Mark.Table<- hline(Mark.Table, i=8)
}
# Import Lynx Observation Data
Data.opp<-data.table::fread(LynxObs,select=(c("lynx_ID","date","type","x","y","canton","SCALP","fk_case_id")))
Data.opp$date<-as.POSIXct(Data.opp$date,  format= '%Y-%m-%d' )

#Subset data Session window
Data.opp<-Data.opp[Data.opp$date>Start & Data.opp$date<Stop,]

#Observations as spatial points
pts <- sf::st_as_sf(Data.opp, coords = c("x","y"))
pts <-sf::st_set_crs(pts, 21781) 

#Import shapefile reference area
Refareashp <- sf::read_sf(dsn = "MAP_Data/Reference_areas_all_new2022.shp")

#Subset observation within the reference area
#Warnings ok
suppressWarnings(pts<-sf::st_intersection(pts,Refareashp$geometry[Refarea]))


#Sort data Opp. Mon.
pts<-as.data.table(pts)
pts<-pts[,c("fk_case_id","date","type","lynx_ID","canton","SCALP")]

Opp.table<-data.table(ID=sort(unique(pts$lynx_ID)))

if(length(Opp.table$ID)>0){
  for(i in 1:length(Opp.table$ID)){
    Opp.table[i,"Event"]<-length(unique(pts[pts$lynx_ID==Opp.table[i,ID],fk_case_id]))
    Opp.table[i,"Min.date"]<-min(pts[pts$lynx_ID==Opp.table[i,ID],date])
    Opp.table[i,"Max.date"]<-max(pts[pts$lynx_ID==Opp.table[i,ID],date])
    Opp.table[i,"Canton"]<-paste(unique(pts[pts$lynx_ID==Opp.table[i,ID],canton]),collapse=", ")
    Opp.table[i,"SCALP"]<-paste(unique(pts[pts$lynx_ID==Opp.table[i,ID],SCALP]),collapse=", ")
  }
}else{warning("No opp. mon. observations")}

Opp.table$Max.date<-as.character(Opp.table$Max.date)
Opp.table$Min.date<-as.character(Opp.table$Min.date)

#Any lynx not seen in det. mon. but seen in opp. mon:
Opp.table$line<-1:length(Opp.table$ID)
Not.seen<-Opp.table[!grepl(paste(unique(Data[Data$animal_species=="Lynx lynx" & Data$TIME>Start & Data$TIME<Stop,id_individual]), collapse="|"),
      Opp.table$ID),which = FALSE]
Opp.table$line<-NULL

#Any lynx not seen in opp. mon but seen in det. mon

Not.seen.opp<-unique(Data[Data$animal_species=="Lynx lynx" & 
              Data$TIME>Start & 
              Data$TIME<Stop &
              !(Data$ID %in% Opp.table$ID),id_individual]) 

#Change name to fit language
names(Opp.table)<-c("ID","Event","Date min.","Date max.","Canton","SCALP")

#Transform data.table into flextable with legend and font
Opp.table<-flextable(Opp.table)
Opp.table<-autofit(Opp.table)
Opp.table <- set_caption(Opp.table, caption = " Opportunistic observations that took place during the monitoring period and were transmitted to KORA. An event can be represented by several images taken at the same place and time.")
Opp.table<-font(Opp.table,fontname = "Calibri", part = "all")
Opp.table<-fontsize(Opp.table,size = 12, part = "body")
Opp.table<-height_all(Opp.table, height = .26)
Opp.table<- hrule(Opp.table, rule = "exact")
Opp.table <- bold(Opp.table, bold = TRUE, i=Not.seen[,line])
Opp.table <- bold(Opp.table, bold = TRUE, part = "header")

#Create footnote:
if(length(Not.seen$line)>0){
  Opp.table<- footnote( Opp.table, i = Not.seen[,line],j=1,
            value = as_paragraph("Lynx observed opportunistically but not observed in the deterministic session are in bold."),ref_symbols = c("1"))
}

if(length(Not.seen$line)==0){
  Opp.table<- footnote( Opp.table, i = 1,j=1,
            value = as_paragraph("All opportunistically observed lynxes were also observed in the deterministic session."),ref_symbols = c("1"),part = "header")
}
Opp.table<-fontsize(Opp.table, size = 8, part = "footer")
# ---- Compute Density 

# - Needed values:
N<-Mark.Table.inline[which.best.model.selected,N]
SE<-Mark.Table.inline[which.best.model.selected,SE]
Surface<-as.numeric(round(sf::st_area(Refareashp[Refarea,])/1000000,0))

# Density per 100 km2
Density.surface<-100*(N/Surface)
Density.surface.min95<-round(Density.surface-((SE/Surface)*100)*1.962,2)
Density.surface.max95<-round(Density.surface+((SE/Surface)*100)*1.962,2)
Density.surface<-round(Density.surface,2)

# Density per 100 km2 suitable habitat
Density.surface.favorable<-100*(N/favorable)
Density.surface.favorable.min95<-round(Density.surface.favorable-((SE/favorable)*100)*1.962,2)
Density.surface.favorable.max95<-round(Density.surface.favorable+((SE/favorable)*100)*1.962,2)
Density.surface.favorable<-round(Density.surface.favorable,2)
#rounding function
foo <- function(x)
{
    round(x+5,-1)
}


jpeg(filename = "Plot_cumulative.jpeg", width = 16.5, height = 8,
     units = "cm", quality = 100, bg = "white", res=1000)

par(mfrow=c(1,2))


par(mar=c(4,3,1,0))
  plot(Occ.table.plot.3$SUM.ID.cumulative, type="o", 
     xlab=NA,ylab=NA,col="blue",
     xaxt='n',las=1,cex.lab=0.8,cex.axis=0.8,pch=19,cex=1.5,lwd=1.5,
     xlim=c(1,20),
     #ylim based on period 5 in case higher:
     ylim=c(0,foo(max(Occ.table.plot.5$SUM.occasion.cumulative))+10)) 
  axis(side=1,at=c(1,5,10,15,20),labels = c(1,5,10,15,20),cex.axis=0.8)
  mtext(side=1,"Occasions", line= 2,cex=0.8)
  mtext(side=2,"Captures", line= 2,cex=0.8)
  par(new=TRUE)
  plot(Occ.table.plot.3$SUM.occasion.cumulative, type="o", 
     xlab = NA,ylab=NA,pch=17,
     axes=FALSE,lwd=1.5,cex=1.5,
     xlim=c(1,20),ylim=c(0,foo(max(Occ.table.plot.5$SUM.occasion.cumulative))+10))
  legend("topleft",paste("Closure Test: Otis et al. (1978), p=",round(closure.result.3$Otis$p,3),
"
Occasion = 3 jours
"),cex=0.6)

par(mar=c(4,0.5,1,0.5))
  plot(Occ.table.plot.5$SUM.ID.cumulative, type="o", 
     xlab=NA,ylab="Captures",col="blue",
     xaxt='n',las=1,cex.lab=0.8,cex.axis=0.8,pch=19,cex=1.5,lwd=1.5,
     axes=FALSE,
     xlim=c(1,12),ylim=c(0,foo(max(Occ.table.plot.5$SUM.occasion.cumulative))+10))
  box()
  axis(side=1,at=c(1:12),labels = c(1:12),cex.axis=0.8)
  mtext(side=1,"Occasions", line= 2,cex=0.8)
  par(new=TRUE)
  plot(Occ.table.plot.5$SUM.occasion.cumulative, type="o", 
     xlab = NA,ylab=NA,pch=17,
     axes=FALSE,lwd=1.5,cex=1.5,
     xlim=c(1,12),ylim=c(0,foo(max(Occ.table.plot.5$SUM.occasion.cumulative))+10))
  legend("topleft",paste("Closure Test: Otis et al. (1978), p=",round(closure.result.5$Otis$p,3),
"
Occasion = 5 jours
"),cex=0.6)


invisible(dev.off())


### Select model Plot in Report

if(Period.selected==3){

  jpeg(filename = "Plot_cumulative_report.jpeg", width = 16.5, height = 8.5,
     units = "cm", quality = 100, bg = "white", res=1000)

  par(mar=c(4,2.7,0.5,0.5),mgp = c(0, 0.2, 0))
  plot(Occ.table.plot.3$SUM.ID.cumulative, type="o", 
     xlab=NA,ylab=NA,col="blue",
     xaxt='n',las=1,cex.lab=0.5,cex.axis=1,pch=19,cex=1.5,lwd=1.5,
     tck=-0.01,
     xlim=c(1,20),
     #ylim based on period 5 in case higher:
     ylim=c(0,foo(max(Occ.table.plot.5$SUM.occasion.cumulative))+20)) 

  axis(side=1,at=c(1:20),labels = c(1:20),cex.axis=0.85,tck=-0.01)
  par(new=TRUE)
  plot(Occ.table.plot.3$SUM.occasion.cumulative, type="o", 
     xlab = NA,ylab=NA,pch=17,
     axes=FALSE,lwd=1.5,cex=1.5,
     xlim=c(1,20),ylim=c(0,foo(max(Occ.table.plot.5$SUM.occasion.cumulative))+20))
  grid(NA,NULL,lwd=1)

  if(Language=="FR"){
     mtext(side=1,"Occasions", line= 2.4,cex=1)
     mtext(side=2,"Captures", line= 1.8,cex=1)
     legend("topleft", legend=c("Captures cumulées", "Nombre de lynx capturés"),
       col=c("black", "blue"),pch=c(17,19), cex=0.85) 
  }

  if(Language=="DE"){
     mtext(side=1,"Fanggelegenheiten", line= 2.4,cex=1)
     mtext(side=2,"Erfassungen", line= 1.8,cex=1)
     legend("topleft", legend=c("Kumulierte Erfassungen", "Anzahl verschiedene Luchse"),
       col=c("black", "blue"),pch=c(17,19), cex=0.85) 
  }

  invisible(dev.off())
}

if(Period.selected==5){

  jpeg(filename = "Plot_cumulative_report.jpeg", width = 14, height = 8.5,
     units = "cm", quality = 100, bg = "white", res=1000)

  par(mar=c(4,2.7,0.5,0.5),mgp = c(0, 0.2, 0))
  plot(Occ.table.plot.5$SUM.ID.cumulative, type="o", 
     xlab=NA,ylab=NA,col="blue",
     xaxt='n',las=1,cex.lab=0.5,cex.axis=1,pch=19,cex=1.5,lwd=1.5,
     tck=-0.01,
     xlim=c(1,12),
     #ylim based on period 5 in case higher:
     ylim=c(0,foo(max(Occ.table.plot.5$SUM.occasion.cumulative))+10)) 
  axis(side=1,at=c(1:12),labels = c(1:12),cex.axis=1,tck=-0.01)
  par(new=TRUE)
  plot(Occ.table.plot.5$SUM.occasion.cumulative, type="o", 
     xlab = NA,ylab=NA,pch=17,
     axes=FALSE,lwd=1.5,cex=1.5,
     xlim=c(1,12),ylim=c(0,foo(max(Occ.table.plot.5$SUM.occasion.cumulative))+10))
  grid(NA,NULL,lwd=1)

  if(Language=="FR"){
     mtext(side=1,"Occasions", line= 2.4,cex=1)
     mtext(side=2,"Captures", line= 1.8,cex=1)
     legend("topleft", legend=c("Captures cumulées", "Nombre de lynx capturés"),
       col=c("black", "blue"),pch=c(17,19), cex=0.85) 
  }

  if(Language=="DE"){
     mtext(side=1,"Fanggelegenheiten", line= 2.4,cex=1)
     mtext(side=2,"Erfassungen", line= 1.8,cex=1)
     legend("topleft", legend=c("Kumulierte Erfassungen", "Anzahl verschiedene Luchse"),
       col=c("black", "blue"),pch=c(17,19), cex=0.85) 
  }



  invisible(dev.off())
}
History<-as.matrix(fread("Density_KORA_Historic_Values.csv",header = TRUE))

Komp.DE<-c("Nordostschweiz","Misox-Südtessin","Tessin","Surselva","Zentralschweiz Ost","Rhone-Nord","Unterwallis Süd","Simme-Saane","Engadin","Jura Nord","Jura Süd","Mittelbünden", "Zentralschweiz West","Zentralschweiz Mitte","Berner Oberland Ost","Oberwallis")

Komp.FR<-c("Nord-est de la Suisse","Val Mesolcina-Sud du Tessin","Tessin","Surselva","Est de la Suisse centrale","Nord du Rhône","Sud du Bas-Valais", "Simme-Saane", "Engadine", "Nord du Jura","Sud du Jura", "Centre des Grisons","Ouest de la Suisse centrale", "Centre de la Suisse centrale","Est de l'Oberland-Bernois","Haut-Valais")

if(Language=="FR"){
  History.comparaison<-data.table(Compartment=Komp.FR)
}
if(Language=="DE"){
  History.comparaison<-data.table(Compartment=Komp.DE)
}

History.comparaison$COMPID<-c("II","Vb","Va","Vc","IIIc","IVc","IVd","IVa","Ve","Ib","Ia","Vd","IIIa","IIIb","IVb","IVe")

for(i in 1:length(History.comparaison$Compartment)){

  list.comp<-tail(as.character(na.omit(History[,i+1])),n=1)
  if(length(list.comp)>0){
  History.comparaison[i,"Winter"]<-History[which(History == list.comp, arr.ind=T)[,"row"],1]
  History.comparaison[i,"Density"]<-as.numeric(substr(list.comp, 1, 4))
  History.comparaison[i,"IC95"]<-paste("(",as.numeric(substr(list.comp, 7, 10)),"-",as.numeric(substr(list.comp, 12, 15)),")",sep="")

  }else{
    if(Language=="FR"){History.comparaison[i,"Winter"]<-"-"}
    if(Language=="DE"){History.comparaison[i,"Winter"]<-"-"}
    }

}

#Compute current values

Current<-data.table(Year=paste(juv.year,max(year(occasions)),sep="/"),
                    Density=Density.surface.favorable, 
                    Min95=Density.surface.favorable.min95,
                    Max95=Density.surface.favorable.max95)


#Add new value of this session
History.comparaison<-as.data.frame(History.comparaison)
History.comparaison[Compartment,"Winter"]<-Current[,Year]
History.comparaison[Compartment,"Density"]<-Current[,Density]
History.comparaison[Compartment,"IC95"]<-paste("(",Current[,Min95],"-",Current[,Max95],")", sep="")
History.comparaison<-as.data.table(History.comparaison)

#Set order COMPID (alphabetic order does not work->manual)
History.comparaison<-History.comparaison[c(11,10,1,13,14,5,8,15,6,7,16,3,2,4,12,9),]

#Set order Density
History.comparaison<-History.comparaison[order(History.comparaison$Density,decreasing=TRUE,na.last=TRUE),]

#Get position actual session:
line.session<-History.comparaison[COMPID==names(Text.table)[Refarea+1],which=TRUE]
#Get position for (NA-NA) sessions:
line.NA<-History.comparaison[COMPID=="IVe"|COMPID=="IVd",which=TRUE]

#change "." to ","
for(i in 1:nrow(History.comparaison)){
  History.comparaison[i,"Density_text"]<-prettyNum(History.comparaison[i,Density],decimal.mark = ",")
  History.comparaison[i,"IC95"]<-gsub("[.]", ",", History.comparaison[i,IC95])
}

#reset order column:
History.comparaison$Density<-NULL
names(History.comparaison)[5]<-"Density"
History.comparaison<-History.comparaison[,c("COMPID","Compartment","Winter","Density","IC95")]

#Change NA to "-" in density column:
History.comparaison[History.comparaison$Density=="NA","Density"]<-"-"

#Set name columns (FR/DE)
if(Language=="FR"){
  names(History.comparaison)<-c("ID Comp.","Aire de référence","Hiver","Densité","IC 95% ")
}
if(Language=="DE"){
  names(History.comparaison)<-c("ID Komp.","Referenzgebiet","Winter","Dichte","95% KI")
}

#Transform data.table into flextable with legend and font
History.comparaison<-flextable(History.comparaison)
History.comparaison<-autofit(History.comparaison)

if(Language=="FR"){
  History.comparaison <- set_caption(History.comparaison, caption = "Pour chacune des 16 aires de référence officielles figure la dernière estimation de la densité de lynx (lynx indépendants pour 100 km^2^ d’habitat favorable) avec l'intervalle de confiance de 95%, à l'exception de celles où il n'y a pas encore eu de session. Les densités sont classées par ordre décroissant. En gras\ua0: les valeurs de la session traitée dans le présent rapport. Les sous-compartiments sont présentés à la Fig. 1.",autonum = officer::run_autonum(
  seq_id = "table",
  pre_label = "Tab. ",
  post_label = ". ",
  bkm = NULL,
  bkm_all = FALSE,
  prop = officer::fp_text(bold= TRUE, font.size = 12, font.family= "Calibri") ,
  start_at = ""
))
}

if(Language=="DE"){
  History.comparaison <- set_caption(History.comparaison, caption = "Für jedes der 16 offiziellen Referenzgebiete wird die letzte Schätzung der Luchsdichte (selbständige Luchse pro 100 km^2^ geeignetem Habitat) mit dem 95 %-Konfidenzintervall angegeben, ausser für die Gebiete, in denen noch kein Durchgang stattgefunden hat. Die Dichten sind in absteigender Reihenfolge aufgeführt. **Fett**: Werte des in diesem Bericht behandelten Durchgangs. Die Teil-Kompartimente sind in der Abb. 1 kartografisch dargestellt.",autonum = officer::run_autonum(
  seq_id = "table",
  pre_label = "Tab. ",
  post_label = ". ",
  bkm = NULL,
  bkm_all = FALSE,
  prop = officer::fp_text(bold= TRUE, font.size = 12, font.family= "Calibri") ,
  start_at = ""
))
}

#Create footnote:
if(Language=="FR"){
  History.comparaison<- footnote( History.comparaison, i = line.NA,j=5,
            value = as_paragraph("NA figure lorsque l’abondance et par la même la précision n’ont pas pu être estimées au moyen de la méthode de capture-recapture photographique, car le nombre de lynx est trop faible."),ref_symbols = c("*"),part = "body")
}

if(Language=="DE"){
  History.comparaison<- footnote( History.comparaison, i = line.NA,j=5,
            value = as_paragraph("NA: Abundanz und damit das Konfidenzintervall konnten aufgrund der zu geringen Anzahl Luchse nicht mit der Fang-Wiederfang Methode geschätzt werden."),ref_symbols = c("*"),part = "body")
}

History.comparaison<-bold(History.comparaison, i = line.session, j = NULL, bold = TRUE, part = "body")
History.comparaison <- bold(History.comparaison, bold = TRUE, part = "header")

History.comparaison<-font(History.comparaison,fontname = "Calibri", part = "all")
History.comparaison<-fontsize(History.comparaison,size = 12, part = "body")
History.comparaison<-height_all(History.comparaison, height = .27)
History.comparaison<- hrule(History.comparaison, rule = "exact")
History.comparaison<-fontsize(History.comparaison, size = 8, part = "footer")
#History Table

History<-fread("Density_KORA_Historic_Values.csv",header = TRUE,select=c("Year",as.character(Compartment)))

History<-na.omit(History,cols = as.character(Compartment))

#Split historic values into columns:
#Warnings are produced if first session ever in area (normal)
suppressWarnings(
for(i in 1:length(History$Year)){
  History[i,"Density"]<-as.numeric(substr(as.character((History[i,2])), 1, 4))
  History[i,"Min95"]<-as.numeric(substr(as.character((History[i,2])), 7, 10))
  History[i,"Max95"]<-as.numeric(substr(as.character((History[i,2])), 12, 15))
})

# Add current Monitoring values to history plot:
#Object current is created in Density comparaison

History[,2]<-NULL

History<-rbind(History,Current)

# Error bar
History.plot<-ggplot2::ggplot(History) +
              ggplot2::geom_bar( ggplot2::aes(x=Year, y=Density), stat="identity",
                                 fill=c(rep("#009E73",length(History$Year)-1),"#E69F00"),
                                 alpha=1) +
              ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45,hjust = 1),
                             axis.title.x=ggplot2::element_blank())+
              ggplot2::geom_errorbar( ggplot2::aes(x=Year, ymin=Min95, ymax=Max95),
                                      width=0.08, colour="black", alpha=0.9, size=0.5)+ggplot2::theme_bw()+
              ggplot2::geom_text(ggplot2::aes(y = 0.4,x=1:length(History$Year)),
                                              label=prettyNum(History$Density,decimal.mark = ","),
                                              vjust = 0, nudge_y = 0,size=3)+
              ggplot2::geom_text(ggplot2::aes(y = 0,x=1:length(History$Year)),
                                              label=paste("(",prettyNum(History$Min95,decimal.mark = ",")
                                              ,"-",
                                              prettyNum(History$Max95,decimal.mark = ","),")",sep=""),
                                              vjust = 0, nudge_y = .2,size=3)+
              ggplot2::scale_y_continuous(breaks = seq(0,6,by=0.5),
                                          minor_breaks=NULL) 


#Add label
if(Language=="FR"){
History.plot<-History.plot+ggplot2::xlab("Hiver")+ggplot2::ylab("Densité")
}

if(Language=="DE"){
History.plot<-History.plot+ggplot2::xlab("Winter")+ggplot2::ylab("Dichte")
}


ggplot2::ggsave("History_denisty_plot.jpeg",plot=History.plot,
                    units = "cm",
                    width = 16.5,
                    height = 8.5)
#Import shapefile compartment
Komp <- sf::read_sf(dsn = "MAP_Data/Luchs_Komp_21_07_2015_new.shp")

# -- Altitude layer #
Alt <- raster::raster("MAP_Data/CHHS.TIF")
# indicate CRS
raster::crs(Alt) <- CRS 
# convert to a df for plotting in two steps:
# First, to a SpatialPointsDataFrame
Alt_pts <- raster::rasterToPoints(Alt, spatial = TRUE)

# -- Use ingfo to compute Study Area 
study_area<-sp::bbox(Alt_pts)
study_area<-rgeos::readWKT(paste("POLYGON((",
                                   study_area[1,1]," ",study_area[2,1],",",
                                   study_area[1,1]," ",study_area[2,2],",",
                                   study_area[1,2]," ",study_area[2,2],",",
                                   study_area[1,2]," ",study_area[2,1],",",
                                   study_area[1,1]," ",study_area[2,1],"))",sep=""),p4s= CRS)

# Then convert SpatialPointsDataFrame to a 'conventional' dataframe
Alt_df  <- data.frame(Alt_pts)
rm(Alt_pts, Alt)

# -- Suitable habitat layer #
suppressWarnings(Suitable <- raster::raster("MAP_Data/suitable_CH1.tif"))
# indicate CRS
raster::crs(Suitable) <- CRS 
# crop 
Suitable<-raster::crop(Suitable,study_area)
# convert to a df for plotting in two steps:
# First, to a SpatialPointsDataFrame
Suitable_pts <- raster::rasterToPoints(Suitable, spatial = TRUE)
# Then to a 'conventional' dataframe
Suitable_df  <- data.frame(Suitable_pts)
rm(Suitable_pts, Suitable)

# -- Import Suitable habitat Liechtenstein

#Import
Suitable_LIE <- rgdal::readOGR("MAP_Data/Shape_model_Alps.shp", verbose=FALSE)
Suitable_LIE <-raster::buffer(Suitable_LIE, width=0) 
Suitable_LIE <- sp::spTransform(Suitable_LIE,CRS)
#Crop FL
suppressWarnings(LIE <- raster::shapefile("MAP_Data/FL.shp"))
LIE <-raster::buffer(LIE, width=0) 
Suitable_LIE<-raster::crop(Suitable_LIE,LIE)
#Set as sf
Suitable_LIE<-sf::st_as_sf(Suitable_LIE)

# -- Import shapefile lakes
suppressWarnings(Lakes <- raster::shapefile("MAP_Data/grandlacs.shp"))
Lakes<-raster::crop(Lakes,study_area)
Lakes<-sf::st_as_sf(Lakes)

#Import shapefile white background
Mask <- suppressWarnings(raster::shapefile("MAP_Data/Mask_CH_Liechtenstein.shp"))
Mask<-raster::crop(Mask,study_area)
Mask<-sf::st_as_sf(Mask)

# --- Import shapefile Study area Polygon
suppressWarnings(Rcompartment <- raster::shapefile("MAP_Data/Reference_areas_all_new2022.shp"))
Rcompartment<-sf::st_as_sf(Rcompartment)

map<-ggplot2::ggplot() +
  #Alt Raster
  ggplot2::geom_raster(data = Alt_df , ggplot2::aes(x = x, y = y, alpha = CHHS))+ 
  ggplot2::scale_alpha(name = "", range = c(0.6, 0))+
  #Background white
  ggplot2::geom_sf(data=Mask,fill="white")+
  #Suitable Raster
  ggplot2::geom_tile(data = Suitable_df , ggplot2::aes(x = x, y = y),fill="yellowgreen", alpha = 0.4)+
  #Suitable Raster LIE
  ggplot2::geom_sf(data=Suitable_LIE$geometry,fill="yellowgreen",colour=NA, alpha = 0.4)+
  #Polygon Compartment all without Rohne Nord
  ggplot2::geom_sf(data=Komp$geometry,colour = "black",fill="NA",lwd=0.6)+
  #Polygon Compartment
  ggplot2::geom_sf(data=Komp$geometry[Compartment],fill="white", alpha = 0.6)+
  #Lakes
  ggplot2::geom_sf(data=Lakes,fill="#56B4E9")+
  #Reference Area all
  ggplot2::geom_sf(data=Rcompartment,col="darkblue",fill="deepskyblue1",lwd=0.5,alpha=0.1)+
  #Reference Area
  ggplot2::geom_sf(data=Rcompartment[Refarea,],col="darkblue",fill="blue",lwd=1.3,alpha=0.2)+
  #Polygon Compartment all text
  ggplot2::geom_sf_text(data=Komp,ggplot2::aes(label =  Nummer))




  #Theme
  map<-map+ggplot2::theme(plot.title =ggplot2::element_blank(),
        strip.background = ggplot2::element_blank(),
        panel.background = ggplot2::element_blank(),
        panel.grid = ggplot2::element_blank(),
        axis.title =ggplot2::element_blank(),
        axis.text = ggplot2::element_blank(),
        axis.ticks= ggplot2::element_blank(),
        legend.position="none",
        panel.border = ggplot2::element_rect(colour = "black", fill=NA, size=0.1))+
        ggplot2::scale_x_continuous(expand = c(0, 0)) +
        ggplot2::scale_y_continuous(expand = c(0, 0))



  # --- Add KORA GIS Logo:
  img <- png::readPNG("KORAlogo.png")#KoraGis_transp
  g <- grid::rasterGrob(img, interpolate=TRUE)

  map<-map+
  ggplot2::annotation_custom(g, xmin=study_area[1]@bbox[1,2]-40000, xmax=study_area[1]@bbox[1,2],
                               ymin=study_area[1]@bbox[2,1], ymax=study_area[1]@bbox[2,1]+40000)






  ggplot2::ggsave("Plot_Introduction_Compartment.jpeg",plot=map,
                    units = "cm",
                    width = 16.5,
                    height = 10.3)
#Import shapefile compartment
Komp <- sf::read_sf(dsn = "MAP_Data/Luchs_Komp_21_07_2015_new.shp")

# -- Study Area 
study_area<-sp::bbox(as.matrix(Komp$geometry[Compartment][[1]]))
study_area.origin<-sp::bbox(as.matrix(Komp$geometry[Compartment][[1]]))
#Add x% around the BBox to have some extra map area

  study_area[1,1]<-study_area[1,1]-round(study_area[1,1]*0.02) #x min
  study_area[2,1]<-study_area[2,1]-round(study_area[2,1]*6*0.02) #y min
  study_area[1,2]<-study_area[1,2]+round(study_area[1,2]*0.02) #x max
  study_area[2,2]<-study_area[2,2]+round(study_area[2,2]*0.02) #y max

study_area<-rgeos::readWKT(paste("POLYGON((",
                                   study_area[1,1]," ",study_area[2,1],",",
                                   study_area[1,1]," ",study_area[2,2],",",
                                   study_area[1,2]," ",study_area[2,2],",",
                                   study_area[1,2]," ",study_area[2,1],",",
                                   study_area[1,1]," ",study_area[2,1],"))",sep=""),p4s= CRS)

# -- Suitable habitat layer #
suppressWarnings(Suitable <- raster::raster("MAP_Data/suitable_CH1.tif"))
# indicate CRS
raster::crs(Suitable) <- CRS 
# crop 
Suitable<-raster::crop(Suitable,study_area)
# convert to a df for plotting in two steps:
# First, to a SpatialPointsDataFrame
Suitable_pts <- raster::rasterToPoints(Suitable, spatial = TRUE)
# Then to a 'conventional' dataframe
Suitable_df  <- data.frame(Suitable_pts)
rm(Suitable_pts, Suitable)

# --- Import shapefile Study area Polygon
suppressWarnings(Rcompartment <- raster::shapefile("MAP_Data/Reference_areas_all_new2022.shp"))
Rcompartment<-Rcompartment[Refarea,]

#Import shapefile compartment
suppressWarnings(Komp <- raster::shapefile("MAP_Data/Luchs_Komp_21_07_2015_new.shp"))
Komp<-Komp[Compartment,]

# ------ Open Source Map#### !!! Study area comes fromr Plot CH (Introduction)

#study area BBox in lat long
study_area_lat_long <- sp::spTransform(study_area, sp::CRS("+init=epsg:4326"))

  LAT1 = study_area_lat_long@bbox[2,1] ; LAT2 = study_area_lat_long@bbox[2,2]
  LON1 = study_area_lat_long@bbox[1,1] ; LON2 = study_area_lat_long@bbox[1,2]

#Get the map
#warnings OK
suppressWarnings(map <- OpenStreetMap::openmap(c(LAT2,LON1), c(LAT1,LON2), zoom = 10, #can be replaced by NULL
                                                 type = c("stamen-terrain")[1],
                                                 mergeTiles = TRUE))

#Correct projection
#warnings OK
suppressWarnings(map <- OpenStreetMap::openproj(map, projection = "+init=epsg:21781"))

# ------ Build map 

# --- Open Source Map in ggplot:
map<-OpenStreetMap::autoplot.OpenStreetMap(map)+
ggplot2::theme(axis.title=ggplot2::element_blank(),
                   axis.text=ggplot2::element_blank(),
                   axis.ticks=ggplot2::element_blank(),
                   panel.border = ggplot2::element_rect(colour = "white", fill=NA, size=12))



# --- Add other layers:

map<-map+
#Suitable Raster
ggplot2::geom_tile(data = Suitable_df , ggplot2::aes(x = x, y = y),fill="yellowgreen", alpha = 0.6)+
#KOMP
ggplot2::geom_polygon(data = Komp@polygons[[1]], ggplot2::aes(x=long, y=lat,group = group),
                               colour="black", fill="white", alpha=0.4, lwd=1.1)+
#Refarea
ggplot2::geom_polygon(data = Rcompartment@polygons[[1]], ggplot2::aes(x=long, y=lat,group = group),
                              col="darkblue",fill="blue",lwd=1.1,alpha=0.2)+
#Points Sites
ggplot2::geom_point(data=DataFULL,ggplot2::aes(x,y), col="white", pch=19,cex=2)+
      ggplot2::geom_point(data=DataFULL,ggplot2::aes(x,y),col="black", pch=1,cex=2)


# --- Add Scale:

  # distance on x axes:
  dist.scale<-plyr::round_any(round(((study_area@bbox[1,2]-study_area@bbox[1,1])/1000)/4), 10, f = ceiling)

  # scale thickness:
  s.thick<-(0.01*(study_area@bbox[2,2]-study_area@bbox[2,1]))

  # scale black rectangle
  xleft<-study_area@bbox[1,1]+((study_area.origin[1,1]-study_area@bbox[1,1])/3)#left
  xright<-xleft+dist.scale*1000#right
  ybottom<-study_area.origin[2,1]-4*s.thick#bottom
  ytop<-study_area.origin[2,1]-3*s.thick#top

  map<-map+
    ggplot2::geom_rect(mapping=ggplot2::aes(xmin=xleft, xmax=xright, ymin=ybottom, ymax=ytop),
                       fill=c("black"),
                       inherit.aes = FALSE)+
    ggplot2::geom_text(x=xright+(2/5)*dist.scale*1000, y=ytop, label=paste(dist.scale,"Km",sep=" "), cex=4, color ="black")

#Save Plot  

ggplot2::ggsave("Plot_Introduction.jpeg",plot=map,
                    units = "cm",
                    width = 16.5,
                    height = 12)
#check if any lynx of the table (Event.report) was once seen far away (> 150 km)

Dispersal<-data.table(ID=Lynx.Table.Adult$ID)

distfunction<-function(x1,x2,y1,y2){
  ((x1-x2)^2+(y1-y2)^2)^0.5
}

for(i in 1:nrow(Dispersal)){
  Dispersal[i,"Lastx"]<-tail(Data[Data$id_individual==Dispersal[i,ID],],1)[,31]
  Dispersal[i,"Lasty"]<-tail(Data[Data$id_individual==Dispersal[i,ID],],1)[,32]

  #Check longest distance in Photo Predator:
  loop1<-PhotoPredator.table[PhotoPredator.table$id_individual==Dispersal[i,ID],]
  for(j in 1:nrow(loop1)){
    loop1[j,"Dist"]<-distfunction(Dispersal[i,Lastx],loop1[j,x],Dispersal[i,Lasty],loop1[j,y])
  }
  Dispersal[i,"PP"]<-round(max(loop1$Dist),0)
  Dispersal[i,"PP.x"]<-unique(loop1[loop1$Dist==max(loop1$Dist),x])
  Dispersal[i,"PP.y"]<-unique(loop1[loop1$Dist==max(loop1$Dist),y])


  #Check longest distance in Lynx Obs:
  loop2<-LynxObs.table[LynxObs.table$lynx_ID==Dispersal[i,ID],]

  if(length(loop2$lynx_ID)>0){

  for(k in 1:nrow(loop2)){
  loop2[k,"Dist"]<-distfunction(Dispersal[i,Lastx],loop2[k,x],Dispersal[i,Lasty],loop2[k,y])
  }
  Dispersal[i,"LO"]<-round(max(loop2$Dist),0)
  Dispersal[i,"LO.x"]<-unique(loop2[loop2$Dist==max(loop2$Dist),x])
  Dispersal[i,"LO.y"]<-unique(loop2[loop2$Dist==max(loop2$Dist),y])

  }

  if(length(loop2$lynx_ID)>0){
  Dispersal[i,"MaxDist"]<-round(max(max(loop1$Dist),max(loop2$Dist)),0)
  }else{Dispersal[i,"MaxDist"]<-round(max(loop1$Dist),0)}


}

Dispersal<-Dispersal[order(Dispersal$MaxDist,decreasing=TRUE),]
Dispersal<-Dispersal[,c("ID","Lastx","Lasty","PP","PP.x","PP.y","LO","LO.x","LO.y","MaxDist")]

#Dispersal as a flextable:
Dispersal<-flextable(Dispersal)
Dispersal <- set_caption(Dispersal, caption = "Dispersal table. This table give an indication of the longest distance between the last observation of a given lynx in the session and the entire KORA data base (Photo Predator + Lynx observations). If the distance is big, please check if it is a dispersal case. PP= Photo Predator's table, LO = Lynx Observations' table. Decreasing order from column MaxDist. [Unit = meters]")
Dispersal<-font(Dispersal,fontname = "Calibri", part = "all")
Dispersal <- bold(Dispersal, bold = TRUE, part = "header")
Dispersal<-autofit(Dispersal)










cat(paste("#",Title,sep=" "))


cat(paste("##",Authors, sep=" "))

\newpage

#Prepare sentence difference CI:
if(length(History$Year)==1){
   DifferentFR<-paste("**Il s'agit de la première estimation de densité pour cette aire de référence.**")
   DifferentDE<-paste("**Dies ist die erste Dichteschätzung für dieses Referenzgebiet.**")
   DifferentEN<-paste("**This is the first density estimate for this reference area.**")

  }else{
    smaller = (History[nrow(History)-1,Min95]>History[nrow(History),Max95])
    bigger = (History[nrow(History)-1,Max95]<History[nrow(History),Min95])

    if(smaller == FALSE & bigger == FALSE){
         DifferentFR<-paste("**La densité estimée n’a pas changé significativement par rapport à la valeur estimée lors de la dernière session.**")
         DifferentDE<-paste("**Die geschätzte Dichte hat sich im Vergleich, zu dem im letzten Durchgang geschätzten Wert nicht signifikant verändert.**")
         DifferentEN<-paste("**The estimated density has not changed significantly from the value estimated in the last session.**")
    }

    if(smaller == TRUE){
         DifferentFR<-paste("**La densité estimée est significativement plus basse par rapport à la valeur estimée lors de la dernière session.**")
         DifferentDE<-paste("**Die geschätzte Dichte ist signifikant niedriger als der im letzten Durchgang geschätzte Wert.**")
         DifferentEN<-paste("**The estimated density is significantly lower than the value estimated in the last session.**")
    }

    if(bigger == TRUE){
         DifferentFR<-paste("**La densité estimée est significativement plus haute par rapport à la valeur estimée lors de la dernière session.**")
         DifferentDE<-paste("**Die geschätzte Dichte ist signifikant höher als der im letzten Durchgang geschätzte Wert.**")
         DifferentEN<-paste("**The estimated density is significantly higher than the value estimated in the last session.**")
    }
  }

#Prepare Winter Name
if(month(min(Data$TIME))<5){
  Wintername<-paste(year(min(Data$TIME))-1, substr(year(max(Data$TIME)),3,4),sep='/')
}
if(month(min(Data$TIME))>5){
  Wintername<-paste(year(min(Data$TIME)), substr(year(max(Data$TIME)),3,4),sep='/')
}



#Build abstract (in 3 languages):
AbstractDE<-capture.output(cat("Das  Fotofallen-Monitoring  des  Luchses  (*Lynx lynx*) im  Referenzgebiet ", Text.table[1,Refarea+1],"  wurde im Winter ", Wintername,"  während  60  Nächten, vom ", format(as.Date(Start), "%d.%m.%Y")," bis ", format(as.Date(Stop), "%d.%m.%Y")," durchgeführt. Die Fotofallen an den ", Sites," Standorten funktionierten während ", prettyNum(effective.TN, big.mark="'",decimal.mark=",")," der potentiellen ", prettyNum(potentially.TN, big.mark="'",decimal.mark=",")," Fallennächte (", gsub("[.]",",",sprintf("%.1f",round((effective.TN*100)/potentially.TN,1))),"%). Im Durchgang wurden während ", Event[ID=="Total",Event]," Ereignissen ", N.ind.lynx," selbständige Luchse an ", Event[ID=="Total",`N site`]," Standorten fotografiert. **Darüber hinaus wurden ", N.juv," Jungtier(e) aus mindestens ", length(unique(Lynx.Table[Lynx.Table$YOB==juv.year & Lynx.Table$mother!="",mother]))," Wurf/Würfen nachgewiesen.** Die Fang-Wiederfang Schätzung der Abundanz (95% Konfidenzintervall) nach dem Modell M~", substr(Mark.Table.inline[which.best.model.selected,Model], 2, 3),"~ ergab ", Mark.Table.inline[which.best.model.selected,N]," (",sub("\\-.*", "", Mark.Table.inline[which.best.model.selected,IC])," - ",sub("..-*", "", Mark.Table.inline[which.best.model.selected,IC]),") selbständige Luchse im Referenzgebiet, was einer Dichte von ", prettyNum(Density.surface.favorable,decimal.mark = ",")," (", prettyNum(Density.surface.favorable.min95,decimal.mark = ",")," - ", prettyNum(Density.surface.favorable.max95,decimal.mark = ","),") selbständigen Luchsen pro 100 km^2^ geeignetem Habitat entspricht. ",DifferentDE,sep = ""))

AbstractFR<-capture.output(cat("Le monitoring du lynx (*Lynx lynx*) par piège-photographique dans l'aire de référence ", Text.table[7,Refarea+1]," durant l’hiver ", Wintername," a été effectué durant 60 nuits du ", format(as.Date(Start), "%d.%m.%Y")," au ",format(as.Date(Stop), "%d.%m.%Y"),". Les pièges-photos placés auprès des ", Sites," sites ont fonctionné pendant ", prettyNum(effective.TN, big.mark="'",decimal.mark=",")," des ", prettyNum(potentially.TN, big.mark="'",decimal.mark=",")," nuits potentielles (",  gsub("[.]",",",sprintf("%.1f",round((effective.TN*100)/potentially.TN,1))),"%). Pendant la session, ", Event[ID=="Total",Event]," événements auprès de ", Event[ID=="Total",`N site`]," sites correspondant à ", N.ind.lynx," lynx indépendants ont été répertoriés. **De plus, ", N.juv," juvénile(s) d'au moins ", length(unique(Lynx.Table[Lynx.Table$YOB==juv.year & Lynx.Table$mother!="",mother]))," portée(s) ont également été détectés.** L’estimation de l’abondance (intervalle de confiance de 95%) par le modèle de capture-recapture M~", substr(Mark.Table.inline[which.best.model.selected,Model], 2, 3),"~ était de ", Mark.Table.inline[which.best.model.selected,N]," (",sub("\\-.*", "", Mark.Table.inline[which.best.model.selected,IC])," - ",sub("..-*", "", Mark.Table.inline[which.best.model.selected,IC]),") lynx indépendants ce qui correspond à une densité de ", prettyNum(Density.surface.favorable,decimal.mark = ",")," (", prettyNum(Density.surface.favorable.min95,decimal.mark = ",")," - ", prettyNum(Density.surface.favorable.max95,decimal.mark = ","),") lynx indépendants pour 100 km^2^ d’habitat favorables. ",DifferentFR,sep = ""))

AbstractEN<-capture.output(cat("The monitoring of the Eurasian  lynx  (*Lynx lynx*)  by means of camera traps in the reference area ", Text.table[8,Refarea+1]," during  winter ", Wintername," was  carried  out  during 60  nights,  from  ", format(as.Date(Start), "%d.%m.%Y"),"  to  ", format(as.Date(Stop), "%d.%m.%Y"),". The camera traps at the ", Sites," locations operated during ", prettyNum(effective.TN, big.mark=",",decimal.mark=".")," of the potential ", prettyNum(potentially.TN, big.mark=",",decimal.mark=".")," trap nights (", sprintf("%.1f",round((effective.TN*100)/potentially.TN,1)),"%). During the session, ", Event[ID=="Total",Event]," events of ", N.ind.lynx," independent lynx at ", Event[ID=="Total",`N site`]," sites  were recorded. **In addition, ", N.juv," juvenile(s) of at least ", length(unique(Lynx.Table[Lynx.Table$YOB==juv.year & Lynx.Table$mother!="",mother]))," litter(s) were detected.** The capture-recapture estimate of abundance (95% confidence interval) under model M~",substr(Mark.Table.inline[which.best.model.selected,Model], 2, 3),"~ was ", Mark.Table.inline[which.best.model.selected,N]," (",sub("\\-.*", "", Mark.Table.inline[which.best.model.selected,IC])," - ",sub("..-*", "", Mark.Table.inline[which.best.model.selected,IC]),") independent lynx, which corresponds to a density of  ", Density.surface.favorable," (", Density.surface.favorable.min95," - ", Density.surface.favorable.max95,") independent lynx per 100 km^2^ of suitable habitat. ",DifferentEN,sep = ""))

#Set Order abstract in the report:
if(Language=="FR"){

cat("### Résumé \n")

cat(AbstractFR)
cat("

\n

")

cat("### Zusammenfassung \n")

cat(AbstractDE)
cat("

\n

")

cat("### Abstract \n")

cat(AbstractEN)
cat("

\n

")

}

if(Language=="DE"){

cat("### Zusammenfassung \n")

cat(AbstractDE)
cat("

\n

")

cat("### Résumé \n")

cat(AbstractFR)
cat("

\n

")

cat("### Abstract \n")

cat(AbstractEN)
cat("

\n

")

}

\newpage

cat("### Introduction \n")

cat("Le monitoring déterministe du lynx par pièges-photos a été développé en Suisse à partir de 1998 (Laass 1999) dans le nord-ouest des Alpes. Dans les sous-compartiments de gestion (voir [Plan Lynx Suisse OFEV 2016](https://www.kora.ch/?action=get_file&id=194&resource_link_id=406)) dont l'habitat favorable pour le lynx est en grande partie colonisé de façon permanente par l’espèce, un suivi déterministe par pièges-photos est effectué dans des aires dites de référence. Ces dernières sont choisies de sorte qu’elles soient représentatives pour les sous-compartiments de gestion donnés (Fig. 1\ua0; voir [ici](https://www.kora.ch/fr/projets/monitoring-grands-carnivores/quest-ce-que-le-monitoring) pour plus d’informations). Au sein de chaque aire de référence, des pièges-photos sont répartis de manière systématique et laissés en place pendant une période donnée, dans notre cas généralement pendant 60 nuits. Ce genre de suivi est conduit tous les trois à quatre ans par le KORA avec l’aide des cantons afin d’estimer la taille et la densité de la population de lynx au moyen de la méthode dite de capture-recapture photographique qui est aujourd’hui une méthode de monitoring standard pour le suivi des espèces aux mœurs cryptiques dont les individus sont reconnaissables au moyen de marques naturelles (voir [ici](https://www.kora.ch/fr/projets/monitoring-grands-carnivores/monitoring-deterministe-par-piegeage-photographique) pour plus d’informations).

Ce rapport présente les résultats de la session déterministe dans l'aire de référence ",Text.table[7,Refarea+1]," au sein du sous-compartiment ",names(Text.table)[Refarea+1]," pour l’hiver ", Wintername,". Les résultats des sessions précédentes sont disponibles [en ligne](https://www.kora.ch/fr/especes/lynx/piegeage-photographique) sur le site du KORA.
\n",sep="")

cat("![**Fig. 1.** Répartition des aires de références (polygones bleus) au sein des 16 sous-compartiments de gestion (polygones noirs). En vert clair, l'habitat favorable du lynx défini par un modèle d'habitat (Zimmermann 2004, Zimmermann & Breitenmoser 2007). L'aire de références (",Text.table[7,Refarea+1],") ainsi que le sous-compartiment (",names(Text.table)[Refarea+1],") de gestion dans lesquels s'est déroulée cette session sont mis en évidence.](Plot_Introduction_Compartment.jpeg) 

\n",sep="")
cat("### Einleitung \n")

cat("Das deterministische Luchsmonitoring mit Fotofallen wurde in der Schweiz ab 1998 (Laass 1999) in den Nordwestalpen entwickelt. In Teil-Kompartimenten (siehe [Konzept  Luchs Schweiz BAFU 2016](https://www.kora.ch/?action=get_file&id=194&resource_link_id=403)), in denen das geeignete Luchshabitat weitgehend und dauerhaft von der Art besiedelt ist, wird ein deterministisches Fotofallen-Monitoring in sogenannten Referenzgebieten durchgeführt. Diese sind so gewählt, dass sie für die jeweiligen Teil-Kompartimente repräsentativ sind (Abb. 1; vgl. [hier](https://www.kora.ch/de/projekte/monitoring-grossraubtiere/was-ist-monitoring) für weitere Informationen). In jedem Referenzgebiet werden die Fotofallen systematisch verteilt und für einen bestimmten Zeitraum, in unserem Fall in der Regel für 60 Nächte, stehen gelassen. Diese Untersuchungen werden alle drei bis vier Jahre von KORA mit Hilfe der Kantone durchgeführt, um die Dichte und Abundanz der Luchspopulation mittels der so genannten fotografischen Fang-Wiederfang-Methode zu schätzen. Dies ist heute eine Standard-Monitoring-Methode für kryptische Arten, deren Individuen durch natürliche Markierungen erkennbar sind (vgl. [hier](https://www.kora.ch/de/projekte/monitoring-grossraubtiere/deterministisches-fotofallen-monitoring) für weitere Informationen)

Dieser Bericht präsentiert die Ergebnisse des deterministischen Durchgangs im Referenzgebiet ",Text.table[1,Refarea+1]," innerhalb des Teil-Kompartiments ",names(Text.table)[Refarea+1]," für den Winter ", Wintername ,". Die Ergebnisse vorangegangener Durchgänge sind [online](https://www.kora.ch/de/arten/luchs/fotofallen-monitoring) auf der KORA-Website verfügbar.
\n",sep="")

cat("![**Abb. 1.** Verteilung der Referenzgebiete (blaue Polygone) innerhalb der 16 Teil-Kompartimente (schwarze Polygone). Das geeignete Luchshabitat (hellgrün) ist definiert durch ein Luchs-Habitat-Modell (Zimmermann 2004). Das in diesem Durchgang untersuchte Referenzgebiet (",Text.table[1,Refarea+1],") sowie dessen Teil-Kompartiment (",names(Text.table)[Refarea+1],"), sind farblich hervorgehoben.](Plot_Introduction_Compartment.jpeg) 

\n",sep="")

\newpage

cat("### Matériel et Méthode \n")

cat("Au total,",Sites," sites choisis avec l’aide des surveillants de la faune ont été échantillonnés durant 60 nuits du ",format(as.Date(Start), "%d.%m.%Y")," au ", format(as.Date(Stop), "%d.%m.%Y")," à l'aide de ",Sites*2," pièges-photos (2 par sites) placés principalement le long de routes forestières et de chemins pédestres. Si un site a été déplacé lors de la session, l'ancien site ainsi que le nouveau sont pris en compte dans la représentation cartographique.",Text.table[6,names(Text.table)[Refarea+1]],"

L’aire de référence présente une superficie de ",  prettyNum(as.numeric(round(sf::st_area(Refareashp[Refarea,]) /1000000,0)), big.mark="'")," km^2^, dont ", prettyNum(favorable, big.mark="'")," km^2^ d'habitat favorable au lynx (Zimmermann 2004) (Fig.2). L’unité de la taille de la population correspond au nombre de lynx âgés de plus d’un an (lynx indépendants), c’est-à-dire les lynx adultes résidents et les lynx subadultes qui se dispersent. Les lynx juvéniles sont aussi identifiés, mais ne sont pas pris en compte dans l’estimation d’abondance et de densité à cause de leur faible taux de détection et de leur fort taux de disparition (mortalité et dispersion). ")

if(JuvasMother==TRUE){
  cat("**Ils sont toutefois comptabilisés comme une capture de leur mère respective dans l'historique de capture.** 

      \n")
}

if(JuvasMother==FALSE){
  cat("**Comme les lynx juvéniles se séparent fréquemment de leur mère entre février et avril, ils ne sont pas comptabilisés comme une capture de leur mère dans cette session.** 

      \n")
}

cat("![**Fig. 2.** Répartitions des sites (disques blancs) dans l'aire de référence ",Text.table[7,names(Text.table)[Refarea+1]]," (polygone bleu) au sein du sous-compartiment de gestion ",names(Text.table)[Refarea+1]," (polygone noir). En vert clair l’habitat favorable du lynx défini par un modèle d’habitat (Zimmermann 2004).](Plot_Introduction.jpeg)")
cat("### Material und Methode \n")

cat("Insgesamt wurden",Sites," Standorte mit Hilfe der Wildhüter ausgewählt und mit ",Sites*2," Fotofallen bestückt (2 pro Standort). Die Fotofallen liefen während 60 Nächten vom ",format(as.Date(Start), "%d.%m.%Y")," bis ",format(as.Date(Stop), "%d.%m.%Y")," und waren hauptsächlich entlang von Forstrassen und Wanderwegen aufgestellt. Wenn ein Standort während dem Durchgang versetzt wurde, werden sowohl der alte als auch der neue Standort in den kartographischen Darstellungen dieses Berichts berücksichtigt.",Text.table[5,names(Text.table)[Refarea+1]], "

Das Referenzgebiet hat eine Fläche von ",prettyNum(as.numeric(round(sf::st_area(Refareashp[Refarea,]) /1000000,0)), big.mark="'")," km^2^, davon sind ",prettyNum(favorable, big.mark="'")," km^2^ geeignetes Luchshabitat (Zimmermann 2004) (Abb.2). Die Einheit der Populationsgrösse ist die Anzahl der Luchse, die älter als ein Jahr sind (selbständige Luchse). Das heisst, residente adulte Luchse und noch nicht sesshafte subadulte Luchse werden in der Analyse berücksichtigt. Jungluchse werden zwar ebenfalls individuell bestimmt, aber aufgrund ihrer geringen Erfassbarkeit und ihrer hohen Verschwinderate (Mortalität und Dispersal) nicht  individuell in die Schätzungen von Abundanz und Dichte mit einbezogen. ")


if(JuvasMother==TRUE){
  cat("**Sie werden jedoch in der Fanggeschichte als Fang der jeweiligen Mutter gezählt.** 

      \n")
}

if(JuvasMother==FALSE){
  cat("**Da sich die Jungluchse zwischen Februar und April häufig von ihren Müttern trennen, werden die Jungtiere bei Durchgängen, die während dieser Zeit (zweite Winterhälfte) stattfinden, nicht als Fang ihrer Mutter gezählt.** 

      \n")
}

cat("![**Abb. 2.** Verteilung der Standorte (weisse Kreise) im Referenzgebiet ",Text.table[1,names(Text.table)[Refarea+1]]," (blaues Polygon) innerhalb des Teil-Kompartiments ",names(Text.table)[Refarea+1]," (schwarzes Polygon). In hellgrün das geeignete Luchshabitat, definiert durch ein Luchs-Habitat-Modell (Zimmermann 2004).](Plot_Introduction.jpeg){width=100%}")

\newpage

cat("### Résultats et Discussion \n")


PercentTN<-round((effective.TN*100)/potentially.TN,1)

if(PercentTN>=95){
  ValuePercentTN<-" **Cette valeur se trouve à la limite supérieure de celles observées dans d’autres études avec les pièges-photos où nous avions des valeurs comprises entre 84,2% (Nord du Jura, hiver 2006/07) et 99,9% (Simme-Saane, hiver 2020/21).**"
}

if(PercentTN<95 & PercentTN>90){
  ValuePercentTN<-" **Cette valeur est raisonnable et se trouve dans la fourchette des valeurs observées dans d’autres études avec les pièges-photographiques où nous avions des valeurs comprises entre 84,2% (Nord du Jura, hiver 2006/07) et 99,9% (Simme-Saane, hiver 2020/21).**"
}

if(PercentTN<90 & PercentTN>=84.2){
  ValuePercentTN<-" **Cette valeur se trouve à la limite inférieure des valeurs observées dans d’autres études avec les pièges-photos où nous avions des valeurs comprises entre 84,2% (Nord du Jura, hiver 2006/07) et 99,9% (Simme-Saane, hiver 2020/21).**"
}

if(PercentTN<84.2){
  ValuePercentTN<-" **Cette valeur est inférieure aux valeurs observées dans d’autres études avec les pièges-photos où nous avions des valeurs comprises entre 84,2% (Nord du Jura, hiver 2006/07) et 99,9% (Simme-Saane, hiver 2020/21).**"
}

cat("La durée d’échantillonnage potentielle était de ",prettyNum(potentially.TN, big.mark="'")," nuits de captures. **Le vol, le vandalisme, des problèmes techniques, des erreurs de manipulation et la neige** ont ramené l’effort d’échantillonnage à ", prettyNum(effective.TN, big.mark="'")," nuits-pièges effectives, soit ", gsub("[.]",",",sprintf("%.1f",PercentTN)) ,"% du potentiel.",ValuePercentTN,"**Les sites positifs sont répartis de façon uniforme sur l'ensemble de la zone d'étude (Fig. 3) (Describe repartition).** 

\n")


cat("![**Fig. 3.** Aire de référence ",Text.table[7,names(Text.table)[Refarea+1]]," (polygone bleu) avec la distribution spatiale des lynx (Polygone Convexe Minimum + zone tampon) photographiés durant l'ensemble de la session déterministe. Bleu\ua0: mâle (zone tampon de 1,4 km), rose\ua0: femelle (1,2 km), noir\ua0: sexe inconnu (1 km). Disques blancs avec un point noir\ua0: sites où au moins une photo d’un individu indépendant a été prise\ua0; **disques blancs avec un point rouge\ua0: sites où au moins une photo de lynx a été prise, mais l'individu n'a pas pu être identifié avec certitude\ua0; disques blancs avec un point vert\ua0: sites où au moins une photo de lynx a été prise, mais l'individu était un juvénile**\ua0; disques blancs\ua0: pas de photo de lynx.](Map_Report.jpeg)",sep="")
cat("### Resultate und Diskussion \n")

PercentTN<-round((effective.TN*100)/potentially.TN,1)

if(PercentTN>=95){
  ValuePercentTN<-" **Dieser Wert liegt im oberen Bereich anderer Fotofallen-Untersuchungen, wo wir Werte zwischen 84,2% (Jura Nord, Winter 2006/07) und 99,9% (Zentralschweiz West, Winter 2020/21) hatten.**"
}

if(PercentTN<95 & PercentTN>90){
  ValuePercentTN<-" **Dieser Wert liegt im Bereich anderer Fotofallen-Untersuchungen, wo wir Werte zwischen 84,2% (Jura Nord, Winter 2006/07) und 99,9% (Zentralschweiz West, Winter 2020/21) hatten.**"
}

if(PercentTN<90 & PercentTN>=84.2){
  ValuePercentTN<-" **Dieser Wert liegt im unteren Bereich anderer Fotofallen-Untersuchungen, wo wir Werte zwischen 84,2% (Jura Nord, Winter 2006/07) und 99,9% (Zentralschweiz West, Winter 2020/21) hatten.**"
}

if(PercentTN<84.2){
  ValuePercentTN<-" **Dieser Wert liegt unter den Bereich anderer Fotofallen-Untersuchungen, wo wir Werte zwischen 84,2% (Jura Nord, Winter 2006/07) und 99,9% (Zentralschweiz West, Winter 2020/21) hatten.**"
}

cat("Die potenzielle Anzahl der Fallennächte lag bei ",prettyNum(potentially.TN, big.mark="'"),". **Diebstahl, Vandalismus, technische Probleme, Bedienungsfehler und Schneefall** reduzierten den Aufwand der Datenerhebung auf ",prettyNum(effective.TN, big.mark="'")," tatsächliche Fallennächte, was ",gsub("[.]",",",sprintf("%.1f",PercentTN)),"% des Potenzials entspricht.",ValuePercentTN," **Die positiven Standorte sind gleichmässig über das gesamte Untersuchungsgebiet verteilt (Abb. 3) (Beschreibung der Verteilung).** 

\n",sep="")

cat("![**Abb. 3.** Referenzgebiet ",Text.table[1,names(Text.table)[Refarea+1]]," (blaues Polygon) mit räumlicher Verteilung von während dem Durchgang fotografierten Luchsen (kleinste Konvexpolygone + Pufferzone). Blau: Männchen (1,4 km Pufferzone), rosa: Weibchen (1,2 km), schwarz: unbekanntes Geschlecht (1 km). Weisse Kreise mit einem schwarzen Punkt: Standorte, an denen mindestens ein Foto eines selbständigen Luchses gemacht wurde; **weisse Kreise mit einem roten Punkt: Standorte, an denen lediglich mindestens ein Foto eines nicht individuell identifizierbaren Luchses gemacht wurde; weisse Kreise mit einem grünen Punkt: Standorte, an denen lediglich mindestens ein Foto eines juvenilen Luchses gemacht wurde**; weisse Kreise ohne Punkt: Standorte ohne Luchsfotos.](Map_Report.jpeg)")

\newpage

cat("#### Nombre minimum de lynx \n")

cat("**En tout, ", N.ind.lynx," lynx indépendants et ", N.juv," juvéniles de ", length(unique(Lynx.Table[Lynx.Table$YOB==juv.year & Lynx.Table$mother!="",mother]))," portées ont été détectés (voir Tab. 1)**. Ces lynx ont été photographiés sur ", as.numeric(Event[ID=="Total",4])," des ",Sites," sites (", prettyNum(as.numeric(round((Event[ID=="Total",4]*100)/Sites,1)),decimal.mark = ",") ,"%) échantillonnés dans le cadre du monitoring. ",sep="")




if(length(Not.seen$ID)>0){
cat("**Lors de cette session (60 nuits),",
length(Not.seen$ID),
" lynx ",
as.character(paste('(',paste(Not.seen$ID,collapse = ', '),')',sep='')),
"ont été observés dans le cadre du monitoring opportuniste sans avoir été photographié dans le cadre du monitoring déterministe.**

\n")
}else{cat("

\n")}

if(length(Not.seen.opp)>0){
cat("**De plus,",
length(Not.seen.opp),
" lynx ",
as.character(paste('(',paste(sort(Not.seen.opp),collapse = ', '),')',sep='')),
"ont été observés dans le cadre de ce monitoring sans avoir été photographié dans le cadre du monitoring opportuniste. Les indices de présence de chaque individu peuvent être consultés dans le [KORA Monitoring Center](https://www.koracenter.ch).**

\n")
}

cat("Un total de ",Event[ID=="Total",Event]," événements ont été répertoriés au cours de la période d’échantillonnage de 60 nuits dans l’ensemble de l’aire de référence (Tab. 1). Ceux-ci ont été ramenés à",
ifelse(Period.selected==3,Occ.table.plot.3[20,SUM.occasion.cumulative], Occ.table.plot.5[12,SUM.occasion.cumulative])
," par le groupement des détections par occasion (est égale à", ifelse(Period.selected==3,3,5) ,"nuits consécutives lors de cette étude). Si un lynx est photographié à plusieures reprises pendant une même occasion, il compte comme une seule capture. **Vu que le nombre de captures cumulées augmente constamment et que les pièges-photos ont fonctionnés normalement, nous pouvons raisonnablement conclure que les lynx n'ont pas commencé à éviter les sites au cours de la période d'échantillonnage (Fig. 4). Cette augmentation continue est pour le monitoring un indice de bonne qualité. De plus la majeure partie des lynx présents dans l’aire de référence ont été détectés vu que le nombre de lynx capturés (courbe bleu, Fig. 4) se stabilise déjà à partir de la XXXXième occasion.** \n")

if(length(Seen.external)==0){
cat("**Lors de cette session, il n’y a pas eu de lynx qui ont été détectés en dehors de la période de 60 jours. Tous les lynx détectés ont donc été comptabilisés dans le calendrier de capture en vue des analyses de capture-recapture.** 

\n")
}else{
cat("**Lors de cette session, le(s) lynx ", paste(Seen.external, collapse = ", "),"a/ont uniquement été détecté(s) en dehors de la période de 60 jours. Ce(s) lynx apparaisse(nt) sur la carte (Fig. 3.) mais ne sont pas pris en compte dans le reste des analyses. Cependant, ces lynx sont toujours pris en compte dans l'abondance estimée par capture-recapture au moyen de l'intervalle de confiance (voir chapitre «Estimation de l'abondance et de la densité»).** 

\n") 

}

cat("![**Fig. 4.** Evolution du nombre de captures de lynx cumulées et du nombre de lynx capturés cumulés dans l'aire de référence ",Text.table[7,Refarea+1],".

\n](Plot_cumulative_report.jpeg)",sep="")
cat("#### Minimale  Anzahl  Luchse \n")

cat("**Insgesamt wurden ", N.ind.lynx," selbständige Luchse und ", N.juv," Jungtier(e) von ", length(unique(Lynx.Table[Lynx.Table$YOB==juv.year & Lynx.Table$mother!="",mother]))," Wurf/Würfen nachgewiesen (Tab. 1)**. Diese Luchse wurden an ", as.numeric(Event[ID=="Total",4])," der ",Sites," aufgestellten Standorte (", prettyNum(as.numeric(round((Event[ID=="Total",4]*100)/Sites,1)),decimal.mark = ","),"%) fotografiert. ",sep="")

if(length(Not.seen$ID)>0){
cat("**Während diesem Durchgang (60 Tage) wurde(n)",
sum(!grepl("XXXX",Not.seen$ID),na.rm = TRUE),
" Luchs(e) ",
as.character(paste('(',paste(Not.seen$ID,collapse = ', '),')',sep='')),
"im Rahmen des opportunistischen Monitorings im untersuchten Referenzgebiet fotografiert, welche/r im Rahmen des deterministischen Monitorings nicht fotografiert wurde(n).**

\n")
}else{cat("

\n")}

if(length(Not.seen.opp)>0){
cat("**Zusätzlich wurde(n)",
length(Not.seen.opp),
"Luchs(e)",
as.character(paste('(',paste(sort(Not.seen.opp),collapse = ', '),')',sep='')),
"während dem Monitoring erfasst, ohne dass sie im selben Zeitraum im Rahmen des opportunistischen Monitorings fotografiert wurden. Die Nachweise von jedem einzelnen Luchsindividuum können im [KORA Monitoring Center](https://www.koracenter.ch) aufgerufen werden.**

\n")
}

cat("Im gesamten Referenzgebiet wurden während der 60 Nächte des deterministischen Fotofallen-Monitorings bei ",Event[ID=="Total",Event]," Ereignissen Luchse fotografiert (Tab. 1). Diese Ereignisse werden für die Analyse in Fanggelegenheiten eingeteilt (in dieser Studie ", ifelse(Period.selected==3,3,5) , "aufeinanderfolgende Nächte) woraus ",ifelse(Period.selected==3,Occ.table.plot.3[20,SUM.occasion.cumulative], Occ.table.plot.5[12,SUM.occasion.cumulative])," Erfassungen hervorgehen. Wird derselbe Luchs während einer Fanggelegenheit mehrmals fotografiert, zählt dies als eine Erfassung. **Die Zahl der kumulativen Erfassungen nimmt ständig zu (Abb. 4). Dieser kontinuierliche Anstieg ist ein Indiz für einen erfolgreichen Durchgang. Da es keine grossflächigen Einbrüche in der Funktionalität der Fotofallen gab, können wir daraus schliessen, dass Luchse keine Scheu vor den Fotofallen entwickelten. Des Weiteren stabilisiert sich die Anzahl der verschiedenen Luchse (blaue Kurve, Abb. 4) bereits ab der XXXXX Fanggelegenheit, weshalb davon ausgegangen werden kann, dass der Grossteil der Luchse im Referenzgebiet nachgewiesen werden konnte.** \n")


if(length(Seen.external)==0){
cat("**Es gab in diesem Durchgang keine Luchse, welche lediglich ausserhalb des 60-Tage-Zeitraums fotografiert wurden, alle erfassten Individuen werden somit in der Fang-Wiederfang-Analyse berücksichtigt.**

\n")
}else{
cat("**Während diesem Durchgang wurde(n) der/die Luchs(e) ", paste(Seen.external, collapse = ", "), "nur ausserhalb des 60-Tage-Zeitraums entdeckt. Diese Luchse erscheinen auf der Karte (Abb. 3), werden aber in den übrigen Analysen nicht berücksichtigt. Mittels Konfidenzintervall werden diese Luchse aber immer noch in der mittels Fang-Wiederfang geschätzten Abundanz berücksichtigt (siehe Kapitel «Schätzung der Abundanz und Dichte»).**

\n") 

}

cat("![**Abb. 4.** Entwicklung der kumulierten Anzahl Erfassungen und der kumulierten Anzahl verschiedener fotografierter Luchse im Referenzgebiet ",Text.table[1,Refarea+1],". 

\n](Plot_cumulative_report.jpeg)",sep="")

\newpage

Event.report

\newpage

cat("#### Estimation de l'abondance et de la densité \n")


#Description du modèle:
model.descritption<-data.table(Model=c("0","b","h","t","bh"),
                               Description=c("**qui considère que la probabilité de capture est constante**", #M0
                      "**qui considère que la probabilité de capture varie dû à une réponse biologique**", #Mb
                      "**qui considère que la probabilité de capture varie d'un individu à l'autre**", #Mh
                      "**qui considère que la probabilité de capture varie en fonction du temps**",#Mt
                      "**qui considère que les individus ont des probabilités de capture différentes, qui changent après leur première capture**")) #Mbh

model.descritption<-model.descritption[model.descritption$Model==substring(Mark.Table.inline[which.best.model.selected,Model], 2),Description]


cat("Le modèle ",paste("M~",substr(Mark.Table.inline[which.best.model.selected,Model], 2, 3),"~",sep=""),", ",model.descritption,", explique le mieux les données. L’estimation de l’abondance (intervalle de confiance de 95%) par capture-recapture était de ",
paste(Mark.Table.inline[which.best.model.selected,N])," (",sub("\\-.*", "", Mark.Table.inline[which.best.model.selected,IC])," - ",sub("..-*", "", Mark.Table.inline[which.best.model.selected,IC]),") lynx indépendants. Ainsi, ", prettyNum(as.numeric(sprintf("%.1f",round((length(Lynx.Table.Adult$ID)*100)/Mark.Table.inline[which.best.model.selected,N],1))),decimal.mark = ","),"% des lynx estimés ont effectivement été photographiés. 
\n",sep="")
cat("#### Schätzung der Abundanz und Dichte \n")

#Description du modèle:
model.descritption<-data.table(Model=c("0","b","h","t","bh"),
                               Description=c("**bei dem davon ausgegangen wird, dass die Wahrscheinlichkeit des Fangens konstant ist**", #M0
                      "**bei dem davon ausgegangen wird, dass die Fangwahrscheinlichkeit aufgrund einer biologischen Reaktion variiert**", #Mb
                      "**welches unterschiedliche individuelle Fangwahrscheinlichkeiten zulässt**", #Mh
                      "**bei dem davon ausgegangen wird, dass die Wahrscheinlichkeit des Fangens mit der Zeit variiert**",#Mt
                      "**bei dem davon ausgegangen wird, dass Individuen unterschiedliche Fangwahrscheinlichkeiten haben, die sich nach ihrem ersten Fang ändern**")) #Mbh

model.descritption<-model.descritption[model.descritption$Model==substring(Mark.Table.inline[which.best.model.selected,Model], 2),Description]


cat("Das Modell ",paste("M~",substr(Mark.Table.inline[which.best.model.selected,Model], 2, 3),"~",sep=""),", ",model.descritption,", erklärt die Daten am besten. Die resultierende geschätzte Abundanz (95% Konfidenzintervall) anhand dieses Fang-Wiederfang-Modells ist ",
paste(Mark.Table.inline[which.best.model.selected,N])," (",sub("\\-.*", "", Mark.Table.inline[which.best.model.selected,IC])," - ",sub("..-*", "", Mark.Table.inline[which.best.model.selected,IC]),")  selbständige  Luchse. ", prettyNum(as.numeric(sprintf("%.1f",round((length(Lynx.Table.Adult$ID)*100)/Mark.Table.inline[which.best.model.selected,N],1))),decimal.mark = ","),"% der geschätzten Luchse wurden demnach tatsächlich fotografiert. 
\n", sep="")
#Plot history plot or (if first session in area) write just the text.
if(length(History$Year)>1){
cat("La densité (intervalle de confiance de 95%) dans l’aire de référence était de ",
    prettyNum(Density.surface,decimal.mark = ",")," (",prettyNum(Density.surface.min95,decimal.mark = ","),"-",prettyNum(Density.surface.max95,decimal.mark = ","),") lynx indépendants pour 100 km^2^ ou ",
      prettyNum(Density.surface.favorable,decimal.mark = ",")," (",prettyNum(Density.surface.favorable.min95,decimal.mark = ","),"-",prettyNum(Density.surface.favorable.max95,decimal.mark = ","),") lynx indépendants pour 100 km^2^ d’habitat favorable (Fig.5). 

\n", sep="")


cat("![**Fig. 5.** Evolution des densités de lynx indépendants pour 100 km^2^ d’habitat favorable (avec intervalle de confiance de 95%) dans l’aire de référence ",Text.table[7,Refarea+1],". En orange la session actuelle.

\n](History_denisty_plot.jpeg)",sep="")
}

if(length(History$Year)==1){
cat("La densité (intervalle de confiance de 95%) dans l’aire de référence était de ",
    prettyNum(Density.surface,decimal.mark = ",")," (",prettyNum(Density.surface.min95,decimal.mark = ","),"-",prettyNum(Density.surface.max95,decimal.mark = ","),") lynx indépendants pour 100 km^2^ ou ",
      prettyNum(Density.surface.favorable,decimal.mark = ",")," (",prettyNum(Density.surface.favorable.min95,decimal.mark = ","),"-",prettyNum(Density.surface.favorable.max95,decimal.mark = ","),") lynx indépendants pour 100 km^2^ d’habitat favorable. 

\n", sep="")

}
#Plot history plot or (if first session in area) write just the text.
if(length(History$Year)>1){

cat("Die Dichte (95% Konfidenzintervall) im Referenzgebiet betrug ",
    prettyNum(Density.surface,decimal.mark = ",")," (",prettyNum(Density.surface.min95,decimal.mark = ","),"-",prettyNum(Density.surface.max95,decimal.mark = ","),") selbständige Luchse pro 100 km^2^ oder ",
      prettyNum(Density.surface.favorable,decimal.mark = ",")," (",prettyNum(Density.surface.favorable.min95,decimal.mark = ","),"-",prettyNum(Density.surface.favorable.max95,decimal.mark = ","),") selbständige Luchse pro 100 km^2^ geeignetem Habitat (Abb. 5). 

\n",sep="")

cat("![**Abb. 5.** Entwicklung der Luchsdichte pro 100 km^2^ geeignetem Habitat (mit 95 % Konfidenzintervall) im Referenzgebiet ",Text.table[1,Refarea+1],". In orange der aktuelle Durchgang.

\n](History_denisty_plot.jpeg)",sep="")

}

if(length(History$Year)==1){
cat("Die Dichte (95% Konfidenzintervall) im Referenzgebiet betrug ",
    prettyNum(Density.surface,decimal.mark = ",")," (",prettyNum(Density.surface.min95,decimal.mark = ","),"-",prettyNum(Density.surface.max95,decimal.mark = ","),") selbständige Luchse pro 100 km^2^ oder ",
      prettyNum(Density.surface.favorable,decimal.mark = ",")," (",prettyNum(Density.surface.favorable.min95,decimal.mark = ","),"-",prettyNum(Density.surface.favorable.max95,decimal.mark = ","),") selbständige Luchse pro 100 km^2^ geeignetem Habitat. 

\n",sep="")

}

\newpage

cat("#### Comparaison avec les densités estimées dans les autres aires de référence 

\n")

cat("**La densité estimée pour 100 km^2^ d’habitat favorable dans le ",Text.table[7,Refarea+1]," (",
      paste(prettyNum(History[length(History$Year),Density],decimal.mark=",")," (",
      prettyNum(History[length(History$Year),Min95],decimal.mark=","),
      "-",
      prettyNum(History[length(History$Year),Max95],decimal.mark=","),")",sep=""),") se situe parmi les (XXXX-check) plus haute estimées dans l’ensemble de la Suisse (Tab. 2).**",sep="")
cat("#### Vergleich mit geschätzten Dichten in den anderen Referenzgebieten 

\n")

cat("**Die geschätzte Dichte pro 100 km^2^ geeignetem Habitat im Referenzgebiet ",Text.table[1,Refarea+1]," (",
      paste(prettyNum(History[length(History$Year),Density],decimal.mark=",")," (",
      prettyNum(History[length(History$Year),Min95],decimal.mark=","),
      "-",
      prettyNum(History[length(History$Year),Max95],decimal.mark=","),")",sep=""),") liegt im (XXXXX-check) mittleren Bereich der zuletzt geschätzten Werte in den Referenzgebieten (Tab. 2).**",sep="")
History.comparaison

\newpage

cat("### Remerciements \n")

cat("Nous remercions vivement tous ceux qui d’une manière ou d’une autre nous ont aidés et
soutenus. En particulier\ua0: 

\n")

cat("- tous les responsables des institutions cantonales et fédérales notamment", 
    paste(Text.table[2,names(Text.table)[Refarea+1]]),
    paste(Text.table[4,names(Text.table)[Refarea+1]]),"pour leur soutien professionnel\ua0;

    \n")

cat("- tous les surveillants de  la faune et volontaires qui nous ont  aidés lors du choix  des sites ainsi  que  lors  de  la mise  en place  des pièges-photos, les contrôles et le  démontage,  en particulier\ua0:", paste(gsub(",","",sort(unique(DataFULL[,site_controller])[!unique(DataFULL[,site_controller]) %in% "KORA"])),collapse=", "),"\ua0;

\n")




cat("- tous les civilistes, stagiaires et collaborateurs au KORA qui ont participé au projet\ua0: **XXXXX**.

\n")

cat("### Références 

\n")

cat("Laass J.1999. *Evaluation von Photofallen für ein quantitatives Monitoring einer Luchspopulation in den Schweizer Alpen.* Universität Wien, Wien. 

\n")

cat("R Core Team. 2021. *A Language and Environment for Statistical Computing*. R Foundation for Statistical
  Computing, Vienna, Austria. URL\ua0: https://www.R-project.org/. 

\n")

cat("Zimmermann F. 2004.  *Conservation  of  the  Eurasian  lynx  (Lynx  lynx) in  a  fragmented  landscape  –  habitat  models, dispersal,  and  potential distribution*.  PhD  Thesis,  Department  of  Ecology  and  Evolution,  University  of  Lausanne, Switzerland. 

\n")

cat("**Zimmermann F. & Breitenmoser U. 2007. *Potential distribution and population size of the Eurasian lynx (Lynx lynx) in the Jura Mountains and possible corridors to adjacent ranges.* Wildlife Biology 13\ua0: 406–416.** 

\n")


cat("**Citation proposée\ua0:** 
\n")

#Change Authors name system Last name F. (Le Grand L.)
Authors<-unlist(strsplit(Authors, ", "))

Authors.abr<-Authors
for(i in 1:length(Authors)){

if(length(unlist(strsplit(Authors[i]," ")))>2){
Authors.abr[i]<-paste(paste(unlist(strsplit(Authors[i]," "))[1:2],collapse=" ")," ", abbreviate(unlist(strsplit(Authors[i]," "))[3],1),".",sep="")
}else{
Authors.abr[i]<-paste(unlist(strsplit(Authors[i]," "))[1]," ", abbreviate(unlist(strsplit(Authors[i]," "))[2],1),".",sep="")
}
}
Authors<-paste(Authors.abr,collapse=", ")

cat(paste(Authors,", ",data.table::year(Sys.time()),". ",Title," KORA Bericht xx, 12pp. 
\n",sep=""))


cat("**Géodonnées\ua0:** 
\n")

cat("Toutes les analyses et le traitement des données ont été effectués à l'aide du langage et de l'environnement de programmation statistique R 4.1.0 (R Core Team, 2021). Les données utilisées pour créer le fond de carte ont été extraites de l'Open Street Map (https://www.openstreetmap.org/). La Figure 1 contient la couche altitudinale et la couche des lacs de GEOSTAT, © Bundesamt für Statistik; Euromaps, © Bartholomew.
\n")
cat("### Danksagung \n")

cat("Wir danken allen ganz herzlich, die uns bei der Durchführung des deterministischen Fotofallen-Durchgangs in irgendeiner Form unterstützt haben. Besonders danken wir: 

\n")

cat("- allen  Verantwortlichen der beteiligten kantonalen  und eidgenössischen  Institutionen, namentlich", 
    paste(Text.table[2,names(Text.table)[Refarea+1]]),
    paste(Text.table[3,names(Text.table)[Refarea+1]]),"für ihre  professionelle  Unterstützung;

    \n")

cat("- allen Wildhütern und Freiwilligen, die  uns bei der Wahl der  Standorte,  sowie  bei den Kontrollen und  dem Abbau der Fotofallen geholfen haben, insbesondere: ", paste(gsub(",","",sort(unique(DataFULL[,site_controller])[!unique(DataFULL[,site_controller]) %in% "KORA"])),collapse=", "),";

\n")

cat("- allen Zivildienstleistenden, Praktikant/innen und Mitarbeiter/innen von KORA, die an dem Projekt teilgenommen haben: **XXXXX**.

\n")


cat("### Referenzen 

\n")

cat("Laass, J. 1999. *Evaluation von Photofallen für ein quantitatives Monitoring einer Luchspopulation in den Schweizer Alpen.* Universität Wien, Wien. 

\n")

cat("R Core Team.  2021. *A Language and Environment for Statistical Computing*. R Foundation for Statistical
  Computing, Vienna, Austria. URL: https://www.R-project.org/. 

\n")

cat("Zimmermann,  F.  2004.  *Conservation  of  the  Eurasian  lynx  (Lynx  lynx) in  a  fragmented  landscape  –  habitat  models, dispersal,  and  potential distribution*.  PhD  Thesis,  Department  of  Ecology  and  Evolution,  University  of  Lausanne, Switzerland. 

\n")

cat("**Zimmermann, F. & Breitenmoser, U. 2007. *Potential distribution and population size of the Eurasian lynx (Lynx lynx) in the Jura Mountains and possible corridors to adjacent ranges.* Wildlife Biology 13: 406–416.** 

\n")

cat("**Vorgeschlagene Zitierung:** 
\n")

#Change Authors name system Last name F. (Le Grand L.)
Authors<-unlist(strsplit(Authors, ", "))

Authors.abr<-Authors
for(i in 1:length(Authors)){

if(length(unlist(strsplit(Authors[i]," ")))>2){
Authors.abr[i]<-paste(paste(unlist(strsplit(Authors[i]," "))[1:2],collapse=" ")," ", abbreviate(unlist(strsplit(Authors[i]," "))[3],1),".",sep="")
}else{
Authors.abr[i]<-paste(unlist(strsplit(Authors[i]," "))[1]," ", abbreviate(unlist(strsplit(Authors[i]," "))[2],1),".",sep="")
}
}
Authors<-paste(Authors.abr,collapse=", ")

cat(paste(Authors,", ",data.table::year(Sys.time()),". ",Title," KORA Bericht xx, 12pp.

\n",sep=""))


cat("**Digitale geografische Daten:** 

\n")

cat("Alle Datenanalysen und die Datenverarbeitung wurden mit der statistischen Programmiersprache und Umgebung R 4.1.0 (R Core Team, 2021) durchgeführt. Die Daten, die für die Erstellung des Kartenhintergrunds verwendet wurden, stammen aus der Open Street Map (https://www.openstreetmap.org/). Abbildung 1 enthält die Höhen- und Seenschicht von GEOSTAT, © Bundesamt für Statistik; Euromaps, © Bartholomew.
\n")

\newpage

cat("### Annexes \n")
cat("### Anhänge \n")
Species.other.table<-data.table(Latin=c("Canis lupus","Felis silvestris","Canis aureus"),
                                FR=c("loup","chat sauvage","chacal doré"),
                                FRs=c("loups","chats sauvages","chacals dorés"))

if(length(Species.other)>0){

for(i in 1:length(Species.other)){

Plot_KORA_Report(Species.other[i])

}

for(i in 1:length(Species.other)){

cat("![**Fig. A.",i,".** Détections de ",Species.other.table[Species.other.table$Latin==Species.other[i],FRs]," (*",Species.other[i],"*) dans l’aire de référence ",Text.table[7,Refarea+1]," (polygone bleu) dans le sous-compartiment ",names(Text.table)[Refarea+1]," (polygone noir). Disques blancs avec un point noir\ua0: sites où au moins une photo de ",Species.other.table[Species.other.table$Latin==Species.other[i],FR]," a été prise\ua0; disques blancs\ua0: pas de photo de ",Species.other.table[Species.other.table$Latin==Species.other[i],FR],".

\n](",paste(Species.other[i],".jpeg",sep=""),")

\n",sep="") 
}
}
Species.other.table<-data.table(Latin=c("Canis lupus","Felis silvestris","Canis aureus",
                                        "Vulpes vulpes", "Capreolus capreolus",
                                        "Rupicapra rupicapra",
                                        "Meles meles",
                                        "Lepus europaeus",
                                        "Cervus elaphus",
                                        "Capra ibex",
                                        "Martes foina",
                                        "Felis catus",
                                        "Lepus timidus",
                                        "Sciurus vulgaris",
                                        "Martes martes",
                                        "Erinaceus europaeus"),
                                DE=c("Wolf","Wildkatze","Goldschakal",
                                     "Fuchs","Reh",
                                     "Gämse",
                                     "Dachs",
                                     "Feldhase",
                                     "Hirsch",
                                     "Steinbock",
                                     "Steinmarder",
                                     "Hauskatze",
                                     "Schneehase",
                                     "Eichhörnchen",
                                     "Baummarder",
                                     "Igel"))

if(length(Species.other)>0){for(i in 1:length(Species.other)){

Plot_KORA_Report(Species.other[i])


 cat("![**Abb. A.",i,".** Nachweise ",Species.other.table[Species.other.table$Latin==Species.other[i],DE]," (*",Species.other[i],"*) im Referenzgebiet ",Text.table[1,Refarea+1]," (blaues Polygon) im Teil-Kompartiment ",names(Text.table)[Refarea+1]," (schwarzes Polygon).  Weisse Kreise mit einem schwarzen Punkt: Standorte mit ",Species.other.table[Species.other.table$Latin==Species.other[i],DE],"; weisse Kreise ohne Punkt: Standorte ohne ",Species.other.table[Species.other.table$Latin==Species.other[i],DE],".

\n](",paste(Species.other[i],".jpeg",sep=""),")

\n",sep="") 

}}

\newpage

KORA Appendix
This part can be removed from the report. It is produced for KORA people writing the report. When finishing the report, one should 1) check that all sites had at least one picture, 2) choose model and re-run the script, 3) check all sentences written in bold, 4) adapt the number of pages in the citation, 5) adapt the XX of the KORA Bericht xx. Always consider the output as a first draft and not as a final version. Every session could present its own exceptions chich will ask for an improvement/correction of the script.

\n

\n

Mark.Table

\n

\n

Opp.table

\newpage

Fig. KORA.1 Cumulative number of occasion (black) and cummulative number of lynx ID (blue). Left = 20 occasions, right = 12 occasions. Closure test should be > 0.05 if the population is closed. The blue line should flatten itself. This shows that we captured most of the lynx ID in the area. The black line should grow constantly. This shows that lynx do not avoid camera traps after their first capture.

\newpage

Fig. KORA.2 The opportunistic monitoring represents here the opportunistic observations that took place during the monitoring period and that were transmitted to KORA. Only observations made in the study area are selected. The blue line delimits the study area. In colour: areas of lynx activity (Minimum Convex Polygon + 2 km buffer zone); black dots: sites were lynx were seen; orange dots: sightings of unidentifiable lynx (U).

\newpage

Dispersal

\newpage

Table.Sites


CrocutaLupus/KORAtool documentation built on April 13, 2025, 10:03 p.m.