Skip to content
Snippets Groups Projects
Commit 080b1e42 authored by Arvid Svensson's avatar Arvid Svensson
Browse files

Upload New File

parent 4e85b865
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 487 additions and 0 deletions
script_NI_LSK.Rmd 0 → 100644
+ 487
− 0
View file @ 080b1e42
---
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")
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment