Plotting
Her viser jeg hvordan vi aggreggerer indikatorverdiene til en verdi for økologisk tilstandsverdi og plotter denne. Det meste er scriptet er skrevet av Somon Jacobsson.
# Number of simulations
nsim <- 10000 # should be run with 10 000 as the end
# Short names for ecosystem characteristics
characteristics <- c("ppr", "isf", "bmb", "fgw", "bio", "lsp", "abf")Step 1:
indicators, characteristics, years and period summaries
Specify reporting year key to periods
periods <- data.frame(period = c(1,1,1,1,1,1,1,1,
2,2,2,2,2,
3,3,3,3,3,
4,4,4,4,4,
5,5,5,5,5,
6,6,6,6,6),
year = c(1988, 1989, 1990, 1991, 1992, 1993, 1994, 1995,
1996, 1997, 1998, 1999, 2000,
2001, 2002, 2003, 2004, 2005,
2006, 2007, 2008, 2009, 2010,
2011, 2012, 2013, 2014, 2015,
2016, 2017, 2018, 2019, 2020))Step 2:
load and adjust indicator data
List files
List files with bootstrapped indicator values
list.files(path = "../output/indicator_values/",
pattern=".csv")## [1] "alien.csv" "breareal.csv" "ell_L_high.csv"
## [4] "ell_L_low.csv" "ell_N_high.csv" "ell_N_low.csv"
## [7] "ellenberg_alien.csv" "fjellindeks.csv" "fjellrev.csv"
## [10] "fjellrype.csv" "fragmentering.csv" "heat_req.csv"
## [13] "inon.csv" "jerv.csv" "kongeorn.csv"
## [16] "lirype.csv" "ndvi_lower.csv" "ndvi_upper.csv"
## [19] "rein.csv" "smågnagere.csv" "snodekke.csv"
## [22] "snodybde.csv" "varmekrav.csv" "vinterregn.csv"Load and adjust data headings etc.
manual steps required due to different data structures
Alien
alien <- read.csv("../output/indicator_values/alien.csv", header=T)
alien$period <-
periods$period[match(alien$year, periods$year)]
alien$X <- "alien"Breareal
breareal <- read.csv("../output/indicator_values/breareal.csv", header=T)
breareal$period <-
periods$period[match(breareal$year, periods$year)]
breareal$X <- "breareal"Ellenberg L
Denne indikatoren er tosidig (se her).
EllLHigh <- read.csv("../output/indicator_values/ell_L_high.csv", header=T)
EllLLow <- read.csv("../output/indicator_values/ell_L_low.csv", header=T)
EllLHigh$period <-
periods$period[match(EllLHigh$year, periods$year)]
EllLLow$period <-
periods$period[match(EllLLow$year, periods$year)]
EllLHigh$X <- "EllLHigh"
EllLLow$X <- "EllLLow"Ellenberg N
Denne indikatoren er også tosidig (se her).
EllNHigh <- read.csv("../output/indicator_values/ell_N_high.csv", header=T)
EllNLow <- read.csv("../output/indicator_values/ell_N_low.csv", header=T)
EllNHigh$period <-
periods$period[match(EllNHigh$year, periods$year)]
EllNLow$period <-
periods$period[match(EllNLow$year, periods$year)]
EllNHigh$X <- "EllNHigh"
EllNLow$X <- "EllNLow"Fjellindeks
fjellindeks <- read.csv("../output/indicator_values/fjellindeks.csv", header=T)
fjellindeks$period <-
periods$period[match(fjellindeks$year, periods$year)]
fjellindeks$X <- "fjellindeks"Fjellrev
fjellrev <- read.csv("../output/indicator_values/fjellrev.csv", header=T)
fjellrev$period <-
periods$period[match(fjellrev$year, periods$year)]
fjellrev$X <- "fjellrev"
table(fjellrev$reg, fjellrev$period)##
## 1 2 4 5 6
## C 1000 1000 1000 1000 1000
## E 1000 1000 1000 1000 1000
## N 1000 1000 1000 1000 1000
## Norge 1000 1000 1000 1000 1000
## S 1000 1000 1000 1000 1000
## W 1000 1000 1000 1000 1000Fjellrype
Her må jeg skjekke litt.
fjellrype <- read.csv("../output/indicator_values/fjellrype.csv", header=T)
fjellrype$X <- "fjellrype"
fjellrype$period <- periods$period[match(fjellrype$year, periods$year)]Fragmentering
Her må jeg skjekke litt.
fragmentering <- read.csv("../output/indicator_values/fragmentering.csv", header=T)
fragmentering$year <- 2020
fragmentering$period <- periods$period[match(fragmentering$year, periods$year)]INON
INON har ingen usikkerhet knyttet til seg.
inon <- read.csv("../output/indicator_values/inon.csv", header=T)
inon$X <- "INON"
inon$period <- periods$period[match(inon$year, periods$year)]Er INON og Konnektivitetsindeksen korrelert?
plot(fragmentering$val,
inon$val[inon$period==6],
ylim = c(0,1), xlim=c(0,1),
ylab = "INON", xlab = "Konnektivitet",
col = "black", pch = 21, bg="grey", cex=2, lwd=2)
abline(a=0, b=1)
text(0.6, 0.2, paste("Pearson's r = ",
round(cor(fragmentering$val,
inon$val[inon$period==6]), 2)))
De er ganske korrelerte ja!
Jerv
jerv <- read.csv("../output/indicator_values/jerv.csv", header=T)
jerv$X <- "jerv"
jerv$period <-
periods$period[match(jerv$year, periods$year)]
table(jerv$reg, jerv$period)##
## 1 2 4 5 6
## C 1000 1000 1000 1000 1000
## E 1000 1000 1000 1000 1000
## N 1000 1000 1000 1000 1000
## Norge 1000 1000 1000 1000 1000
## S 1000 1000 1000 1000 1000
## W 1000 1000 1000 1000 1000Kongeørn
kongeorn <- read.csv("../output/indicator_values/kongeorn.csv", header=T)
kongeorn$period <-
periods$period[match(kongeorn$year, periods$year)]
kongeorn$X <- "kongeørn"Lirype
lirype <- read.csv("../output/indicator_values/lirype.csv", header=T)
lirype$year <- 2020
lirype$period <-
periods$period[match(lirype$year, periods$year)]
lirype$X <- "lirype"NDVI
ndvi_upper <- read.csv("../output/indicator_values/ndvi_upper.csv", header=T)
ndvi_upper$year <- 2020
ndvi_upper$period <-
periods$period[match(ndvi_upper$year, periods$year)]
ndvi_upper <- dplyr::select(ndvi_upper, -X.1)
table(ndvi_upper$reg)##
## C E N Norge S W
## 10000 10000 10000 10000 10000 10000ndvi_lower <- read.csv("../output/indicator_values/ndvi_lower.csv", header=T)
ndvi_lower$year <- 2020
ndvi_lower$period <-
periods$period[match(ndvi_lower$year, periods$year)]
ndvi_lower <- dplyr::select(ndvi_lower, -X.1)
table(ndvi_lower$reg)##
## C E N Norge S W
## 10000 10000 10000 10000 10000 10000Rein
rein <- read.csv("../output/indicator_values/rein.csv", header=T)
rein$period <-
periods$period[match(rein$year, periods$year)]
table(rein$reg, rein$period)##
## 6
## C 500
## E 500
## N 500
## Norge 1000
## S 500
## W 500Snødekke
snodekke <- read.csv("../output/indicator_values/snodekke.csv", header=T)
snodekke$X <- "snodekke"
snodekke$period <-
periods$period[match(snodekke$year, periods$year)]
table(snodekke$reg, snodekke$period)##
## 6
## C 10000
## E 10000
## N 10000
## Norge 10000
## S 10000
## W 10000Snødybde
snodybde <- read.csv("../output/indicator_values/snodybde.csv", header=T)
snodybde$period <-
periods$period[match(snodybde$year, periods$year)]
table(snodybde$reg, snodybde$period)##
## 6
## C 10000
## E 10000
## N 10000
## Norge 10000
## S 10000
## W 10000Smågnagere
smaagnagere <- read.csv("../output/indicator_values/smågnagere.csv", header=T)
smaagnagere$X <- "smaagnagere"
smaagnagere$period <-
periods$period[match(smaagnagere$year, periods$year)]
table(smaagnagere$reg, smaagnagere$period)##
## 1 2 4 5 6
## C 10000 10000 10000 10000 10000
## E 10000 10000 10000 10000 10000
## N 10000 10000 10000 10000 10000
## Norge 10000 10000 10000 10000 10000
## S 10000 10000 10000 10000 10000
## W 10000 10000 10000 10000 10000Vegetasjonen varmekrav
Fjernes?
varmekrav <- read.csv("../output/indicator_values/varmekrav.csv", header=T)
varmekrav <- dplyr::select(varmekrav, -X.1)
varmekrav$X <- "varmekrav"
varmekrav$period <-
periods$period[match(varmekrav$year, periods$year)]Vinterregn
vinterregn <- read.csv("../output/indicator_values/vinterregn.csv", header=T, encoding = "UTF-16")
vinterregn$X <- "vinterregn"
vinterregn$period <-
periods$period[match(vinterregn$year, periods$year)]
table(vinterregn$reg, vinterregn$period)##
## 6
## C 10000
## E 10000
## N 10000
## Norge 10000
## S 10000
## W 10000Step 3: Combine datasets
# Periods to use
periods_use <- 6Subset additional data
alien <- alien[alien$period%in%periods_use,]
breareal <- breareal[breareal$period%in%periods_use,]
EllLHigh <- EllLHigh[EllLHigh$period%in%periods_use,]
EllLLow <- EllLLow[EllLLow$period%in%periods_use,]
EllNHigh <- EllNHigh[EllNHigh$period%in%periods_use,]
EllNLow <- EllNLow[EllNLow$period%in%periods_use,]
fjellindeks <- fjellindeks[fjellindeks$period%in%periods_use,]
fjellrev <- fjellrev[fjellrev$period%in%periods_use,]
fjellrype <- fjellrype[fjellrype$period%in%periods_use,]
fragmentering <- fragmentering[fragmentering$period%in%periods_use,]
inon <- inon[inon$period%in%periods_use,]
jerv <- jerv[jerv$period%in%periods_use,]
kongeorn <- kongeorn[kongeorn$period%in%periods_use,]
lirype <- lirype[lirype$period%in%periods_use,]
ndvi_upper <- ndvi_upper[ndvi_upper$period%in%periods_use,]
ndvi_lower <- ndvi_lower[ndvi_lower$period%in%periods_use,]
rein <- rein[rein$period%in%periods_use,]
smaagnagere <- smaagnagere[smaagnagere$period%in%periods_use,]
snodekke <- snodekke[snodekke$period%in%periods_use,]
snodybde <- snodybde[snodybde$period%in%periods_use,]
varmekrav <- varmekrav[varmekrav$period%in%periods_use,]
vinterregn <- vinterregn[vinterregn$period%in%periods_use,]Order data
regOrder = c("C","N","W", "S", "E","Norge")
#perOrder = c(2,4,5,6) # not relevantalien <- alien %>%
arrange(match(reg, regOrder), desc(period)) %>%
select(X, reg, period, val)
breareal <- breareal %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
EllLHigh <- EllLHigh %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
EllLLow <- EllLLow %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
EllNHigh <- EllNHigh %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
EllNLow <- EllNLow %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
fjellindeks <- fjellindeks %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
fjellrev <- fjellrev %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
fjellrype <- fjellrype %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
fragmentering <- fragmentering %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
inon <- inon %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
jerv <- jerv %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
kongeorn <- kongeorn %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
lirype <- lirype %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
ndvi_upper <- ndvi_upper %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
ndvi_lower <- ndvi_lower %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
rein <- rein %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
snodekke <- snodekke %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
snodybde <- snodybde %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
smaagnagere <- smaagnagere %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
varmekrav <- varmekrav %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)
vinterregn <- vinterregn %>%
arrange(match(reg, regOrder), desc(period))%>%
select(X, reg, period, val)Add additional data to allind data
allind <- rbind(
alien,
breareal,
EllLHigh,
EllLLow,
EllNHigh,
EllNLow,
fjellindeks,
fjellrev,
fjellrype,
fragmentering,
inon,
jerv,
kongeorn,
lirype,
ndvi_upper,
ndvi_lower,
rein,
snodekke,
snodybde,
smaagnagere,
varmekrav,
vinterregn)
table(allind$X)##
## alien breareal EllLHigh EllLLow EllNHigh
## 60000 60000 60000 60000 60000
## EllNLow fjellindeks fjellrev fjellrype fragmentering
## 60000 60000 6000 300000 6
## INON jerv kongeørn lirype NDVI_lower
## 6 6000 60000 160000 60000
## NDVI_upper rein smaagnagere snodekke snodybde
## 60000 3500 60000 60000 60000
## varmekrav vinterregn
## 6000 60000The number of re-samplings differs.
Save allind file
write.csv(allind, "../output/allind_temp.csv")Weights
NDVI øvre
ndvihigh <- read.csv("../output/indicator_values/ndvi_upper.csv", header=T)
ndvilow <- read.csv("../output/indicator_values/ndvi_lower.csv", header=T)
ndviW <- select(ndvihigh, reg, high = nCells)
ndviW$low <- ndvilow$nCells
ndviW$tot <- ndviW$high+ndviW$low
ndviW$wHigh <- ndviW$high/ndviW$tot
ndviW$wLow <- ndviW$low/ndviW$tot
ndviW <- select(ndviW, reg, wLow, wHigh)
head(ndviW)## reg wLow wHigh
## 1 N 0.1154746 0.8845254
## 2 N 0.1154746 0.8845254
## 3 N 0.1154746 0.8845254
## 4 N 0.1154746 0.8845254
## 5 N 0.1154746 0.8845254
## 6 N 0.1154746 0.8845254Ellenberg L
Lhigh <- read.csv("../output/indicator_values/ell_L_high.csv", header=T)
Llow <- read.csv("../output/indicator_values/ell_L_low.csv", header=T)
Lhigh <- aggregate(data = Lhigh,
plotsL2_tot~reg,
FUN = mean) # a bit dangerous to use mean here, but I've cheched that all values are the same
Llow <- aggregate(data = Llow,
plotsL1_tot~reg,
FUN = mean)
L <- cbind(Lhigh, Llow$plotsL1_tot)
L$total <- L$plotsL2_tot+L$`Llow$plotsL1_tot`
L$wLow <- L$`Llow$plotsL1_tot`/L$total
L$wHigh <- L$plotsL2_tot/L$total
L <- select(L, reg, wLow, wHigh)
L## reg wLow wHigh
## 1 C 0.1459075 0.8540925
## 2 E 0.3375315 0.6624685
## 3 N 0.1906780 0.8093220
## 4 Norge 0.2085308 0.7914692
## 5 S 0.1714286 0.8285714
## 6 W 0.1400966 0.8599034Ellenberg N
Nhigh <- read.csv("../output/indicator_values/ell_N_high.csv", header=T)
Nlow <- read.csv("../output/indicator_values/ell_N_low.csv", header=T)
Nhigh <- aggregate(data = Nhigh,
plotsN2_tot~reg,
FUN = mean) # a bit dangerous to use mean here, but I've cheched that all values are the same
Nlow <- aggregate(data = Nlow,
plotsN1_tot~reg,
FUN = mean)
N <- cbind(Nhigh, Nlow$plotsN1_tot)
N$total <- N$plotsN2_tot+N$`Nlow$plotsN1_tot`
N$wLow <- N$`Nlow$plotsN1_tot`/N$total
N$wHigh <- N$plotsN2_tot/N$total
N <- select(N, reg, wLow, wHigh)
N## reg wLow wHigh
## 1 C 0.6441281 0.3558719
## 2 E 0.7229219 0.2770781
## 3 N 0.7235169 0.2764831
## 4 Norge 0.7040548 0.2959452
## 5 S 0.7428571 0.2571429
## 6 W 0.6473430 0.3526570Weigths table
Creating a data frame for the weights and setting all weights equal to 1 as defult
weights <- aggregate(data = allind,
val~X+reg,
FUN = length)
weights$val <- 1Then changing the weight for the two-sided indicators
weights$val[weights$X=="EllLHigh"] <- L$wHigh[match(weights$reg[weights$X=="EllLHigh"], L$reg)]
weights$val[weights$X=="EllLLow"] <- L$wLow[match(weights$reg[weights$X=="EllLLow"], L$reg)]
weights$val[weights$X=="EllNHigh"] <- N$wHigh[match(weights$reg[weights$X=="EllNHigh"], L$reg)]
weights$val[weights$X=="EllNLow"] <- N$wLow[match(weights$reg[weights$X=="EllNLow"], L$reg)]
weights$val[weights$X=="NDVI_upper"] <- ndviW$wHigh[match(weights$reg[weights$X=="NDVI_upper"], ndviW$reg)]
weights$val[weights$X=="NDVI_lower"] <- ndviW$wLow[match(weights$reg[weights$X=="NDVI_lower"], ndviW$reg)]Write weights file
write.csv(weights, "../output/weights_temp.csv")Step 4: Fill indicator <-> characteristics matrix
Empty indicator <-> character matrix
# Specify total number of indicators (incl. double counting of two-sided)
(nind <- length(unique(allind$X)))## [1] 22ind_char <- data.frame(matrix(nrow=nind, ncol=8))
names(ind_char) <- c("IND", characteristics)Start filling indicator <-> characteristic matrix
ind_char$IND <- unique(allind$X)# Primary productivity
ind_char$ppr <- ifelse(
ind_char$IND=="NDVI_upper"|
ind_char$IND=="NDVI_lower",
1,0)
# Biomass distribution among trophic levels
ind_char$bmb <- ifelse(
ind_char$IND=="NDVI_upper"|
ind_char$IND=="NDVI_lower"|
ind_char$IND=="fjellrev"|
ind_char$IND=="jerv"|
ind_char$IND=="rein"|
ind_char$IND=="smaagnagere"|
ind_char$IND=="lirype"|
ind_char$IND=="fjellrype"|
ind_char$IND=="kongeørn"
,
1,0)
# Functional groups within trophic levels
ind_char$fgw <- ifelse(
ind_char$IND=="fjellrev"|
ind_char$IND=="jerv"|
ind_char$IND=="rein"|
ind_char$IND=="smaagnagere"|
ind_char$IND=="lirype"|
ind_char$IND=="fjellrype"|
ind_char$IND=="kongeørn"
,
1,0)
# Structurally important species and biophysical structures
ind_char$isf <- ifelse(
ind_char$IND=="smaagnagere"|
ind_char$IND=="rein"|
ind_char$IND=="alien"
,
1,0)
# Biodiversity
ind_char$bio <- ifelse(
ind_char$IND=="fjellindeks"
,
1,0)
# Landscape ecological patterns
ind_char$lsp <- ifelse(
ind_char$IND=="INON"|
ind_char$IND=="fragmentering"
,
1,0)
# Abiotic factors
ind_char$abf <- ifelse(
ind_char$IND=="EllNHigh"|
ind_char$IND=="EllNLow"|
ind_char$IND=="EllLHigh"|
ind_char$IND=="EllLLow"|
ind_char$IND=="breareal"|
ind_char$IND=="snodybde"|
ind_char$IND=="snodekke"|
ind_char$IND=="vinterregn"|
ind_char$IND=="varmekrav"
,
1,0)Write indicator-characteristic matrix
write.csv(ind_char, "../output/ind_char.csv")Step 5: Indicator quantiles per region and period
#unique(allind$X[allind$reg=="A"])
periods <- unlist(unique(allind$period))
regions <- unlist(unique(allind$reg))
per_reg <- length(periods)*length(regions)all_Q <- data.frame()
for(n in unique(allind$X)){
temp <- allind[allind$X == n,]
td <- aggregate(data = temp,
val~X+reg+period,
FUN = function(x) c(quantile(x, c(0.025, .5, 0.975))))
td <- do.call(data.frame, td)
names(td) <- c("ind", "reg", "per", "low", "med", "upp")
all_Q <- rbind(all_Q, td)
}Step 6: Aggregated characteristic calculations
agg_ind <- data.frame()
w_boot <- data.frame()for (n in 1:length(characteristics)){
# Select by ind_char data frame
temp_char <- cbind(ind_char$IND, ind_char[names(ind_char)%in%characteristics[n]])
temp_char <- temp_char[temp_char[,2]==1,]
temp_ind <- allind[allind$X %in% temp_char$`ind_char$IND`,]
# Here Simon had a for loop for subsetting period, but I only ise one period
for (k in 1:length(regions)){
temp_ind3 <- subset(temp_ind, reg==regions[k])
# Simon use a wide format, but thats not so good when the number of rows differ between indicators
# weights
temp_weights <- weights[weights$reg==regions[k],]
temp_weights <- temp_weights[temp_weights$X %in% unique(temp_ind3$X),]
#temp_weights2 <- decostand(temp_weights$val, method="max", margin =2)[,1]
# check?
#(names(temp_ind3)==names(temp_weights))
# Empty sample values vector
temp_x <- NULL
temp_x2 <- NULL
for (i in 1:length(unique(temp_ind3$X))) {
print(paste(
characteristics[n],
regions[k],
unique(temp_ind3$X)[i]))
# Sample from indicator i
start <- length(temp_x)
temp_x[(start+1):(start+nsim)] <-
tryCatch(
sample(
temp_ind3$val[temp_ind3$X == unique(temp_ind3$X)[i]],
nsim, replace = T),
error=function(e){})
temp_x2[(start+1):(start+nsim)] <- unique(temp_ind3$X)[i]
}
temp_x3 <- as.data.frame(cbind(val=temp_x,
ind =temp_x2))
setDT(temp_x3)
temp_x3$sim <- rep(1:nsim, times=length(unique(temp_ind3$X)))
temp_x3$val <- as.numeric(temp_x3$val)
temp_x3 <- data.table::dcast(temp_x3,
sim~ind,
value.var="val")
mat <- as.matrix(temp_x3[,-1])
x_boot <- rowWeightedMeans(mat, temp_weights$val)
# Info data
chr <- characteristics[n]
reg <- regions[k]
#per <- periods[j]
w_boot <- rbind(w_boot, data.frame(chr, reg, x_boot))
# Quantiles
temp_Q <- tryCatch(quantile(x_boot,c(0.025, 0.5, 0.975)), error=function(e){})
# Compile data
agg_ind <- rbind(agg_ind, data.frame(chr, reg,
ifelse(is.null(temp_Q), NA, temp_Q[1]),
ifelse(is.null(temp_Q), NA, temp_Q[2]),
ifelse(is.null(temp_Q), NA, temp_Q[3])))
}
}
# Rename columns
names(agg_ind) <- c("ind", "reg", "low", "med", "upp")
names(w_boot)[3] <- c("val")Step 7: Aggregated total calculations
regions <- unlist(unique(allind$reg))
char_w <- unlist(unique(w_boot$chr))
tot_ind <- data.frame()
tot_ind_w <- data.frame()
#per_reg <- length(periods)*length(regions)
weights_all <- data.frame()
# Subset period
temp_ind2 <- allind
temp_chr2 <- w_boot
for (k in 1:length(regions)){
# Subset region
temp_ind3 <- subset(temp_ind2, reg==regions[k])
temp_chr3 <- subset(temp_chr2, reg==regions[k])
# temp_ind3 <- temp_ind3[,-(1:3)]
# weights
temp_weights <- weights[weights$reg==regions[k],]
#temp_weights <- temp_weights[temp_weights$X %in% unique(temp_ind3$X),]
#temp_weights2 <- decostand(temp_weights$val, method="max", margin =2)[,1]
#
# Empty sample mean vector
temp_x_boot <- NULL
temp_chr_boot <- NULL
# Empty sample values vector
temp_x <- NULL
temp_x2 <- NULL
for (i in 1:length(unique(temp_ind3$X))) {
# Sample from indicator i
print(paste(regions[k],
unique(temp_ind3$X)[i]))
start <- length(temp_x)
temp_x[(start+1):(start+nsim)] <-
tryCatch(
sample(
temp_ind3$val[temp_ind3$X == unique(temp_ind3$X)[i]],
nsim, replace = T),
error=function(e){})
temp_x2[(start+1):(start+nsim)] <- unique(temp_ind3$X)[i]
}
temp_x3 <- as.data.frame(cbind(val=temp_x,
ind =temp_x2))
setDT(temp_x3)
temp_x3$sim <- rep(1:nsim, times=length(unique(temp_ind3$X)))
temp_x3$val <- as.numeric(temp_x3$val)
temp_x3 <- data.table::dcast(temp_x3,
sim~ind,
value.var="val")
mat <- as.matrix(temp_x3[,-1])
x_boot <- rowWeightedMeans(mat, temp_weights$val)
# Empty sample values vector
temp_chr <- NULL
temp_chrX <- NULL
for (n in 1:length(char_w)) {
temp_chr4 <- subset(temp_chr3, chr==char_w[n])
start2 <- length(temp_chr)
temp_chr[(start2+1):(start2+nsim)] <-
tryCatch(
sample(
temp_chr4$val,
nsim, replace = T),
error=function(e){})
temp_chrX[(start2+1):(start2+nsim)] <- char_w[n]
}
temp_chr3 <- as.data.frame(cbind(val=temp_chr,
chr =temp_chrX))
setDT(temp_chr3)
temp_chr3$sim <- rep(1:nsim, times=length(char_w))
temp_chr3$val <- as.numeric(temp_chr3$val)
temp_chr3 <- data.table::dcast(temp_chr3,
sim~chr,
value.var="val")
temp_chr_boot <- rowMeans(temp_chr3[,-1])
# Characteristic mean (unweighted)
#temp_chr_boot[m] <- mean(na.omit(temp_chr))
# Info data
tot <- "Total"
totw <- "Total_w"
reg <- regions[k]
#per <- periods[j]
weights_all <- rbind(weights_all, data.frame(reg, temp_weights))
# Quantiles
temp_Q <- tryCatch(quantile(x_boot,c(0.025, 0.5, 0.975)), error=function(e){})
temp_Qw <- tryCatch(quantile(temp_chr_boot,c(0.025, 0.5, 0.975)), error=function(e){})
# Compile data
tot_ind <- rbind(tot_ind, data.frame(tot, reg,
ifelse(is.null(temp_Q), NA, temp_Q[1]),
ifelse(is.null(temp_Q), NA, temp_Q[2]),
ifelse(is.null(temp_Q), NA, temp_Q[3])))
tot_ind_w <- rbind(tot_ind_w, data.frame(totw, reg,
ifelse(is.null(temp_Qw), NA, temp_Qw[1]),
ifelse(is.null(temp_Qw), NA, temp_Qw[2]),
ifelse(is.null(temp_Qw), NA, temp_Qw[3])))
}
# Rename columns
names(tot_ind) <- c("ind", "reg", "low", "med", "upp")
names(tot_ind_w) <- c("ind", "reg", "low", "med", "upp")Step 8: Rename indicator labels key
temp_names <- data.frame(c(unique(allind$X), names(ind_char[,-1]),
"Total", "Total_w"))
names(temp_names) <- "work"
temp_names$plot <- c("Fravær av fremmede arter",
"Areal av isbreer",
"Ellenberg L (øvre)",
"Ellenberg L (nedre)",
"Ellenberg N (øvre)",
"Ellenberg N (nedre)",
"Naturindeks for fjell (mod.)",
"Fjellrev",
"Fjellrype",
"Konnektivitet", # Fragmentering
"Areal uten tekniske inngrep",
"Jerv",
"Kongeørn",
"Lirype",
"NDVI (øvre)",
"NDVI (nedre)",
"Rein",
"Snødekkets varighet",
"Snødybde",
"Smågnagere",
"Vegetasjonens varmekrav",
"Vinterregn",
"Primærproduksjon",
"Funksjonelt viktige arter og biofysiske strukturer",
"Fordeling av biomasse i ulike trofiske nivå",
"Funksjonelle grupper innen trofiske nivåer",
"Biologisk mangfold",
"Landskapsøkologiske mønstre",
"Abiotiske forhold",
"Total",
"Total (vektet pr. egenskap)")Engelske navn
temp_names2 <- data.frame(c(unique(allind$X), names(ind_char[,-1]),
"Total", "Total_w"))
names(temp_names2) <- "work"
temp_names2$plot <- c("Absence of alien species",
"Area of glaciers",
"Ellenberg L (upper)",
"Ellenberg L (lower)",
"Ellenberg N (upper)",
"Ellenberg N (lower)",
"Nature Index for mountains (mod.)",
"Artic fox",
"Ptarmigan",
"Connectivity", # Fragmentering
"Area without technical infrastructure",
"Wolverine",
"Golden eagle",
"Willow grouse",
"NDVI (upper)",
"NDVI (lower)",
"Reindeer",
"Snow cover duration",
"Snow depth",
"Small rodents",
"Vegetation heat requirement",
"Winter rain",
"Primary production",
"Functionally important species and structures",
"Biomass composition between trophic levels",
"Functional composition within trophic levels",
"Biological diversity",
"Landscape ecological patterns",
"Abiotic conditions",
"Overall Ecological Condition",
"Overall Ecological Condition (weighted)")Saving
Extra column for name, matched from temp_names
all_Q$names_plot <- temp_names$plot[match(all_Q$ind, temp_names$work)]write.csv(all_Q, "../output/tables/all_Q.csv")
write.csv(tot_ind, "../output/tables/tot_ind.csv")
write.csv(tot_ind_w, "../output/tables/tot_ind_w.csv")
write.csv(agg_ind, "../output/tables/agg_ind.csv")
#write.csv(weights_all, "../output/tables/weights_all.csv")all_Q <- read.csv("../output/tables/all_Q.csv", encoding = "UTF-16")
all_Q <- all_Q[,-1]
tot_ind <- read.csv("../output/tables/tot_ind.csv")
tot_ind <- tot_ind[,-1]
tot_ind_w <- read.csv("../output/tables/tot_ind_w.csv")
tot_ind_w <- tot_ind_w[,-1]
agg_ind <- read.csv("../output/tables/agg_ind.csv")
agg_ind <- agg_ind[,-1]
ind_char <- read.csv("../output/ind_char.csv")
ind_char <- ind_char[,-1]
allind <- read.csv("../output/allind_temp.csv")
allind <- allind[,-1]
weights <- read.csv("../output/weights_temp.csv")
weights <- weights[,-1]Omit the extra name column (only for better understanding of backup csv file)
all_Q <- dplyr::select(all_Q,
-names_plot)Step 9: Plotting
Ignorer per(iodene)
all_Q <- dplyr::select(all_Q, -per)Define plotting order
myOrder <- c(
"NDVI (øvre)",
"NDVI (nedre)",
"Rein",
"Smågnagere",
"Lirype",
"Fjellrype",
"Fjellrev",
"Jerv",
"Kongeørn",
"Fravær av fremmede arter",
"Areal uten tekniske inngrep",
"Konnektivitet",
"Naturindeks for fjell (mod.)",
"Ellenberg N (øvre)",
"Ellenberg N (nedre)",
"Ellenberg L (øvre)",
"Ellenberg L (nedre)",
"Vegetasjonens varmekrav",
"Areal av isbreer",
"Snødybde",
"Snødekkets varighet",
"Vinterregn",
"Total",
"Total (vektet pr. egenskap)",
"Primærproduksjon",
"Fordeling av biomasse i ulike trofiske nivå",
"Funksjonelle grupper innen trofiske nivåer",
"Funksjonelt viktige arter og biofysiske strukturer",
"Landskapsøkologiske mønstre",
"Biologisk mangfold",
"Abiotiske forhold"
)Norsk figure
regions <- unlist(unique(allind$reg))
for (j in 1:length(regions)){
TEMP_Q <- rbind(
all_Q[all_Q$reg==regions[j],],
tot_ind[tot_ind$reg==regions[j],],
tot_ind_w[tot_ind_w$reg==regions[j],],
agg_ind[agg_ind$reg==regions[j],])
TEMP_Q$ind <- temp_names$plot[match(TEMP_Q$ind, temp_names$work)]
TEMP_Q <- TEMP_Q[order(match(TEMP_Q$ind, myOrder)),]
TEMP_Q$ind_2 <- c(1:nrow(TEMP_Q))
png(paste("../output/aggregated_plots/","region_", regions[j],"_", ".png", sep=""),
units="in", width=8, height=14, res=600)
par(xaxt="n", mfrow=c(1,1), par(mar=c(25,5,1,1)))
plot(TEMP_Q$med~TEMP_Q$ind_2, ylab="Skalerte indikatorverdier", xlab="",
type="n", ylim=c(0,1.15), cex.lab=1.2)
abline(h=1, col="blue", lwd=2)
abline(h=0.6, col="red", lwd=2, lty=2)
abline(v=length(unique(allind$X))+0.5, col="black", lwd=2)
abline(v=length(unique(allind$X))+1.5, col="grey", lwd=1, lty=2)
abline(v=length(unique(allind$X))+2.5, col="black", lwd=2)
pex <- 1.8 # median point cex
lar <- 0.0 # length of CI arrows
war <- 5 # width of CI arrows
lcol <- "dark grey"
bg1 <- "white" # ind point col
bg2 <- "black" # char point col
# add values and CI's
arrows(TEMP_Q$ind_2[1:(nind)],TEMP_Q$med[1:(nind)],TEMP_Q$ind_2[1:(nind)],
TEMP_Q$upp[1:(nind)], angle=90, length=lar, lwd= war, col=lcol)
arrows(TEMP_Q$ind_2[1:(nind)],TEMP_Q$med[1:(nind)],TEMP_Q$ind_2[1:(nind)],
TEMP_Q$low[1:(nind)], angle=90, length=lar, lwd= war, col=lcol)
points(TEMP_Q$ind_2[1:(nind)],TEMP_Q$med[1:(nind)],
pch=21, bg=bg1, cex=pex)
# unweighted
arrows(TEMP_Q$ind_2[nind+1],TEMP_Q$med[nind+1],TEMP_Q$ind_2[nind+1],
TEMP_Q$upp[nind+1], angle=90, length=lar, lwd= war, col=lcol)
arrows(TEMP_Q$ind_2[nind+1],TEMP_Q$med[nind+1],TEMP_Q$ind_2[nind+1],
TEMP_Q$low[nind+1], angle=90, length=lar, lwd= war, col=lcol)
points(TEMP_Q$ind_2[nind+1],TEMP_Q$med[nind+1],
pch=23, bg=bg1, cex=pex)
# weighted
arrows(TEMP_Q$ind_2[nind+2],TEMP_Q$med[nind+2],TEMP_Q$ind_2[nind+2],
TEMP_Q$upp[nind+2], angle=90, length=lar, lwd= war, col=lcol)
arrows(TEMP_Q$ind_2[nind+2],TEMP_Q$med[nind+2],TEMP_Q$ind_2[nind+2],
TEMP_Q$low[nind+2], angle=90, length=lar, lwd= war, col=lcol)
points(TEMP_Q$ind_2[nind+2],TEMP_Q$med[nind+2],
pch=23, bg=bg2, cex=pex)
arrows(TEMP_Q$ind_2[nind+3:nrow(TEMP_Q)],TEMP_Q$med[nind+3:nrow(TEMP_Q)],TEMP_Q$ind_2[nind+3:nrow(TEMP_Q)],
TEMP_Q$upp[nind+3:nrow(TEMP_Q)], angle=90, length=lar, lwd= war, col=lcol)
arrows(TEMP_Q$ind_2[nind+3:nrow(TEMP_Q)],TEMP_Q$med[nind+3:nrow(TEMP_Q)],TEMP_Q$ind_2[nind+3:nrow(TEMP_Q)],
TEMP_Q$low[nind+3:nrow(TEMP_Q)], angle=90, length=lar, lwd= war, col=lcol)
points(TEMP_Q$ind_2[nind+3:nrow(TEMP_Q)],TEMP_Q$med[nind+3:nrow(TEMP_Q)],
pch=21, bg=bg2, cex=pex)
# Name axes #
lablist.x<-as.vector(TEMP_Q$ind)
axis(1, at=seq(1, nrow(TEMP_Q), by=1), labels = FALSE)
# Change first command number to align labels
text(TEMP_Q$ind_2+0.40, par("usr")[3]-0.02, labels = lablist.x,pos=2,srt = 90, xpd = TRUE, cex=1.2)
dev.off()
}
#rm(list=ls())
## tempdir()
#dir(tempdir())
#unlink(paste0(normalizePath(tempdir()), "/", dir(tempdir())), recursive = TRUE)
#dir(tempdir())
Engelsk figur
regions <- unlist(unique(allind$reg))
for (j in 1:length(regions)){
TEMP_Q <- rbind(
all_Q[all_Q$reg==regions[j],],
tot_ind[tot_ind$reg==regions[j],],
tot_ind_w[tot_ind_w$reg==regions[j],],
agg_ind[agg_ind$reg==regions[j],])
TEMP_Q$temp <- temp_names$plot[match(TEMP_Q$ind, temp_names$work)]
TEMP_Q$ind <- temp_names2$plot[match(TEMP_Q$ind, temp_names2$work)]
TEMP_Q <- TEMP_Q[order(match(TEMP_Q$temp, myOrder)),]
TEMP_Q$ind_2 <- c(1:nrow(TEMP_Q))
png(paste("../output/aggregated_plots/englishPlots/fjell/","region_", regions[j],"_", ".png", sep=""),
units="in", width=8, height=14, res=600)
par(xaxt="n", mfrow=c(1,1), par(mar=c(25,5,1,1)))
plot(TEMP_Q$med~TEMP_Q$ind_2, ylab="Scaled indicator values", xlab="",
type="n", ylim=c(0,1.15), cex.lab=1.2)
abline(h=1, col="blue", lwd=2)
abline(h=0.6, col="red", lwd=2, lty=2)
abline(v=length(unique(allind$X))+0.5, col="black", lwd=2)
abline(v=length(unique(allind$X))+1.5, col="grey", lwd=1, lty=2)
abline(v=length(unique(allind$X))+2.5, col="black", lwd=2)
pex <- 1.8 # median point cex
lar <- 0.0 # length of CI arrows
war <- 5 # width of CI arrows
lcol <- "dark grey"
bg1 <- "white" # ind point col
bg2 <- "black" # char point col
# add values and CI's
arrows(TEMP_Q$ind_2[1:(nind)],TEMP_Q$med[1:(nind)],TEMP_Q$ind_2[1:(nind)],
TEMP_Q$upp[1:(nind)], angle=90, length=lar, lwd= war, col=lcol)
arrows(TEMP_Q$ind_2[1:(nind)],TEMP_Q$med[1:(nind)],TEMP_Q$ind_2[1:(nind)],
TEMP_Q$low[1:(nind)], angle=90, length=lar, lwd= war, col=lcol)
points(TEMP_Q$ind_2[1:(nind)],TEMP_Q$med[1:(nind)],
pch=21, bg=bg1, cex=pex)
# unweighted
arrows(TEMP_Q$ind_2[nind+1],TEMP_Q$med[nind+1],TEMP_Q$ind_2[nind+1],
TEMP_Q$upp[nind+1], angle=90, length=lar, lwd= war, col=lcol)
arrows(TEMP_Q$ind_2[nind+1],TEMP_Q$med[nind+1],TEMP_Q$ind_2[nind+1],
TEMP_Q$low[nind+1], angle=90, length=lar, lwd= war, col=lcol)
points(TEMP_Q$ind_2[nind+1],TEMP_Q$med[nind+1],
pch=23, bg=bg1, cex=pex)
# weighted
arrows(TEMP_Q$ind_2[nind+2],TEMP_Q$med[nind+2],TEMP_Q$ind_2[nind+2],
TEMP_Q$upp[nind+2], angle=90, length=lar, lwd= war, col=lcol)
arrows(TEMP_Q$ind_2[nind+2],TEMP_Q$med[nind+2],TEMP_Q$ind_2[nind+2],
TEMP_Q$low[nind+2], angle=90, length=lar, lwd= war, col=lcol)
points(TEMP_Q$ind_2[nind+2],TEMP_Q$med[nind+2],
pch=23, bg=bg2, cex=pex)
arrows(TEMP_Q$ind_2[nind+3:nrow(TEMP_Q)],TEMP_Q$med[nind+3:nrow(TEMP_Q)],TEMP_Q$ind_2[nind+3:nrow(TEMP_Q)],
TEMP_Q$upp[nind+3:nrow(TEMP_Q)], angle=90, length=lar, lwd= war, col=lcol)
arrows(TEMP_Q$ind_2[nind+3:nrow(TEMP_Q)],TEMP_Q$med[nind+3:nrow(TEMP_Q)],TEMP_Q$ind_2[nind+3:nrow(TEMP_Q)],
TEMP_Q$low[nind+3:nrow(TEMP_Q)], angle=90, length=lar, lwd= war, col=lcol)
points(TEMP_Q$ind_2[nind+3:nrow(TEMP_Q)],TEMP_Q$med[nind+3:nrow(TEMP_Q)],
pch=21, bg=bg2, cex=pex)
# Name axes #
lablist.x<-as.vector(TEMP_Q$ind)
axis(1, at=seq(1, nrow(TEMP_Q), by=1), labels = FALSE)
# Change first command number to align labels
text(TEMP_Q$ind_2+0.40, par("usr")[3]-0.02, labels = lablist.x,pos=2,srt = 90, xpd = TRUE, cex=1.2)
dev.off()
}
#rm(list=ls())
## tempdir()
#dir(tempdir())
#unlink(paste0(normalizePath(tempdir()), "/", dir(tempdir())), recursive = TRUE)
#dir(tempdir())For plotting av indikatorverdier gruppert etter påvirkningsfaktor, se her