diff --git a/script_NI_LSK.Rmd b/script_NI_LSK.Rmd new file mode 100644 index 0000000000000000000000000000000000000000..ed426a6f1c873e6022735ecf69096f7813883836 --- /dev/null +++ b/script_NI_LSK.Rmd @@ -0,0 +1,487 @@ +--- +title: "Script for beregning av Naturindeks-verdier for indikatorer basert på data fra Landsskogtakseringen" +author: "Arvid Svensson & Ken Olaf Storaunet" +date: "Dec. 2023" +output: pdf_document +editor_options: + chunk_output_type: console +--- + + +# Innledning/Bakgrunn + +Ti av 89 indikatorer i Naturindeksen (NI) for skog er basert på Landsskogtakseringens (LSK) data: ‘23 Blåbær’, ‘41 Eldre lauvsuksesjon (MiS)’, ‘60 Gamle trær (MiS)’, ‘108 Liggende død ved (MiS) – arealandel’, ‘183 Stående død ved (MiS) – arealandel’, ‘207 Trær med hengelav (MiS)’, ‘418 Stående død ved – mengde’, ‘419 Liggende død ved – mengde’, ‘433 Rogn-Osp-Selje’, og ‘434 Gammel skog’. +Datauthenting, bearbeiding, beregninger og innlegging i naturindeks-databasen er tidligere gjort manuelt. + +Her har vi automatisert disse prosessene for de nevnte indikatorene. Beregninger av selve indikatorverdiene er ofte enkle gjennomsnittsverdier, men det er avgjørende å holde orden på ulike utvalg/strata av prøveflatene. + +Det vesentlige av kodene i LSKs database er beskrevet i LSKs takstinstruks (Viken 2021, https://nibio.brage.unit.no/nibio-xmlui/handle/11250/2826859). + + +# Innstillinger + +Inneholder diverse settings og options. + +'data_lokal' avgjør om data skal leses lokalt/testdata (default = T), eller '= F' om data skal leses fra LSK-databasen. +'periode' endres ihht. hvilken tidsperiode som skal leses fra LSK-databasen. "1418" betyr 5-årsperioden 2014-2018, og er default for de inkluderte testdata. + + +```{r, setup} + +############## +## Settings ## +############## + +# Empty Global Environment +rm(list = ls()) + +# Language setting +Sys.setlocale(locale = "no_NB.utf8") + +# Setting RMD output +knitr::opts_chunk$set(echo = T, fig.cap = FALSE,message=FALSE,warning = FALSE,comment = FALSE) + +# Variable for local-/test-data or update from LSK-db +data_lokal=T + +# Years LSK-data input +periode<-"1418" + +# Set Location R-libraries +# .libPaths("") + +# Set Working directory to source file location +setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) + +# Libraries +library(tidyverse) +library(boot) + +``` + +# Definerer Bootstrap Funksjonen + +I NI angis usikkerhet for indikatorverdiene som estimater/verdier for 1. og 3. kvartil. For LSK-indikatorene beregnes dette ved enkel ikkeparametrisk bootstrap, med boot- og boot.ci-funksjonene i R-pakken ‘boot’. Disse beregningene gjøres separat for hvert prøveflate-utvalg av indikator og 'Area_Name'. + +```{r} + + +################################### +## Define Function for Bootstrap ## +################################### + +my.function = function(data,index){ + + d = data[index,] # Create bootstrap sample of all columns of original data + return(weighted.mean(d$var_of_interest, d$AU_AREAL)) # Calculate weighted mean using 'counts' and 'weights' columns + +} + +# Bootstrap for all indicators/NI-regions +beregninger.boot <-function(p){ +sp <- split(p, p$region) +y <- lapply(sp, function(x){ + wtd.avg <- weighted.mean(x$var_of_interest, x$AU_AREAL) + perc <- boot.ci(boot(x, my.function, R = 1000), type = "perc", conf = 0.5)$perc + CI.LL <- perc[4] + CI.UL <- perc[5] + SESONG<-round(mean(as.numeric(x$SESONG))) + n=nrow(x) + data.frame(wtd.avg, CI.LL, CI.UL,SESONG,n ) +}) + +return(rownames_to_column(do.call(rbind, y))) +} + + +``` + +# Lese data lokalt eller fra LSK-databasen + +Dersom en skal hente data fra LSK-databasen trenger en R-pakken 'ROracle', samt brukertilgang. Dette gjøres ved NIBIO. + +Det defineres ulike strata basert på region/fylke (DISTRIKT), hvilket prøveflatenett (LS_NETT), produktiv og uproduktiv skog (AREALTYPE), og arealanvendelse (AREAL_ANV, der 1 = skog/utmark, 5 = naturreservat/nasjonalpark, 9 = friluftsområde). + +NB: Merk at dersom data for tidligere år (<2010?) skal hentes fra LSK-databasen, er koding for enkelte av disse kategoriene endret. + +Definisjonen/alderskravet for 'Biologisk gammel skog' (BGS) avhenger av dominerende treslag og bonitet. + +AREA_ID i NI er unik for hver enkelt indikator, men alle LSK-baserte indikatorer benytter den samme areal-inndelingen (hele eller delte fylker). Koblingen mellom FLATEID i LSK og AREA_ID i NI går gjennom to trinn og hentes fra vedlagte .csv-filer (1: LSKs 'FLATEID' og NIs 'Area_Name', 2: NIs 'Area_Name' og 'Area_ID' for den enkelte indikator). + +Det gjennomføres to kontroller spesifikt for MiS-dataene som potensielt stanser lesingen av data: +1) at alle prøveflater i skog har MiS-data +2) at det bare finnes 1 rad med data pr prøveflate og MiS-livsmiljø +Hvis dette stanser koden, må det feilsøkes i LSK-databasen. + + +```{r} + +# Use local-/test-data or update from LSK-db + +if(data_lokal==T){ + load("data_til_naturindeks.RData") +}else{ +# source("oppdater_data.R") +} + +``` + +# Beregning av dødved-mengde-indikatorene (ind_ID 418 og 419) + +```{r} +####################################################### +## Deadwood (volume) indicators (ind_ID 418 and 419) ## +####################################################### + +fla<-merge(dat_dod%>%select(-SESONG), + dta2[c("AU_AREAL","SkogIO","SESONG","HEL_DELT_FLATE","FLATEID","strata","region")], + by=c( "HEL_DELT_FLATE","FLATEID"), + all.y = T)%>% + filter(strata!="2") + +fla2<-fla%>%mutate_all(~replace(., is.na(.), 0))%>% + select(FLATEID,AU_AREAL,strata,AU_VG_LEGER,AU_VF_LEGER, + AU_VUKJ_LEGER,AU_VL_LEGER,AU_VMG_GADD,AU_VMF_GADD, + AU_VML_GADD,SkogIO,region,SESONG) + +fla_boot<-fla2%>% + mutate( + AU_vol_LEGER=AU_VG_LEGER+AU_VF_LEGER+AU_VUKJ_LEGER+AU_VL_LEGER, + AU_vol_GADD=AU_VMG_GADD+AU_VMF_GADD+AU_VML_GADD)%>% + select(FLATEID,AU_AREAL,strata,AU_vol_LEGER, + AU_vol_GADD,SkogIO,region,SESONG)%>% + mutate(vol_LEGER_ha=AU_vol_LEGER/AU_AREAL, + vol_GADD_ha=AU_vol_GADD/AU_AREAL, + SkogIO_2=case_when(SkogIO!=0~1, + TRUE~0))%>% + + group_by(strata,region, + FLATEID,SESONG,SkogIO_2)%>% + summarise(AU_vol_LEGER=sum(AU_vol_LEGER), + AU_vol_GADD=sum(AU_vol_GADD), + AU_AREAL=sum(AU_AREAL))%>% + mutate(vol_LEGER_ha=AU_vol_LEGER/AU_AREAL, + vol_GADD_ha=AU_vol_GADD/AU_AREAL)%>% + + ungroup(.)%>% + filter(SkogIO_2!=0)%>% as.data.frame() + + +# Change var-of-interest to ind_ID 419 'Fallen deadwood, m3/ha' +fla_boot<-fla_boot%>%mutate(var_of_interest=vol_LEGER_ha)%>% as.data.frame() +# Result ind_ID 419: 'Fallen deadwood, m3/ha' +Res_ind_419<-cbind("Liggende død ved – mengde",beregninger.boot(fla_boot)) + +# Change var-of-interest to ind_ID 418 'Standing deadwood, m3/ha' +fla_boot<-fla_boot%>%mutate(var_of_interest=vol_GADD_ha)%>% as.data.frame() +# Result ind_ID 418: 'Standing deadwood, m3/ha' +Res_ind_418<-cbind("Stående død ved – mengde",beregninger.boot(fla_boot)) + +``` + +# Beregning av indikatoren Gammel skog (ind_ID 434) + +```{r} +############################################# +## Old Forest (BGS) indicator (ind_ID 434) ## +############################################# + +fla<-dta2%>%mutate_all(~replace(., is.na(.),0))%>%filter(strata!="2") +fla2<-fla%>%group_by(strata,region,FLATEID,SESONG)%>% + filter(prodSkogIO!=0)%>% + summarise(prodSkogIO=sum(prodSkogIO),gammelSkogIO=sum(gammelSkogIO), + AU_AREAL=sum(AU_AREAL))%>% + ungroup(.)%>% + + mutate( + + andel_gammelskog=case_when(prodSkogIO!=0~100*(gammelSkogIO/prodSkogIO), + TRUE~0), + + prodSkogIO_2=case_when(prodSkogIO!=0~1, + TRUE~0) + + ) + +# Change var-of-interest to ind_ID 434 'Old Forest' +fla_boot<-fla2%>%filter(prodSkogIO_2==1)%>%mutate(var_of_interest=andel_gammelskog)%>% as.data.frame() +# Result ind_ID 434: 'Old Forest' +Res_ind_434<-cbind("Gammel skog",beregninger.boot(fla_boot)) + +``` + +# Beregning av indikatoren Rogn-Osp-Selje (ind_ID 433) + +```{r} +################################################# +## ROS-indicator (Rogn-Osp-Selje) (ind_ID 433) ## +################################################# + +fla<-merge(dta2,dat_tre,by=c("FLATEID","HEL_DELT_FLATE"),all.x = T)%>% + mutate_all(~replace(., is.na(.),0))%>% + mutate(VMT=`SUM(VMT)`)%>% + filter(SkogIO!=0)%>% + group_by(strata,region,FLATEID,SESONG)%>% + reframe(VMT=sum(VMT),AU_AREAL=sum(AU_AREAL))%>% + mutate(VMT_ha=VMT/AU_AREAL)%>% + filter(strata!="2") + +# Change var-of-interest to ind_ID 433 'ROS (Rogn-Osp-Selje)' (VOLUM m.b. pr ha, trees >10cm DBH, tree species Rogn-Osp-Selje) +fla_boot<-fla%>%mutate(var_of_interest=VMT_ha)%>% as.data.frame() +# Result ROS indicator (ind_ID 433) +Res_ind_433<-cbind("Rogn-Osp-Selje",beregninger.boot(fla_boot)) + +``` + +# Beregning av indikatoren Blåbær (ind_ID 23) + +```{r} +################################## +## Blåbær indicator (ind_ID 23) ## +################################## + +fla<-dta2%>%filter(!is.na(BLAABAER_GJSN))%>% + mutate_all(~replace(., is.na(.),0))%>% + filter(SkogIO!=0)%>% + group_by(strata,region,FLATEID,SESONG)%>% + reframe(BLAABAER_GJSN=weighted.mean(BLAABAER_GJSN,AU_AREAL), + AU_AREAL=sum(AU_AREAL))%>%filter(strata!="2") + +# Change var-of-interest to ind_ID 23 'Blåbær' +fla_boot<-fla%>%mutate(var_of_interest=BLAABAER_GJSN)%>% as.data.frame() +# Result 'Blåbær' indicator (ind_ID 23) +Res_ind_23<-cbind("Blåbær",beregninger.boot(fla_boot)) + +``` + +# Beregning av MiS-indikatorene (ind_ID 41, 60, 108, 183 og 207) + +```{r} +####################################################### +## MiS-indicators (ind_ID 41, 60, 108, 183, and 207) ## +####################################################### + +dat_MIS_prod<-MIS%>%filter(strata!="2") + +dat_MIS_prod<-dat_MIS_prod%>%mutate(AREALTYPE=as.numeric(AREALTYPE), + ANDEL_2DAA=case_when(ANDEL_2DAA==99~100, + TRUE~ANDEL_2DAA), + ANDEL_SKOG=case_when(ANDEL_SKOG==99~100, + TRUE~ANDEL_SKOG), + + areal=(ANDEL_2DAA/ANDEL_SKOG*100)) #*AU_AREAL/100 + +dat_MIS_prod<-merge(dat_MIS_prod,regioner,by=c("FLATEID"),all.x = T) +dat_MIS_prod<-dat_MIS_prod%>%mutate(region=Area.Name) + +dat_MIS_prod$ELEMENT_NAVN[dat_MIS_prod$ELEMENT_NR=="1"]<-"Stående død ved (MiS) – arealandel" # ind_ID 183 +dat_MIS_prod$ELEMENT_NAVN[dat_MIS_prod$ELEMENT_NR=="2"]<-"Liggende død ved (MiS) – arealandel" # ind_ID 108 +dat_MIS_prod$ELEMENT_NAVN[dat_MIS_prod$ELEMENT_NR=="4"]<-"Trær med hengelav (MiS)" # ind_ID 207 +dat_MIS_prod$ELEMENT_NAVN[dat_MIS_prod$ELEMENT_NR=="5"]<-"Eldre lauvsuksesjon (MiS)" # ind_ID 41 +dat_MIS_prod$ELEMENT_NAVN[dat_MIS_prod$ELEMENT_NR=="6"]<-"Gamle trær (MiS)" # ind_ID 60 + +dat_MIS_prod_wide<-dat_MIS_prod%>% + select(FLATEID,areal,AREALTYPE,region, + AU_AREAL,ELEMENT_NAVN,DISTRIKT, + LS_NETT,SESONG,AU_AREAL_SKOG)%>% + mutate_all(~replace(., is.na(.), 0))%>% + pivot_wider(names_from=ELEMENT_NAVN,values_from = areal) + +fla_boot_all<-dat_MIS_prod_wide%>% + mutate_all(~replace(., is.na(.),0))%>% + mutate(strata = case_when(DISTRIKT<7&LS_NETT=="3x3"~"1", + DISTRIKT<7&LS_NETT!="3x3"~"2", + DISTRIKT==7&LS_NETT=="3x3"~"3", + DISTRIKT==7&LS_NETT!="3x3"~"4", + TRUE~ "NA"))%>% + filter(AU_AREAL_SKOG>0)%>% + group_by(strata,region,FLATEID,SESONG)%>% + reframe(`Stående død ved (MiS) – arealandel`=sum(`Stående død ved (MiS) – arealandel`), + + `Liggende død ved (MiS) – arealandel`=sum(`Liggende død ved (MiS) – arealandel`), + + `Trær med hengelav (MiS)`=sum(`Trær med hengelav (MiS)`), + + `Eldre lauvsuksesjon (MiS)`=sum(`Eldre lauvsuksesjon (MiS)`), + + `Gamle trær (MiS)`=sum(`Gamle trær (MiS)`), + + AU_AREAL=sum(AU_AREAL), + + AU_AREAL_SKOG=sum(AU_AREAL_SKOG))%>% + as.data.frame() + +alle_lm<-c("Stående død ved (MiS) – arealandel", +"Liggende død ved (MiS) – arealandel", +"Trær med hengelav (MiS)", +"Eldre lauvsuksesjon (MiS)", +"Gamle trær (MiS)" +) +names_for_columns<-c("Indicator name","Area name", "wtd.avg", "CI.LL", "CI.UL","SESONG","n") + +for (i in alle_lm){ + +fla_boot_one<-fla_boot_all%>%mutate(var_of_interest=!!sym(i)) + +# Result for all MiS-indicators +res_one<-cbind(paste0(i),beregninger.boot(fla_boot_one)) +colnames(res_one)<-names_for_columns +if(i==alle_lm[1]){ +res_mis<-res_one +}else{ + res_mis<-rbind(res_mis,res_one) +} + +} + +``` + +# Samler resultatene + +```{r} + + colnames(Res_ind_419)<-names_for_columns + colnames(Res_ind_418)<-names_for_columns + colnames(Res_ind_434)<-names_for_columns + colnames(Res_ind_433)<-names_for_columns + colnames(Res_ind_23)<-names_for_columns + + + Res_alle<-rbind(Res_ind_419, + Res_ind_418, + Res_ind_434, + Res_ind_433, + Res_ind_23, + res_mis) + +Res_alle<-merge(Res_alle, Indicator_ID_long_format, by=c("Area name", "Indicator name")) + +``` + +# Henter data fra NI-databasen og lagrer data lokalt + +```{r} +# devtools::install_github("NINAnor/NIcalc", build_vignettes = T) + +library(NIcalc) + +NIcalc::getToken( + username = "", + password = "", + url = "https://www8.nina.no/NaturindeksNiCalc" +) + +all_indicators<-NIcalc::getIndicators() + +Res_alle<-Res_alle%>%mutate(SESONG="2024") +til_loop<-unique(Res_alle$Indicator_ID) + +for(i in til_loop){ + data<-NIcalc::getIndicatorValues(i) +data[["indicatorValues"]]<-merge(data[["indicatorValues"]], + + Res_alle[c("Indicator_ID","Area name","wtd.avg","CI.LL","CI.UL","SESONG")], + + by.x =c("indicatorId","areaName","yearName"), + by.y=c("Indicator_ID","Area name","SESONG"), + all.x=T)%>% + mutate( + + verdi=case_when(is.na(`wtd.avg`)~verdi, + TRUE~`wtd.avg`), + + nedre_Kvartil=case_when(is.na(`CI.LL`)~nedre_Kvartil, + TRUE~`CI.LL`), + + ovre_Kvartil=case_when(is.na(`CI.UL`)~ovre_Kvartil, + TRUE~`CI.UL`), + + datatypeId=case_when(yearName==unique(Res_alle$SESONG)~2, + TRUE~datatypeId), + + datatypeName=case_when(yearName==unique(Res_alle$SESONG)~"Overvåkingsdata", + TRUE~datatypeName))%>% + select(-c(wtd.avg,CI.LL,CI.UL)) + +assign(paste0("final_res_ind_",i),data) + +} + + +``` + +# Figurer for kvalitetssikring av resultatene + +```{r} +# Specifications for plot +COLORS <- c(current = "cornflowerblue") +widths <- c("2000"=0.2, "2010"=0.2, "2014"=0.2, "2019"= 0.2) + +## Make plot for all NI-regions of data already present in NI-db, i.e., change over time +for(i in til_loop){ + res_ind_temp<-get(paste0("final_res_ind_",i)) + +# Change value '2024' to 'current' last value on x-axes, and remove NA-values (i.e., '-1') and ref.values from plot +for_plot<-res_ind_temp[["indicatorValues"]]%>% + mutate(yearName=case_when(yearName=="2024"~"current", + TRUE~yearName))%>% + filter(yearName%in%c("2000", "2010", "2014", "2019","current"), + !is.na(verdi), + !verdi==-1, + !yearName=="Referanseverdi")%>% + data.frame + +# Order on x-axes +for_plot$yearName<-factor(for_plot$yearName, + levels = c("2000", "2010", "2014", "2019","current")) + +# Add value for current estimates, for visual interpretation and quality control + ggplot(for_plot,aes(x=yearName,y=verdi,group=areaName,colour = yearName))+ + geom_path(aes(linewidth = yearName))+ + geom_point()+ + scale_color_manual(values = COLORS)+ + scale_linewidth_manual(values = widths)+ + geom_errorbar(aes(ymin=nedre_Kvartil, ymax=ovre_Kvartil), width=.2, + position=position_dodge(0.05))+ + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ + facet_wrap(vars(areaName), scales = "free_y")+ + theme(legend.position = "none") + +# Save plot as .png +ggsave(filename = paste0("utv_indicator_",i,"_",for_plot$indicatorName,".png"), + scale = 1, + dpi = 1000, + width = 20, + height = 20, + units = "cm") +} + +``` + +# Lagrer nye beregninger i NI-db + +```{r} + +# Loop collecting values for all indicators +for(i in til_loop){ + res_ind_temp<-get(paste0("final_res_ind_",i)) + +## for .csv-files + if(i==til_loop[1]){ + res_og_gamle<-res_ind_temp[["indicatorValues"]]%>%filter(yearName%in%c(2000, 2010, 2014, 2019,as.numeric(unique(Res_alle$SESONG ))),!is.na(verdi),!yearName== +"Referanseverdi") + }else{ + res_og_gamle<-rbind(res_og_gamle,res_ind_temp[["indicatorValues"]]) + } + + +## Update/write new values to NI-db (to avoid unintentional overwriting in database, the next row is commented out) +# NIcalc::writeIndicatorValues(res_ind_temp) +} + +## Save result as .csv-files +write.csv2(Res_alle,file="samlet_resultat.csv") +write.csv2(res_og_gamle,file="samlet_resultat_og_gamle.csv") + + +```