Fjellrype
Klikk her for å se starten på utregningen av indiatorverdien som ble gjort av Markus F Israelsen.
På denne siden gjør jeg de siste stegene i analysen og viser resultatene.
Import
dat <- read_excel("../output/indicator_values/fjellrype.xlsx")
head(dat)## # A tibble: 6 x 12
## `2010` `2011` `2012` `2013` `2014` `2015` `2016` `2017` `2018` `2019` `2020`
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 2.94 1.37 1.98 1.97 1.67 1.08 2.26 1.81 1.71 0.919
## 2 1 0.706 2.27 1.50 3.71 2.13 1.27 1.32 1.18 1.15 1.41
## 3 1 1.31 0.819 3.07 2.16 1.71 0.876 1.85 1.67 1.27 2.51
## 4 1 0.973 1.47 1.86 5.06 1.30 3.16 1.99 1.05 1.87 1.39
## 5 1 2.63 1.20 1.37 1.86 1.66 0.875 4.79 1.07 0.893 0.691
## 6 1 1.50 2.19 1.14 2.39 1.58 2.38 1.66 1.99 1.31 0.678
## # ... with 1 more variable: Region <chr>Dette er bootstrappede indikatorverdier, skalerte mot indeksåret 2010.
La oss få dette over i langt format.
setDT(dat)
dat <- melt(dat,
id.vars = 'Region')
head(dat)## Region variable value
## 1: nord 2010 1
## 2: nord 2010 1
## 3: nord 2010 1
## 4: nord 2010 1
## 5: nord 2010 1
## 6: nord 2010 1names(dat) <- c("reg", "year", "value")
table(dat$reg, dat$year)##
## 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020
## midt 10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 10000
## nord 10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 10000
## sør 10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 10000
## vest 10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 10000
## øst 10000 10000 10000 10000 10000 10000 10000 10000 10000 10000 1000010 000 verdier per kombinasjon.
Gir regionene nytt navn.
unique(dat$reg)## [1] "nord" "midt" "øst" "vest" "sør"dat$reg <- plyr::revalue(dat$reg,
c("nord"="N",
"midt"="C",
"øst" ="E",
"vest"="W",
"sør"="S"))plot(dat$value)
Det er mange høye verdier. En verdi på 2 tilsier 2x så mange ryper som i 2010.
dat$year <- as.numeric(as.character(dat$year))
ggplot(data=dat, aes(x = year, y = value, group= reg, colour = reg))+
geom_smooth()## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Det har vært en bedring i bestandedn i de siste årene. Østlandet har noen store svingninger.
Nå må vi skalere dataene etter en referansetilstand, som vi henter fra Naturindeksen. Der er det oppgitt hvor stor prosentanndel av referasnetilstanden som var realisert i 2010.
Henter referansetilstand
Fyll inn ditt eget passord og brukernavn
myUser <- "anders.kolstad@nina.no"
myPwd <- "" # hemmelig passordImporterer data fra NI-databasen
ref <- NIcalc::importDatasetApi(
username = myUser,
password = myPwd,
eco = "Fjell",
indic = "Fjellrype",
year = c(2010))
saveRDS(ref, "../data/fjellrypeNIexport.rds")ref <- readRDS("../data/fjellrypeNIexport.rds")ref2 <- ref$referenceValues$referenceValues
ref2 <- select(ref2, ICunitName, expectedValue)Det er usikkerhet knyttet til referanseverdiene, men her behandler jeg dem som faste.
comp <- ref$indicatorObservations$indicatorValues
comp <- select(comp, ICunitName, expectedValue)
ref2$comp <- comp$expectedValue[match(ref2$ICunitName, comp$ICunitName)]
head(ref2)## # A tibble: 6 x 3
## ICunitName expectedValue comp
## <chr> <dbl> <dbl>
## 1 Hedmark 49.3 17.2
## 2 Oppland 114. 35.7
## 3 Buskerud 42.3 14.8
## 4 Telemark 38.6 19
## 5 Aust-Agder 19.7 12.7
## 6 Vest-Agder 17.7 13Så må jeg summere bestandsestimate for hver region
ref2$reg <- NA
ref2$reg[ref2$ICunitName == "Hedmark" |
ref2$ICunitName == "Oppland"|
ref2$ICunitName == "Buskerud"] <- "E"
ref2$reg[ref2$ICunitName == "Telemark" |
ref2$ICunitName == "Aust-Agder"|
ref2$ICunitName == "Vest-Agder"] <- "S"
ref2$reg[ref2$ICunitName == "Rogaland" |
ref2$ICunitName == "Hordaland"|
ref2$ICunitName == "Sogn og Fjordane"] <- "W"
ref2$reg[ref2$ICunitName == "Møre og Romsdal" |
ref2$ICunitName == "Sør-Trøndelag"|
ref2$ICunitName == "Nord-Trøndelag"] <- "C"
ref2$reg[ref2$ICunitName == "Nordland" |
ref2$ICunitName == "Troms"|
ref2$ICunitName == "Finnmark"] <- "N"ref3 <- aggregate(data=ref2,
expectedValue~reg,
FUN = sum)
temp <- aggregate(data=ref2,
comp~reg,
FUN = sum)
ref3$comp <- temp[,2]
ref3$realisert <- ref3$comp/ref3$expectedValueSkalerer TRIMverdiene mot referansetilstand
dat$factor <- ref3$realisert[match(dat$reg, ref3$reg)]
dat$val <- dat$value*dat$factor
# trunkerer
dat$val[dat$val >1] <- 1Veiing
Regner ut nasjonale indikatorverider som et veid (etter total fjellareal) gjennomsnitt av vedien i regionene. Veiingen gjøres ved å samle større antall verdier for de regionene som har mest fjell i seg.
wgt <- readRDS("../data/fjellareal.rds")
wgt$Fjellareal2 <- wgt$Fjellareal/max(wgt$Fjellareal)
wgt$reg <- c("N", "C", "E", "W", "S")myYears <- 2010:2020
datNorge <- data.frame(
reg = rep("Norge", 110000),
year = rep(2010:2020, each = 10000),
val = NA
)
for(n in myYears){
temp <- dat[dat$year == n,]
temp2 <- c(
sample(temp$val[temp$reg == "N"], wgt$Fjellareal2[wgt$reg == "N"]*10000, replace =T),
sample(temp$val[temp$reg == "E"], wgt$Fjellareal2[wgt$reg == "E"]*10000, replace =T),
sample(temp$val[temp$reg == "W"], wgt$Fjellareal2[wgt$reg == "W"]*10000, replace =T),
sample(temp$val[temp$reg == "S"], wgt$Fjellareal2[wgt$reg == "S"]*10000, replace =T),
sample(temp$val[temp$reg == "C"], wgt$Fjellareal2[wgt$reg == "C"]*10000, replace =T)
)
datNorge$val[datNorge$year == n] <- sample(temp2, 10000, replace = T)
}dat <- rbind(select(dat, reg, year, val), datNorge)Plotting
eval(parse("indicator_plots2.R", encoding="UTF-8"))png("../output/indicatorPlots/skalert/fjellrype.png",
units="in", width=12, height=7, res=300)
par(mfrow=c(1,1), mar=c(4.5,
5.5,
0,
2))
indicator_plot2(dataset = dat,
yAxisTitle = "Fjellrypebestand skalert mot referanseverdi",
lowYlimit = 0,
upperYlimit = 1.2,
yStep = .2,
minyear = 2009,
maxyear = 2021,
colours = c("#FFB25B", "#2DCCD3", "#004F71", "#7A9A01", "#93328E", "dark grey"),
legendPosition = "top",
legendInset = 0,
move = 0.05,
horizontal = T,
legendTextSize = 1.25)
dev.off()Det blir så stor variasjon her at det er vanskelig å lese figuren: 
finalTbl <- aggregate(data=dat,
val~year+reg,
FUN= function(x) round(
quantile(x, c(.025, .5, .975)), 2))
finalTbl <- do.call(data.frame, finalTbl)
names(finalTbl) <- c("year", "reg", "low", "med", "upp")
DT::datatable(
finalTbl,
extensions = "FixedColumns",
options = list(
scrollX = TRUE,
scrollY=T,
pageLength = 10
))Her er en annen figur som ikke viser variasjonen.
ggplot(data=finalTbl, aes(x = year, y = med, group = reg, linetype=reg, colour=reg))+
geom_line(size=2)+
theme_bw(base_size = 20)+
ylab("Fjellrypebestand\nskalert mor referanseverdi")+
scale_x_continuous(breaks=c(2010, 2015, 2020))+
xlab("")
Og en versjon til
gg <- ggplot(data = finalTbl, aes(x = year, y = med))+
geom_line(size=1.2, colour="grey30")+
geom_errorbar(aes(ymax=upp, ymin=low), size = .8,
colour="grey30",
alpha=0.6,
width=0.5)+
theme_bw(base_size = 20)+
theme(
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)
)+
scale_x_continuous(breaks = c(2010, 2012, 2014, 2016, 2018, 2020))+
ylab("Fjellrypebestand\nskalert mot referanseverdi")+
xlab("")+
facet_wrap(.~reg)png("../output/indicatorPlots/fjellrype_facet.png",
units="in", width=12, height=7, res=300)
gg
dev.off()
Som alltid, errorbars er 95%CI.
Eksporter csv
write.csv(dat, "../output/indicator_values/fjellrype.csv", row.names = F)