6 Functional Plant Indicators - Wetlands

Norwegian name: Vegetasjon og fuktighet/lys/pH/nitrogen

Author and date:

Joachim Töpper

May 2023

Ecosystem	Økologisk egenskap	ECT class
våtmark	Primærproduksjon	Functional state characteristic
våtmark	Abiotiske forhold	Functional state characteristic

Indicators described in this chapter:

Moisture
Light
pH
Nitrogen

6.1 Introduction

Functional plant indicators can be used to describe the functional signature of plant communities by calculating community-weighted means of functional trait values for plant communities (Diekmann 2003). The functional signature of plant communities may be indicative of ecosystem identity, depending on which functional plant indicator we look at (cf. Töpper et al. 2018). For instance, using an indicator for moisture one would find a functional signature with higher moisture values for plant communities in mires compared to e.g. grasslands or forests. Deviations in the functional signature of such an indicator beyond a certain range of indicator values (as there of course is natural variation of functional signatures within an ecosystem type) may be related to a reduction in ecological condition.

Here, we combine functional indicator data with field sampled plant community data from the Norwegian nature monitoring program ANO (Tingstad et al. 2019) for wetland ecosystems. We calculate the functional signature of plant communities in the monitored sites with respect to light, moisture, pH and nitrogen. These functional signatures are then compared to reference distributions of functional signature, separately for each wetland ecosystem type, calculated from ‘generalized species lists’ underlying the Norwegian categorization system for eco-diversity (Halvorsen et al. 2020). These plant functional type indicators are developed following the principles and technical protocol of the IBECA framework (Jakobsson et al. 2021, Töpper & Jakobsson 2021).

Parallel indicator sets are explored and developed for semi-natural ecosystems and naturally open ecosystems.

6.2 About the underlying data

We use three sets of data for building indicators for ecological condition:

as test data we use plant community data from the ANO monitoring scheme (cf. Tingstad et al. 2019)
as reference data we use generalized species lists from the Norwegian categorization system for eco-diversity (NiN) (cf. Halvorsen et al. 2020)
Swedish plant indicator data for light, moisture, pH, and nitrogen from Tyler et al. (2021)

ANO monitoring data: ANO stands for areal-representativ naturovervåking, i.e. area representative nature monitoring. 1000 sites are randomly distributed across mainland Norway and visitied within a 5-year cycle. Each ANO site spans a 500 x 500 m grid cell, and the data collection at each ANO site includes up to 12 evenly distributed vegetation analyses in 1 x 1 m plots (up to 12, because some of these evenly distributed points may be in water or otherwise inaccessible). For the vegetation analyses, the cover of each vascular plant species in the plot is recorded. Every vegetation analysis is accompanied by an assessment of the ecosystem the plot lies in, including ecosystem type and some additional variables related to ecosystem-specific drivers of state. In the wetland-analysis in this document, we only use the plots which were classified as one of the wetland ecosystem types in the Norwegian categorization system for eco-diversity (NiN). In the analysis in this document, we use the data available on Miljødirektoratets kartkatalog (https://kartkatalog.miljodirektoratet.no/Dataset/Details/2054), which comprises data from the first three ANO-years, 2019-2021, and a total of 8887 plots in 498 sites.
NiN reference data: The generalized species lists underlying the ecosystem categorization in NiN represent expert-compiled species lists based on empirical evidence from the literature and expert knowledge of the systems and their species. In these lists, every species is assigned an abundance value on a 6-step scale, with each step representing a range for the ‘expected combination of frequency and cover’ of occurrence
- 1 = < 1/32
- 2 = 1/32 - 1/8
- 3 = 1/8 - 3/8
- 4 = 3/8 - 4/5
- 5 = 3/8 - 4/5 + dominance
- 6 = > 4/5

For the purpose of this project, these steps are simplified to maximum expected combination of frequency and cover, whereby steps 4 & 5 are assigned 0.6 and 0.8, respectively, in order to distinguish between them.

The Swedish plant indicator set published by Tyler et al. (2021) contains a large collection of plant indicators based on the Swedish flora, which is well representative of the Norwegian flora as well. From this set, we decided to include indicator data for light, moisture, pH, and nitrogen for wetlands as these are thought to be subject to potential change due to shrub/tree encroachment, drainage, and pollution.

6.2.1 Representativity in time and space

For wetlands, the ANO data contain 1351 plots in 330 sites, in principle distributed randomly across the country. As wetlands occur more often in certain regions of Norway than in others, the amount of plots and sites is not equal among Norway’s five regions. The 1351 plots are distributed across regions in the following way:

Northern Norway: 416
Central Norway: 361
Eastern Norway: 294
Western Norway: 150
Southern Norway: 126

6.2.2 Temporal coverage

The evaluation data cover the first three years, 2019-2021, of the first 5-year-cycle in the ANO monitoring scheme. Thus there is no actual time series to these data, and the indicator evaluation does therefore not include any temporal analyses.

6.3 Collinearities with other indicators

Collinearity could occur with other indicators that are responding to the same drivers/pressures, but no really problematic issues have been discussed.

6.4 Reference state and values

6.4.1 Reference state

The reference state is defined via the functional signature of the generalized species lists in NiN (see also Töpper et al. 2018). By bootstrapping the species lists (see details further below) and calculating community-weighted means of functional plant indicators for every re-sampled community, we describe the reference state as a distribution of indicator values for each respective plant functional indicator. These distributions are calculated for minor ecosystem types ( grunntyper or kartleggingsenheter at a 1:5000 mapping scale) within the major ecosystem types (hovedtyper) in NiN. A more extensive discussion on the use of reference communities can be found in Jakobsson et al. (2020).

6.4.2 Reference values, thresholds for defining good ecological condition, minimum and/or maximum values

In this analysis, we derive scaling values from statistical (here, non-parametric) distributions (see Jakobsson et al. 2010). For each ecosystem sub-type and plant functional indicator, the reference value is defined as the median value of the indicator value distribution. As in most cases the distributions naturally are two-sided (but see the Heat-requirement indicator in the mountain assessment for an example of a one-sided plant functional indicator, Fremstad et al. 2022), and deviation from the optimal state thus may occur in both direction (e.g. indicating too low or too high pH), we need to define two threshold values for good ecological condition as well as both a minimum and maximum value. In line with previous assessments of ecological condition for Norwegian forests and mountains, we define a lower and an upper threshold value via the 95% confidence interval of the reference distribution, i.e. its 0.025 and 0.975 quantiles. The minimum and maximum values are given by the minimum and maximum of the possible indicator values for each respective plant functional indicator. For details on the scaling principles in IBECA, please see Töpper & Jakobsson (2021).

6.5 Uncertainties

We can calculate a mean indicator value (after scaling) for every region (or any other delimited area of interest) as well as its corresponding standard error as a measure of spatial uncertainty for a geographical area.

6.6 References

Diekmann, M. 2003. Species indicator values as an important tool in applied plant ecology - a review. Basic and Applied Ecology 4: 493-506, doi:10.1078/1439-1791-00185

Framstad, E., Kolstad, A. L., Nybø, S., Töpper, J. & Vandvik, V. 2022. The condition of forest and mountain ecosystems in Norway. Assessment by the IBECA method. NINA Report 2100. Norwegian Institute for Nature Research.

Halvorsen, R., Skarpaas, O., Bryn, A., Bratli, H., Erikstad, L., Simensen, T., & Lieungh, E. (2020). Towards a systematics of ecodiversity: The EcoSyst framework. Global Ecology and Biogeography, 29(11), 1887-1906. doi:10.1111/geb.13164

Jakobsson, S., Töpper, J.P., Evju, M., Framstad, E., Lyngstad, A., Pedersen, B., Sickel, H., Sverdrup-Thygeson, A., Vandvik. V., Velle, L.G., Aarrestad, P.A. & Nybø, S. 2020. Setting reference levels and limits for good ecological condition in terrestrial ecosystems. Insights from a case study based on the IBECA approach. Ecological Indicators 116: 106492.

Jakobsson, S., Evju, M., Framstad, E., Imbert, A., Lyngstad, A., Sickel, H., Sverdrup-Thygeson, A., Töpper, J., Vandvik, V., Velle, L.G., Aarrestad, P.A. & Nybø, S. 2021. An index-based assessment of ecological condition and its links to international frameworks. Ecological Indicators 124: 107252.

Tingstad, L., Evju, M., Sickel, H., & Töpper, J. 2019. Utvikling av nasjonal arealrepresentativ naturovervåking (ANO). Forslag til gjennomføring, protokoller og kostnadsvurderinger med utgangspunkt i erfaringer fra uttesting i Trøndelag. NINA Rapport 1642.

Töpper, J. & Jakobsson, S. 2021. The Index-Based Ecological Condition Assessment (IBECA) - Technical protocol, version 1.0. NINA Report 1967. Norwegian Institute for Nature Research.

Töpper, J., Velle, L.G. & Vandvik, V. 2018. Developing a method for assessment of ecological state based on indicator values after Ellenberg and Grime (revised edition). NINA Report 1529b. Norwegian Institute for Nature Research.

Tyler, T., Herbertsson, L., Olofsson, J., & Olsson, P. A. 2021. Ecological indicator and traits values for Swedish vascular plants. Ecological In-dicators, 120. doi:10.1016/j.ecolind.2020.106923

6.7 Analyses

6.7.1 Data sets

Import ANO data:

# Add data from cache
ANO.sp<-readRDS(paste0(here::here(),"/data/cache/ANO.sp.RDS"))
ANO.geo<-readRDS(paste0(here::here(),"/data/cache/ANO.geo.RDS"))

# st_layers(dsn = "/data/P-Prosjekter2/41201785_okologisk_tilstand_2022_2023/data/Naturovervaking_eksport.gdb")
# ANO.sp <- st_read("/data/P-Prosjekter2/41201785_okologisk_tilstand_2022_2023/data/Naturovervaking_eksport.gdb",
#                    layer="ANO_Art")
# ANO.geo <- st_read("/data/P-Prosjekter2/41201785_okologisk_tilstand_2022_2023/data/Naturovervaking_eksport.gdb",
#                   layer="ANO_SurveyPoint")
# head(ANO.sp)
# head(ANO.geo)

Import Plant indicators from Tyler et al. (2021):

## Tyler indicator data
## data from cache
ind.dat <-readRDS(paste0(here::here(),"/data/cache/ind.dat.RDS"))
# ind.dat <- read.table("/data/P-Prosjekter2/41201785_okologisk_tilstand_2022_2023/data/functional plant indicators/Tyler et al_Swedish plant indicators.txt",
#                         sep = '\t', header=T, quote = '')
# head(ind.dat)

Import Generalized species lists NiN

# data from cache
Eco_State<-readRDS(paste0(here::here(),"/data/cache/Eco_State.RDS"))
# load("/data/P-Prosjekter2/41201785_okologisk_tilstand_2022_2023/functional plant indicators/reference from NiN/Eco_State.RData")
# str(Eco_State)

6.7.1.1 Data handling

Checking for errors
Checking species nomenclature in the different species lists to make species and indicator data possible to merge
Merging indicator data with monitoring data and indicator data with reference data (not shown here, but documented in the code)


### Plant indicator data
names(ind.dat)
names(ind.dat)[1] <- 'species'
ind.dat$species <- as.factor(ind.dat$species)
summary(ind.dat$species)
ind.dat <- ind.dat[!is.na(ind.dat$species),] # removing species-NAs
ind.dat[,'species.orig'] <- ind.dat[,'species'] # make a backup of the original species variable
ind.dat[,'species'] <- word(ind.dat[,'species'], 1,2) # trimming away sub-species & co, and descriptor info

# the trimming above leaves some duplicates that need to be taken care of
ind.dat[duplicated(ind.dat[,'species']),"species"]
ind.dat.dup <- ind.dat[duplicated(ind.dat[,'species']),"species"]
ind.dat[ind.dat$species %in% ind.dat.dup,c("Light","Moisture","Soil_reaction_pH","Nitrogen","species.orig","species")]
ind.dat <- ind.dat %>% filter( !(species.orig %in% list("Ammophila arenaria x Calamagrostis epigejos",
                                                            "Anemone nemorosa x ranunculoides",
                                                            "Armeria maritima ssp. elongata",
                                                            "Asplenium trichomanes ssp. quadrivalens",
                                                            "Calystegia sepium ssp. spectabilis",
                                                            "Campanula glomerata 'Superba'",
                                                            "Dactylorhiza maculata ssp. fuchsii",
                                                            "Erigeron acris ssp. droebachensis",
                                                            "Erigeron acris ssp. politus",
                                                            "Erysimum cheiranthoides L. ssp. alatum",
                                                            "Euphrasia nemorosa x stricta var. brevipila",
                                                            "Galium mollugo x verum",
                                                            "Geum rivale x urbanum",
                                                            "Hylotelephium telephium (ssp. maximum)",
                                                            "Juncus alpinoarticulatus ssp. rariflorus",
                                                            "Lamiastrum galeobdolon ssp. argentatum",
                                                            "Lathyrus latifolius ssp. heterophyllus",
                                                            "Medicago sativa ssp. falcata",
                                                            "Medicago sativa ssp. x varia",
                                                            "Monotropa hypopitys ssp. hypophegea",
                                                            "Ononis spinosa ssp. hircina",
                                                            "Ononis spinosa ssp. procurrens",
                                                            "Pilosella aurantiaca ssp. decolorans",
                                                            "Pilosella aurantiaca ssp. dimorpha",
                                                            "Pilosella cymosa ssp. gotlandica",
                                                            "Pilosella cymosa ssp. praealta",
                                                            "Pilosella officinarum ssp. peleteranum",
                                                            "Poa x jemtlandica (Almq.) K. Richt.",
                                                            "Poa x herjedalica Harry Sm.",
                                                            "Ranunculus peltatus ssp. baudotii",
                                                            "Sagittaria natans x sagittifolia",
                                                            "Salix repens ssp. rosmarinifolia",
                                                            "Stellaria nemorum L. ssp. montana",
                                                            "Trichophorum cespitosum ssp. germanicum")
) )

# checking duplicates again
ind.dat[duplicated(ind.dat[,'species']),"species"]
ind.dat.dup <- ind.dat[duplicated(ind.dat[,'species']),"species"]
ind.dat[ind.dat$species %in% ind.dat.dup,c("Light","Moisture","Soil_reaction_pH","Nitrogen","species.orig","species")]
# getting rid of sect. for Hieracium
ind.dat <- ind.dat %>% mutate(species=gsub("sect. ","",species.orig))
ind.dat[,'species'] <- word(ind.dat[,'species'], 1,2)

ind.dat[duplicated(ind.dat[,'species']),"species"]
ind.dat.dup <- ind.dat[duplicated(ind.dat[,'species']),"species"]
ind.dat[ind.dat$species %in% ind.dat.dup,c("Light","Moisture","Soil_reaction_pH","Nitrogen","species.orig","species")]
# only hybrids left -> get rid of these
ind.dat <- ind.dat[!duplicated(ind.dat[,'species']),]
ind.dat[duplicated(ind.dat[,'species']),"species"]

ind.dat$species <- as.factor(ind.dat$species)
summary(ind.dat$species)
# no duplicates left

head(ind.dat)



### ANO monitoring data

head(ANO.sp)
head(ANO.geo)

## fix NiN information
ANO.geo$hovedtype_rute <- substr(ANO.geo$kartleggingsenhet_1m2,1,3) # take the 3 first characters
ANO.geo$hovedtype_rute <- gsub("-", "", ANO.geo$hovedtype_rute) # remove hyphen
unique(as.factor(ANO.geo$hovedtype_rute))

## fix NiN-variables
colnames(ANO.geo)
colnames(ANO.geo)[42:47] <- c("groeftingsintensitet",
                              "bruksintensitet",
                              "beitetrykk",
                              "slatteintensitet",
                              "tungekjoretoy",
                              "slitasje")
head(ANO.geo)

# remove variable code in the data
ANO.geo$groeftingsintensitet <- gsub("7GR-GI_", "", ANO.geo$groeftingsintensitet) 
unique(ANO.geo$groeftingsintensitet)
ANO.geo$groeftingsintensitet <- gsub("X", "NA", ANO.geo$groeftingsintensitet)
unique(ANO.geo$groeftingsintensitet)
ANO.geo$groeftingsintensitet <- as.numeric(ANO.geo$groeftingsintensitet)
unique(ANO.geo$groeftingsintensitet)

ANO.geo$bruksintensitet <- gsub("7JB-BA_", "", ANO.geo$bruksintensitet) 
unique(ANO.geo$bruksintensitet)
ANO.geo$bruksintensitet <- gsub("X", "NA", ANO.geo$bruksintensitet)
unique(ANO.geo$bruksintensitet)
ANO.geo$bruksintensitet <- as.numeric(ANO.geo$bruksintensitet)
unique(ANO.geo$bruksintensitet)

ANO.geo$beitetrykk <- gsub("7JB-BT_", "", ANO.geo$beitetrykk) 
unique(ANO.geo$beitetrykk)
ANO.geo$beitetrykk <- gsub("X", "NA", ANO.geo$beitetrykk)
unique(ANO.geo$beitetrykk)
ANO.geo$beitetrykk <- as.numeric(ANO.geo$beitetrykk)
unique(ANO.geo$beitetrykk)

ANO.geo$slatteintensitet <- gsub("7JB-SI_", "", ANO.geo$slatteintensitet) 
unique(ANO.geo$slatteintensitet)
ANO.geo$slatteintensitet <- gsub("X", "NA", ANO.geo$slatteintensitet)
unique(ANO.geo$slatteintensitet)
ANO.geo$slatteintensitet <- as.numeric(ANO.geo$slatteintensitet)
unique(ANO.geo$slatteintensitet)

ANO.geo$tungekjoretoy <- gsub("7TK_", "", ANO.geo$tungekjoretoy) 
unique(ANO.geo$tungekjoretoy)
ANO.geo$tungekjoretoy <- gsub("X", "NA", ANO.geo$tungekjoretoy)
unique(ANO.geo$tungekjoretoy)
ANO.geo$tungekjoretoy <- as.numeric(ANO.geo$tungekjoretoy)
unique(ANO.geo$tungekjoretoy)

ANO.geo$slitasje <- gsub("7SE_", "", ANO.geo$slitasje) 
unique(ANO.geo$slitasje)
ANO.geo$slitasje <- gsub("X", "NA", ANO.geo$slitasje)
unique(ANO.geo$slitasje)
ANO.geo$slitasje <- as.numeric(ANO.geo$slitasje)
unique(ANO.geo$slitasje)

## check that every point is present only once
length(levels(as.factor(ANO.geo$ano_flate_id)))
length(levels(as.factor(ANO.geo$ano_punkt_id)))
summary(as.factor(ANO.geo$ano_punkt_id))
# there's many double presences, probably some wrong registrations of point numbers,
# but also double registrations (e.g. ANO0159_55)
# CHECK THIS when preparing ecosystem-datasets for scaling




# fix species names
ANO.sp$Species <- ANO.sp$art_navn
unique(as.factor(ANO.sp$Species))
ANO.sp[,'Species'] <- word(ANO.sp[,'Species'], 1,2) # lose subspecies
ANO.sp$Species <- str_to_title(ANO.sp$Species) # make first letter capital
ANO.sp$Species <- gsub("( .*)","\\L\\1",ANO.sp$Species,perl=TRUE) # make capital letters after hyphen to lowercase
ANO.sp$Species <- gsub("( .*)","\\L\\1",ANO.sp$Species,perl=TRUE) # make capital letters after space to lowercase

## merge species data with indicators
ANO.sp.ind <- merge(x=ANO.sp[,c("Species", "art_dekning", "ParentGlobalID")], 
                y= ind.dat[,c("species", "Light", "Moisture", "Soil_reaction_pH", "Nitrogen")],
                by.x="Species", by.y="species", all.x=T)
summary(ANO.sp.ind)


## checking which species didn't find a match
unique(ANO.sp.ind[is.na(ANO.sp.ind$Moisture),'Species'])

# fix species name issues
ind.dat <- ind.dat %>% 
  mutate(species=str_replace(species,"Aconitum lycoctonum", "Aconitum septentrionale")) %>% 
  mutate(species=str_replace(species,"Carex simpliciuscula", "Kobresia simpliciuscula")) %>%
  mutate(species=str_replace(species,"Carex myosuroides", "Kobresia myosuroides")) %>%
  mutate(species=str_replace(species,"Clinopodium acinos", "Acinos arvensis")) %>%
  mutate(species=str_replace(species,"Artemisia rupestris", "Artemisia norvegica")) %>%
  mutate(species=str_replace(species,"Cherleria biflora", "Minuartia biflora"))

ANO.sp <- ANO.sp %>% 
  mutate(Species=str_replace(Species,"Arctous alpinus", "Arctous alpina")) %>%
  mutate(Species=str_replace(Species,"Betula tortuosa", "Betula pubescens")) %>%
  mutate(Species=str_replace(Species,"Blysmopsis rufa", "Blysmus rufus")) %>%
  mutate(Species=str_replace(Species,"Cardamine nymanii", "Cardamine pratensis")) %>%
  mutate(Species=str_replace(Species,"Carex adelostoma", "Carex buxbaumii")) %>%
  mutate(Species=str_replace(Species,"Carex leersii", "Carex echinata")) %>%
  mutate(Species=str_replace(Species,"Carex paupercula", "Carex magellanica")) %>%
  mutate(Species=str_replace(Species,"Carex simpliciuscula", "Kobresia simpliciuscula")) %>%
  mutate(Species=str_replace(Species,"Carex viridula", "Carex flava")) %>%
  mutate(Species=str_replace(Species,"Chamaepericlymenum suecicum", "Cornus suecia")) %>%
  mutate(Species=str_replace(Species,"Cicerbita alpina", "Lactuca alpina")) %>%
  mutate(Species=str_replace(Species,"Empetrum hermaphroditum", "Empetrum nigrum")) %>%
  mutate(Species=str_replace(Species,"Festuca prolifera", "Festuca rubra")) %>%
  mutate(Species=str_replace(Species,"Galium album", "Galium mollugo")) %>%
  mutate(Species=str_replace(Species,"Galium elongatum", "Galium palustre")) %>%
  mutate(Species=str_replace(Species,"Helictotrichon pratense", "Avenula pratensis")) %>%
  mutate(Species=str_replace(Species,"Helictotrichon pubescens", "Avenula pubescens")) %>%
  mutate(Species=str_replace(Species,"Hieracium alpina", "Hieracium Alpina")) %>%
  mutate(Species=str_replace(Species,"Hieracium alpinum", "Hieracium Alpina")) %>%
  mutate(Species=str_replace(Species,"Hieracium hieracium", "Hieracium Hieracium")) %>%
  mutate(Species=str_replace(Species,"Hieracium hieracioides", "Hieracium umbellatum")) %>%
  mutate(Species=str_replace(Species,"Hieracium murorum", "Hieracium Vulgata")) %>%
  mutate(Species=str_replace(Species,"Hieracium oreadea", "Hieracium Oreadea")) %>%
  mutate(Species=str_replace(Species,"Hieracium prenanthoidea", "Hieracium Prenanthoidea")) %>%
  mutate(Species=str_replace(Species,"Hieracium vulgata", "Hieracium Vulgata")) %>%
  mutate(Species=str_replace(Species,"Hieracium pilosella", "Pilosella officinarum")) %>%
  mutate(Species=str_replace(Species,"Hieracium vulgatum", "Hieracium umbellatum")) %>%
  mutate(Species=str_replace(Species,"Hierochloã« alpina", "Hierochloë alpina")) %>%
  mutate(Species=str_replace(Species,"Hierochloã« hirta", "Hierochloë hirta")) %>%
  mutate(Species=str_replace(Species,"Hierochloã« odorata", "Hierochloë odorata")) %>%
  mutate(Species=str_replace(Species,"Listera cordata", "Neottia cordata")) %>%
  mutate(Species=str_replace(Species,"Leontodon autumnalis", "Scorzoneroides autumnalis")) %>%
  mutate(Species=str_replace(Species,"Loiseleuria procumbens", "Kalmia procumbens")) %>%
  mutate(Species=str_replace(Species,"Mycelis muralis", "Lactuca muralis")) %>%
  mutate(Species=str_replace(Species,"Omalotheca supina", "Gnaphalium supinum")) %>%
  mutate(Species=str_replace(Species,"Omalotheca norvegica", "Gnaphalium norvegicum")) %>%
  mutate(Species=str_replace(Species,"Omalotheca sylvatica", "Gnaphalium sylvaticum")) %>%
  mutate(Species=str_replace(Species,"Oreopteris limbosperma", "Thelypteris limbosperma")) %>%
  mutate(Species=str_replace(Species,"Oxycoccus microcarpus", "Vaccinium microcarpum")) %>%
  mutate(Species=str_replace(Species,"Oxycoccus palustris", "Vaccinium oxycoccos")) %>%
  mutate(Species=str_replace(Species,"Phalaris minor", "Phalaris arundinacea")) %>%
  mutate(Species=str_replace(Species,"Pinus unicinata", "Pinus mugo")) %>%
  mutate(Species=str_replace(Species,"Poa alpigena", "Poa pratensis")) %>%
  mutate(Species=str_replace(Species,"Poa angustifolia", "Poa pratensis")) %>%
  mutate(Species=str_replace(Species,"Pyrola grandiflora", "Pyrola rotundifolia")) %>%
  mutate(Species=str_replace(Species,"Rumex alpestris", "Rumex acetosa")) %>%
  mutate(Species=str_replace(Species,"Syringa emodi", "Syringa vulgaris")) %>%
  mutate(Species=str_replace(Species,"Taraxacum crocea", "Taraxacum officinale")) %>%
  mutate(Species=str_replace(Species,"Taraxacum croceum", "Taraxacum officinale")) %>%
  mutate(Species=str_replace(Species,"Trientalis europaea", "Lysimachia europaea")) %>%
  mutate(Species=str_replace(Species,"Trifolium pallidum", "Trifolium pratense"))

## merge species data with indicators
ANO.sp.ind <- merge(x=ANO.sp[,c("Species", "art_dekning", "ParentGlobalID")], 
                    y= ind.dat[,c("species", "Light", "Moisture", "Soil_reaction_pH", "Nitrogen")],
                    by.x="Species", by.y="species", all.x=T)
summary(ANO.sp.ind)
# checking which species didn't find a match
unique(ANO.sp.ind[is.na(ANO.sp.ind$Moisture),'Species'])
# don't find synonyms for these in the ind lists

## trimming away the points without information on NiN, species or cover
ANO.sp.ind <- ANO.sp.ind[!is.na(ANO.sp.ind$Species),]
ANO.sp.ind <- ANO.sp.ind[!is.na(ANO.sp.ind$art_dekning),]


summary(ANO.sp.ind)
head(ANO.sp.ind)
rm(ANO.sp)


### NiN reference data - data handling

## generalized species lists from NiN
str(Eco_State)

# species
Eco_State$Concept_Data$Species$Species_List$species
# environments
t(Eco_State$Concept_Data$Env$Env_Data)
# abundances
t(Eco_State$Concept_Data$Species$Species_Data)

## transposing abundance data
NiN.sp <- t(Eco_State$Concept_Data$Species$Species_Data)
NiN.sp <- as.data.frame(NiN.sp)
NiN.sp$sp <- as.factor(as.vector(Eco_State$Concept_Data$Species$Species_List$species))
# only genus and species name
NiN.sp$sp <- word(NiN.sp$sp, 1,2)
NiN.sp$spgr <- as.factor(as.vector(Eco_State$Concept_Data$Species$Species_List$art.code))
# if relevant, trimming to desired species groups (for forests e.g. removing trees)
#NiN.sp <- NiN.sp[NiN.sp$spgr!="a1a",]

## environment data
NiN.env <- Eco_State$Concept_Data$Env$Env_Data

## merging with indicator values
NiN.sp.ind <- merge(NiN.sp,ind.dat, by.x="sp", by.y="species", all.x=T)
summary(NiN.sp.ind)

NiN.sp.ind[NiN.sp.ind==999] <- NA

## checking which species didn't find a match
unique(NiN.sp.ind[is.na(NiN.sp.ind$Moisture) & NiN.sp.ind$spgr %in% list("a1a","a1b","a1c"),'sp'])

## fix species name issues
ind.dat <- ind.dat %>% 
  mutate(species=str_replace(species,"Aconitum lycoctonum", "Aconitum septentrionale")) %>% 
  mutate(species=str_replace(species,"Carex simpliciuscula", "Kobresia simpliciuscula")) %>%
  mutate(species=str_replace(species,"Carex myosuroides", "Kobresia myosuroides")) %>%
  mutate(species=str_replace(species,"Clinopodium acinos", "Acinos arvensis")) %>%
  mutate(species=str_replace(species,"Artemisia rupestris", "Artemisia norvegica")) %>%
  mutate(species=str_replace(species,"Cherleria biflora", "Minuartia biflora"))

NiN.sp <- NiN.sp %>% 
  mutate(sp=str_replace(sp,"Aconitum lycoctonum", "Aconitum septentrionale")) %>% 
  mutate(sp=str_replace(sp,"Anagallis minima", "Lysimachia minima")) %>% 
  mutate(sp=str_replace(sp,"Arctous alpinus", "Arctous alpina")) %>%
  mutate(sp=str_replace(sp,"Betula tortuosa", "Betula pubescens")) %>%
  mutate(sp=str_replace(sp,"Blysmopsis rufa", "Blysmus rufus")) %>%
  mutate(sp=str_replace(sp,"Cardamine nymanii", "Cardamine pratensis")) %>%
  mutate(sp=str_replace(sp,"Carex adelostoma", "Carex buxbaumii")) %>%
  mutate(sp=str_replace(sp,"Carex leersii", "Carex echinata")) %>%
  mutate(sp=str_replace(sp,"Carex paupercula", "Carex magellanica")) %>%
  mutate(sp=str_replace(sp,"Carex simpliciuscula", "Kobresia simpliciuscula")) %>%
  mutate(sp=str_replace(sp,"Carex _vacillans", "Carex vacillans")) %>%
  mutate(sp=str_replace(sp,"Carex viridula", "Carex flava")) %>%
  mutate(sp=str_replace(sp,"Chamaepericlymenum suecicum", "Cornus suecia")) %>%
  mutate(sp=str_replace(sp,"Cornus suecia", "Cornus suecica")) %>%
  mutate(sp=str_replace(sp,"Cicerbita alpina", "Lactuca alpina")) %>%
  mutate(sp=str_replace(sp,"Dactylorhiza sphagnicola", "Dactylorhiza majalis")) %>%
  mutate(sp=str_replace(sp,"Diphasiastrum complanatum", "Lycopodium complanatum")) %>%
  mutate(sp=str_replace(sp,"Elymus alaskanus", "Elymus kronokensis")) %>%
  mutate(sp=str_replace(sp,"Empetrum hermaphroditum", "Empetrum nigrum")) %>%
  mutate(sp=str_replace(sp,"Erigeron eriocephalus", "Erigeron uniflorus")) %>%
  mutate(sp=str_replace(sp,"Festuca prolifera", "Festuca rubra")) %>%
  mutate(sp=str_replace(sp,"Galium album", "Galium mollugo")) %>%
  mutate(sp=str_replace(sp,"Galium elongatum", "Galium palustre")) %>%
  mutate(sp=str_replace(sp,"Glaux maritima", "Lysimachia maritima")) %>%
  mutate(sp=str_replace(sp,"Helictotrichon pratense", "Avenula pratensis")) %>%
  mutate(sp=str_replace(sp,"Helictotrichon pubescens", "Avenula pubescens")) %>%
  mutate(sp=str_replace(sp,"Hieracium alpina", "Hieracium Alpina")) %>%
  mutate(sp=str_replace(sp,"Hieracium alpinum", "Hieracium Alpina")) %>%
  mutate(sp=str_replace(sp,"Hieracium aurantiacum", "Pilosella aurantiaca")) %>%
  mutate(sp=str_replace(sp,"Hieracium hieracium", "Hieracium Hieracium")) %>%
  mutate(sp=str_replace(sp,"Hieracium hieracioides", "Hieracium umbellatum")) %>%
  mutate(sp=str_replace(sp,"Hieracium lactucella", "Pilosella lactucella")) %>%
  mutate(sp=str_replace(sp,"Hieracium murorum", "Hieracium Vulgata")) %>%
  mutate(sp=str_replace(sp,"Hieracium oreadea", "Hieracium Oreadea")) %>%
  mutate(sp=str_replace(sp,"Hieracium prenanthoidea", "Hieracium Prenanthoidea")) %>%
  mutate(sp=str_replace(sp,"Hieracium vulgata", "Hieracium Vulgata")) %>%
  mutate(sp=str_replace(sp,"Hieracium pilosella", "Pilosella officinarum")) %>%
  mutate(sp=str_replace(sp,"Hieracium vulgatum", "Hieracium umbellatum")) %>%
  mutate(sp=str_replace(sp,"Hierochloã« alpina", "Hierochloë alpina")) %>%
  mutate(sp=str_replace(sp,"Hierochloã« hirta", "Hierochloë hirta")) %>%
  mutate(sp=str_replace(sp,"Hierochloã« odorata", "Hierochloë odorata")) %>%
  mutate(sp=str_replace(sp,"Huperzia appressa", "Huperzia selago")) %>%
  mutate(sp=str_replace(sp,"Hylotelephium maximum", "Hylotelephium telephium")) %>%
  mutate(sp=str_replace(sp,"Lappula myosotis", "Lappula squarrosa")) %>%
  mutate(sp=str_replace(sp,"Lepidotheca suaveolens", "Matricaria discoidea")) %>%
  mutate(sp=str_replace(sp,"Listera cordata", "Neottia cordata")) %>%
  mutate(sp=str_replace(sp,"Leontodon autumnalis", "Scorzoneroides autumnalis")) %>%
  mutate(sp=str_replace(sp,"Loiseleuria procumbens", "Kalmia procumbens")) %>%
  mutate(sp=str_replace(sp,"Logfia arvensis", "Filago arvensis")) %>%
  mutate(sp=str_replace(sp,"Mentha _verticillata", "Mentha verticillata")) %>%
  mutate(sp=str_replace(sp,"Minuartia rubella", "Sabulina rubella")) %>%
  mutate(sp=str_replace(sp,"Minuartia stricta", "Sabulina stricta")) %>%
  mutate(sp=str_replace(sp,"Mycelis muralis", "Lactuca muralis")) %>%
  mutate(sp=str_replace(sp,"Omalotheca supina", "Gnaphalium supinum")) %>%
  mutate(sp=str_replace(sp,"Omalotheca norvegica", "Gnaphalium norvegicum")) %>%
  mutate(sp=str_replace(sp,"Omalotheca sylvatica", "Gnaphalium sylvaticum")) %>%
  mutate(sp=str_replace(sp,"Ononis arvensis", "Ononis spinosa")) %>%
  mutate(sp=str_replace(sp,"Oreopteris limbosperma", "Thelypteris limbosperma")) %>%
  mutate(sp=str_replace(sp,"Oxycoccus microcarpus", "Vaccinium microcarpum")) %>%
  mutate(sp=str_replace(sp,"Oxycoccus palustris", "Vaccinium oxycoccos")) %>%
  mutate(sp=str_replace(sp,"Phalaris minor", "Phalaris arundinacea")) %>%
  mutate(sp=str_replace(sp,"Phalaroides arundinacea", "Phalaris arundinacea")) %>%
  mutate(sp=str_replace(sp,"Pinus unicinata", "Pinus mugo")) %>%
  mutate(sp=str_replace(sp,"Platanthera montana", "Platanthera chlorantha")) %>%
  mutate(sp=str_replace(sp,"Poa alpigena", "Poa pratensis")) %>%
  mutate(sp=str_replace(sp,"Poa angustifolia", "Poa pratensis")) %>%
  mutate(sp=str_replace(sp,"Poa laxa", "Poa flexuosa")) %>%
  mutate(sp=str_replace(sp,"Poa _herjedalica", "Poa herjedalica")) %>%
  mutate(sp=str_replace(sp,"Poa _jemtlandica", "Poa jemtlandica")) %>%
  mutate(sp=str_replace(sp,"Poa lindebergii", "Poa arctica")) %>%
  mutate(sp=str_replace(sp,"Pyrola grandiflora", "Pyrola rotundifolia")) %>%
  mutate(sp=str_replace(sp,"Rhamnus catharticus", "Rhamnus cathartica")) %>%
  mutate(sp=str_replace(sp,"Rumex alpestris", "Rumex acetosa")) %>%
  mutate(sp=str_replace(sp,"Salix _fragilis", "Salix fragilis")) %>%
  mutate(sp=str_replace(sp,"Saxifraga _opdalensis", "Saxifraga opdalensis")) %>%
  mutate(sp=str_replace(sp,"Spergularia salina", "Spergularia marina")) %>%
  mutate(sp=str_replace(sp,"Syringa emodi", "Syringa vulgaris")) %>%
  mutate(sp=str_replace(sp,"Taraxacum crocea", "Taraxacum officinale")) %>%
  mutate(sp=str_replace(sp,"Taraxacum croceum", "Taraxacum officinale")) %>%
  mutate(sp=str_replace(sp,"Taraxacum erythrospermum", "Taraxacum officinale")) %>%
  mutate(sp=str_replace(sp,"Taraxacum hamatum", "Taraxacum officinale")) %>%
  mutate(sp=str_replace(sp,"Trientalis europaea", "Lysimachia europaea")) %>%
  mutate(sp=str_replace(sp,"Trifolium pallidum", "Trifolium pratense")) %>%
  mutate(sp=str_replace(sp,"Vicia orobus", "Vicia cassubica"))


## merge species data with indicators
NiN.sp.ind <- merge(NiN.sp,ind.dat, by.x="sp", by.y="species", all.x=T)
summary(NiN.sp.ind)

NiN.sp.ind[NiN.sp.ind==999] <- NA

# checking which species didn't find a match
unique(NiN.sp.ind[is.na(NiN.sp.ind$Moisture) & NiN.sp.ind$spgr %in% list("a1a","a1b","a1c"),'sp'])
# ok now

## matching with NiN ecosystem types - for wetlands
# NB! beware of rogue spaces in the 'Nature_type' & 'Sub_Type' variables, e.g. "Spring_Forest "
NiN.wetland <- NiN.sp.ind[,c("sp",paste(NiN.env[NiN.env$Nature_Type=="Mire","ID"]),colnames(ind.dat)[15:18])]   # Light, Moisture, Soil_reaction_pH, Nitrogen
NiN.wetland[1,]
names(NiN.wetland)

cbind(colnames(NiN.wetland),
      c("",
        'V3-C1a','V3-C1b','V3-C1c','V3-C1d','V3-C1e',
        'V1-C1a','V1-C1b','V1-C1c','V1-C1d','V1-C1e',
        'V1-C2a','V1-C2b','V1-C2c','V1-C2d',
        'V1-C3a','V1-C3b','V1-C3c','V1-C3d',
        'V1-C4a','V1-C4b','V1-C4c','V1-C4d',
        'V1-C4e','V1-C4f','V1-C4g','V1-C4h',
        'V3-C2','V1-C5',
        'V1-C6a','V1-C6b',
        'V1-C7a','V1-C7b',
        'V1-C8a','V1-C8b',
        'V2-C1a','V2-C1b',
        'V2-C2a','V2-C2b',
        'V2-C3a','V2-C3b',
        "V4-C2","V4-C3",
        "","",
        'V8-C1','V8-C2','V8-C3',
        rep("",10),
        rep("",4)
      )
)

NiN.wetland <- NiN.wetland[,c(1:43,46:48,59:62)]
colnames(NiN.wetland)[2:46] <- c('V3-C1a','V3-C1b','V3-C1c','V3-C1d','V3-C1e',
                              'V1-C1a','V1-C1b','V1-C1c','V1-C1d','V1-C1e',
                              'V1-C2a','V1-C2b','V1-C2c','V1-C2d',
                              'V1-C3a','V1-C3b','V1-C3c','V1-C3d',
                              'V1-C4a','V1-C4b','V1-C4c','V1-C4d',
                              'V1-C4e','V1-C4f','V1-C4g','V1-C4h',
                              'V3-C2','V1-C5',
                              'V1-C6a','V1-C6b',
                              'V1-C7a','V1-C7b',
                              'V1-C8a','V1-C8b',
                              'V2-C1a','V2-C1b',
                              'V2-C2a','V2-C2b',
                              'V2-C3a','V2-C3b',
                              'V4-C2','V4-C3',
                              'V8-C1','V8-C2','V8-C3'
)
head(NiN.wetland)


# translating the abundance classes into %-cover
coverscale <- data.frame(orig=0:6,
                         cov=c(0,1/32,1/8,3/8,0.6,4/5,1)
)

NiN.wetland.cov <- NiN.wetland
colnames(NiN.wetland.cov)
for (i in 2:46) {
  NiN.wetland.cov[,i] <- coverscale[,2][ match(NiN.wetland[,i], 0:6 ) ]
}

NiN.wetland.cov$sp <- as.factor(NiN.wetland.cov$sp)

This leaves us with the monitoring data including plant indicators (ANO.sp.ind) and the reference data including plant indicators (NiN.wetland.cov):

Table 6.1: ANO species occurence data set with attached plant trait data.
Species	ParentGlobalID	Light	Moisture	Soil_reaction_pH	Nitrogen
Abies alba	{CB1796B9-01F5-4109-B44E-4582CA855F93}	2	5	5	5
Abies alba	{AB9ED5C2-E906-4C73-B543-EC6CB28B39D5}	2	5	5	5
Abies alba	{32A9B462-5483-4D47-ADAF-78F11AF201AA}	2	5	5	5
Abies alba	{004C000D-459B-4244-96F4-4FF8B06454D4}	2	5	5	5
Abies alba	{B7DD61EE-A113-4486-A4B8-D50ACAAC648B}	2	5	5	5
Abies alba	{0431743B-F268-4819-98F7-FFB7006E55BA}	2	5	5	5

Table 6.2: Reference data set
sp	Light	Moisture	Soil_reaction_pH	Nitrogen
Abietinella abietina	NA	NA	NA	NA
Acer platanoides	4	4	6	5
Achillea millefolium	6	2	5	5
Achillea ptarmica	5	6	4	4
Acinos arvensis	7	1	7	3
Aconitum septentrionale	4	5	7	7

For each ecosystem type with a NiN species list, we can calculate a community weighted mean (CWM) for the relevant functional plant indicators. For wetland ecosystem we are testing “Light”, “Moisture”, “Soil_reaction_pH”, and “Nitrogen”. In order to get distributions of CWMs rather than one single value (for comparison with the empirical testing data) the NiN lists can be bootstrapped.

6.7.1.1.1 Bootstrap function for frequency abundance

function to calculate community weighted means of selected indicator values (ind)
for species lists (sp) with given abundances in percent (or on a scale from 0 to 1) in one or more ‘sites’ (abun)
with a given number of iterations (iter),
with species given a certain minimum abundance occurring in all bootstraps (obl), and
with a given re-sampling ratio of the original species list (rat)
in every bootstrap iteration the abundance of the sampled species can be randomly changed by a limited amount if wished by introducing a re-sampling of abundance values from adjacent abundance steps with a certain probability (var.abun)

indBoot.freq <- function(sp,abun,ind,iter,obl,rat=2/3,var.abun=F) {
  
  ind.b <- matrix(nrow=iter,ncol=length(colnames(abun)))
  colnames(ind.b) <- colnames(abun)
  ind.b <- as.data.frame(ind.b)  
  
  ind <- as.data.frame(ind)
  ind.list <- as.list(1:length(colnames(ind)))
  names(ind.list) <- colnames(ind)
  
  for (k in 1:length(colnames(ind)) ) {
    ind.list[[k]] <- ind.b }
  
  for (j in 1:length(colnames(abun)) ) {
    
    dat <- cbind(sp,abun[,j],ind)
    dat <- dat[dat[,2]>0,]            # only species that are present in the ecosystem
    dat <- dat[!is.na(dat[,3]),]      # only species that have indicator values
    
    for (i in 1:iter) {
      
      speciesSample <- sample(dat$sp[dat[,2] < obl], size=round( (length(dat$sp)-length(dat$sp[dat[,2]>=obl])) *rat,0), replace=F)  
      dat.b <- rbind(dat[dat[,2] >= obl,],
                     dat[match(speciesSample,dat$sp),]
      )
      
      if (var.abun==T) {
        for (m in 1:nrow(coverscale[-1,]) ) {
          xxx <- dat.b[dat.b[,2]==coverscale[-1,][m,2],2]
          if ( m==1 ) { dat.b[dat.b[,2]==coverscale[-1,][m,2],2] <- sample( c(0.01,coverscale[2:7,2]), prob = c(0.5, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0) ,size=length(xxx),replace=T) }
          if ( m==2 ) { dat.b[dat.b[,2]==coverscale[-1,][m,2],2] <- sample( c(0.01,coverscale[2:7,2]), prob = c(0.2, 0.3, 0.5, 0.0, 0.0, 0.0, 0.0) ,size=length(xxx),replace=T) }
          if ( m==3 ) { dat.b[dat.b[,2]==coverscale[-1,][m,2],2] <- sample( c(0.01,coverscale[2:7,2]), prob = c(0.0, 0.2, 0.3, 0.5, 0.0, 0.0, 0.0) ,size=length(xxx),replace=T) }
          if ( m==4 ) { dat.b[dat.b[,2]==coverscale[-1,][m,2],2] <- sample( c(0.01,coverscale[2:7,2]), prob = c(0.0, 0.0, 0.2, 0.3, 0.5, 0.0, 0.0) ,size=length(xxx),replace=T) }
          if ( m==5 ) { dat.b[dat.b[,2]==coverscale[-1,][m,2],2] <- sample( c(0.01,coverscale[2:7,2]), prob = c(0.0, 0.0, 0.0, 0.2, 0.3, 0.5, 0.0) ,size=length(xxx),replace=T) }
          if ( m==6 ) { dat.b[dat.b[,2]==coverscale[-1,][m,2],2] <- sample( c(0.01,coverscale[2:7,2]), prob = c(0.0, 0.0, 0.0, 0.0, 0.2, 0.3, 0.5) ,size=length(xxx),replace=T) }
        }
        dat.b[!is.na(dat.b[,2]) & dat.b[,2]<=(0),2] <- 0.01
        dat.b[!is.na(dat.b[,2]) & dat.b[,2]>1,2] <- 1
      }
      
      for (k in 1:length(colnames(ind))) {
        
        if ( nrow(dat.b)>2 ) {
          
          ind.b <- sum(dat.b[!is.na(dat.b[,2+k]),2] * dat.b[!is.na(dat.b[,2+k]),2+k] , na.rm=T) / sum(dat.b[!is.na(dat.b[,2+k]),2],na.rm=T)
          ind.list[[k]][i,j] <- ind.b
          
        } else {ind.list[[k]][i,j] <- NA}
        
      }
 
#      print(paste(i,"",j)) 
    }
    
  }
  return(ind.list)
}

colnames(NiN.wetland)
#>  [1] "sp"               "V3-C1a"          
#>  [3] "V3-C1b"           "V3-C1c"          
#>  [5] "V3-C1d"           "V3-C1e"          
#>  [7] "V1-C1a"           "V1-C1b"          
#>  [9] "V1-C1c"           "V1-C1d"          
#> [11] "V1-C1e"           "V1-C2a"          
#> [13] "V1-C2b"           "V1-C2c"          
#> [15] "V1-C2d"           "V1-C3a"          
#> [17] "V1-C3b"           "V1-C3c"          
#> [19] "V1-C3d"           "V1-C4a"          
#> [21] "V1-C4b"           "V1-C4c"          
#> [23] "V1-C4d"           "V1-C4e"          
#> [25] "V1-C4f"           "V1-C4g"          
#> [27] "V1-C4h"           "V3-C2"           
#> [29] "V1-C5"            "V1-C6a"          
#> [31] "V1-C6b"           "V1-C7a"          
#> [33] "V1-C7b"           "V1-C8a"          
#> [35] "V1-C8b"           "V2-C1a"          
#> [37] "V2-C1b"           "V2-C2a"          
#> [39] "V2-C2b"           "V2-C3a"          
#> [41] "V2-C3b"           "V4-C2"           
#> [43] "V4-C3"            "V8-C1"           
#> [45] "V8-C2"            "V8-C3"           
#> [47] "Light"            "Moisture"        
#> [49] "Soil_reaction_pH" "Nitrogen"
# 1st column is the species
# 2nd-46th column is the abundances of sp in different ecosystem types
# 47th-50th column is the indicator values of the respective species
# we choose 1000 iterations
# species with abundance 1 (i.e. a max of 100%, must be included in each sample)
# each sample re-samples 2/3 of the number of species
# the abundance of the re-sampled species may vary (see bootstrap function for details)

Running the bootstraps:

wetland.ref.cov <- indBoot.freq(sp=NiN.wetland.cov[,1],abun=NiN.wetland.cov[,2:46],ind=NiN.wetland.cov[,47:50],
                          iter=1000,obl=1,rat=2/3,var.abun=T)

# fixing NaNs
for (i in 1:length(wetland.ref.cov) ) {
  for (j in 1:ncol(wetland.ref.cov[[i]]) ) {
    v <- wetland.ref.cov[[i]][,j]
    v[is.nan(v)] <- NA
    wetland.ref.cov[[i]][,j] <- v
  }
}

head(wetland.ref.cov[[1]]) %>%
    kable("html", caption = "Table showing the first 6 rows of the boostrapped data set.") %>% kable_styling("striped") %>% scroll_box(width = "100%")

Table 6.3: Table showing the first 6 rows of the boostrapped data set.
V3-C1a	V3-C1b	V3-C1c	V3-C1d	V3-C1e	V1-C1a	V1-C1b	V1-C1c	V1-C1d	V1-C1e	V1-C2a	V1-C2b	V1-C2c	V1-C2d	V1-C3a	V1-C3b	V1-C3c	V1-C3d	V1-C4a	V1-C4b	V1-C4c	V1-C4d	V1-C4e	V1-C4f	V1-C4g	V1-C4h	V3-C2	V1-C5	V1-C6a	V1-C6b	V1-C7a	V1-C7b	V1-C8a	V1-C8b	V2-C1a	V2-C1b	V2-C2a	V2-C2b	V2-C3a	V2-C3b	V4-C2	V4-C3	V8-C1	V8-C2	V8-C3
6.336842	5.702703	5.891986	5.407407	5.284661	5.682713	5.723592	5.291436	5.353811	5.171986	5.823009	5.360731	5.617886	5.255973	5.879501	5.451771	5.150421	5.454386	6.036180	5.361275	5.287129	5.246544	5.826291	6.096953	5.435681	4.918367	4.994775	4.727076	4.996255	4.630435	5.005076	4.927380	4.901750	4.761329	4.395683	4.030844	4.660793	3.813056	4.138365	3.871470	4.859013	4.983607	4.767647	4.595033	4.620053
6.162896	5.603411	5.233202	5.353871	5.077465	5.908840	5.667160	5.929448	5.263294	5.245902	6.299618	5.416517	5.716802	5.514939	6.344234	5.197347	5.304290	5.319415	5.798913	5.616487	5.851117	5.579498	6.095665	5.593711	5.651720	5.438084	5.015432	4.652138	5.033333	4.773573	4.801927	4.846557	5.127968	4.896014	4.188559	4.155105	4.615487	4.059268	4.056991	3.942257	5.042409	5.249798	4.660854	4.547687	4.647333
5.807692	6.099119	5.567628	5.296296	5.269912	6.087065	6.065979	5.575145	5.536041	5.169304	5.894265	5.898420	5.455844	5.415879	5.713128	5.320140	5.617048	5.505593	6.293651	5.648649	5.619326	5.577703	5.855586	5.934783	5.580581	5.851128	4.773709	4.726855	4.904824	4.773680	4.572449	4.504233	5.070941	4.808081	4.110694	4.091158	4.479858	3.668596	4.187162	4.057981	5.047048	4.992795	4.679577	4.609576	4.564313
6.775086	5.789572	5.493739	5.475034	5.269912	6.352518	5.483051	5.570014	5.728624	5.022152	5.849462	5.354067	5.591304	5.561783	5.722591	5.573379	5.204703	5.590814	5.889374	5.685484	5.371715	5.093677	5.851802	5.616374	5.151030	5.386395	4.870244	4.667209	4.994338	4.894553	4.884762	4.812344	4.881336	5.094451	4.260394	4.165308	4.733990	3.914540	4.233787	4.077534	4.738647	5.027458	4.835616	4.512004	4.539221
5.903114	5.789572	5.641694	5.492308	5.007634	6.364508	5.517544	5.720264	5.458237	5.034188	5.867383	5.418605	5.108579	5.505747	5.830116	5.317019	5.466457	5.370787	5.776280	5.713755	5.469660	5.717724	5.981818	5.697906	5.566054	5.235374	4.886615	4.814668	4.880233	4.836431	4.783305	4.893412	4.988858	4.848684	4.273535	4.234497	4.586890	3.797063	4.404070	3.836165	4.844905	4.734485	4.500463	4.502172	4.762053
6.336842	5.726054	5.668407	5.569579	5.178284	5.616307	5.634868	5.710968	5.466031	4.925319	5.891344	5.359712	5.739508	5.353218	5.906937	5.390173	5.092166	5.493927	6.177054	5.651599	5.235872	5.522613	5.779343	5.812823	5.526852	5.490647	4.876712	4.749831	5.161597	4.823015	4.814433	4.869048	4.859482	4.702290	4.174721	4.005499	4.520000	4.033181	4.377852	3.859554	4.999081	4.966232	4.805742	4.491489	4.608948

This results in an R-list with a slot for every selected indicator, and in every slot there’s a data frame with as many columns as there are NiN species lists and as many rows as there were iterations in the bootstrap.

Next, we need to derive scaling values from these bootstrap-lists (the columns) for every mapping unit in NiN. Here, we define things in the following way:

Median = reference values
0.025 and 0.975 quantiles = lower and upper limit values
min and max of the respective indicator’s scale = min/max values


# every NiN-type is represented by one 'generalisert artsliste'
# some NiN-types are represented by two such species lists
# in some cases two NiN-types are represented by the same species list
head(wetland.ref.cov[[1]])
wetland.ref.cov[[1]][0,]

# checking the actual NiN-types in the wetland lists
wetland.NiNtypes <- colnames(wetland.ref.cov[["Light"]])
wetland.NiNtypes <- substr(wetland.NiNtypes,1,5)
unique(wetland.NiNtypes)

# 4 indicator-value indicators: Tyler's Light, Moisture, Soil_reaction_pH, "Nitrogen"
indEll.n=4
# creating a table to hold:
# Tyler: the 0.5 quantile (median), 0.05 quantile and  0.95 quantile for each NiN-type
# for every nature type (nrows)
tab <- matrix(ncol=3*indEll.n, nrow=length(unique(wetland.NiNtypes)) ) # 43 basic ecosystem types
# coercing the values into the table
# NiN-types where each type is represented by one species list (including when one species list represents two NiN-types)
names(wetland.ref.cov[["Light"]])
x <- c(27,28,41:45)

for (i in 1:length(x) ) {
  tab[i,1:3] <- quantile(as.matrix(wetland.ref.cov[["Light"]][,x[i]]),probs=c(0.025,0.5,0.975),na.rm=T)
  tab[i,4:6] <- quantile(as.matrix(wetland.ref.cov[["Moisture"]][,x[i]]),probs=c(0.025,0.5,0.975),na.rm=T)
  tab[i,7:9] <- quantile(as.matrix(wetland.ref.cov[["Soil_reaction_pH"]][,x[i]]),probs=c(0.025,0.5,0.975),na.rm=T)
  tab[i,10:12] <- quantile(as.matrix(wetland.ref.cov[["Nitrogen"]][,x[i]]),probs=c(0.025,0.5,0.975),na.rm=T)
}

tab <- as.data.frame(tab)
tab$NiN <- NA
tab$NiN[1:length(x)] <- names(wetland.ref.cov[[1]])[x]

# NiN-types represented by several species lists
wetland.NiNtypes2 <- wetland.NiNtypes[-x]
unique(wetland.NiNtypes2)
grep(pattern=unique(wetland.NiNtypes2)[1], x=wetland.NiNtypes) # finds columns in e.g. colnames(wetland.ref.cov[["Continentality"]]) that match the first NiN-type


for (i in 1:length(unique(wetland.NiNtypes2)) ) {
  tab[length(x)+i,1:3] <- quantile(as.matrix(wetland.ref.cov[["Light"]][,grep(pattern=unique(wetland.NiNtypes2)[i], x=wetland.NiNtypes)]),probs=c(0.025,0.5,0.975),na.rm=T)
  tab[length(x)+i,4:6] <- quantile(as.matrix(wetland.ref.cov[["Moisture"]][,grep(pattern=unique(wetland.NiNtypes2)[i], x=wetland.NiNtypes)]),probs=c(0.025,0.5,0.975),na.rm=T)
  tab[length(x)+i,7:9] <- quantile(as.matrix(wetland.ref.cov[["Soil_reaction_pH"]][,grep(pattern=unique(wetland.NiNtypes2)[i], x=wetland.NiNtypes)]),probs=c(0.025,0.5,0.975),na.rm=T)
  tab[length(x)+i,10:12] <- quantile(as.matrix(wetland.ref.cov[["Nitrogen"]][,grep(pattern=unique(wetland.NiNtypes2)[i], x=wetland.NiNtypes)]),probs=c(0.025,0.5,0.975),na.rm=T)
  tab$NiN[length(x)+i] <- unique(wetland.NiNtypes2)[i]
  
}

tab

# making it a proper data frame
round(tab[,1:12],digits=2)

colnames(tab) <- c("Light_q2.5","Light_q50","Light_q97.5",
                   "Moist_q2.5","Moist_q50","Moist_q97.5",
                   "pH_q2.5","pH_q50","pH_q97.5",
                   "Nitrogen_q2.5","Nitrogen_q50","Nitrogen_q97.5",
                   "NiN")
summary(tab)
tab$NiN <- gsub("C", "C-", tab$NiN) # add extra hyphon after C for NiN-types
tab



# restructuring into separate indicators for lower (q2.5) and higher (q97.5) than reference value (=median, q50)
y.Light <- numeric(length=nrow(tab)*2)
y.Light[((1:dim(tab)[1])*2)-1] <- tab$Light_q2.5 
y.Light[((1:dim(tab)[1])*2)] <- tab$Light_q97.5 

y.Moist <- numeric(length=nrow(tab)*2)
y.Moist[((1:dim(tab)[1])*2)-1] <- tab$Moist_q2.5 
y.Moist[((1:dim(tab)[1])*2)] <- tab$Moist_q97.5 

y.pH <- numeric(length=nrow(tab)*2)
y.pH[((1:dim(tab)[1])*2)-1] <- tab$pH_q2.5 
y.pH[((1:dim(tab)[1])*2)] <- tab$pH_q97.5 

y.Nitrogen <- numeric(length=nrow(tab)*2)
y.Nitrogen[((1:dim(tab)[1])*2)-1] <- tab$Nitrogen_q2.5 
y.Nitrogen[((1:dim(tab)[1])*2)] <- tab$Nitrogen_q97.5 

# creating final objects holding the reference and limit values for all indicators

# ref object for indicators
wetland.ref.cov.val <- data.frame(N1=rep('wetland',(nrow(tab)*2*indEll.n)),
                              hoved=c(rep('NA',(nrow(tab)*2*indEll.n))),
                              grunn=c(rep(rep(tab$NiN,each=2),indEll.n)),
                              county=rep('all',(nrow(tab)*2*indEll.n)),
                              region=rep('all',(nrow(tab)*2*indEll.n)),
                              Ind=c(rep(c('Light1','Light2'),nrow(tab)),
                                    rep(c('Moist1','Moist2'),nrow(tab)),
                                    rep(c('pH1','pH2'),nrow(tab)),
                                    rep(c('Nitrogen1','Nitrogen2'),nrow(tab))
                              ),
                              Rv=c(rep(tab$Light_q50,each=2),
                                   rep(tab$Moist_q50,each=2),
                                   rep(tab$pH_q50,each=2),
                                   rep(tab$Nitrogen_q50,each=2)
                              ),
                              Gv=c(y.Light,y.Moist,y.pH,y.Nitrogen),
                              maxmin=c(rep(c(1,7),nrow(tab)), # 7 levels of light
                                       rep(c(1,12),nrow(tab)), # 12 levels of moisture
                                       rep(c(1,8),nrow(tab)), # 8 levels of soil reaction pH
                                       rep(c(1,9),nrow(tab))  # 9 levels of nitrogen
                              )
)

wetland.ref.cov.val
wetland.ref.cov.val$grunn <- as.factor(wetland.ref.cov.val$grunn)
wetland.ref.cov.val$Ind <- as.factor(wetland.ref.cov.val$Ind)
summary(wetland.ref.cov.val)

wetland.ref.cov.val %>%
      kable("html", caption = "Reference values (Rv and threshold values (Gv) for each indicator and nature type combination") %>% kable_styling("striped") %>% scroll_box(width = "100%", height = "300px")

Table 6.4: Reference values (Rv and threshold values (Gv) for each indicator and nature type combination
N1	hoved	grunn	county	region	Ind	Rv	Gv	maxmin
wetland	NA	V3-C-2	all	all	Light1	4.900265	4.736150	1
wetland	NA	V3-C-2	all	all	Light2	4.900265	5.120000	7
wetland	NA	V1-C-5	all	all	Light1	4.760764	4.597880	1
wetland	NA	V1-C-5	all	all	Light2	4.760764	4.902329	7
wetland	NA	V4-C-2	all	all	Light1	4.901870	4.610423	1
wetland	NA	V4-C-2	all	all	Light2	4.901870	5.232242	7
wetland	NA	V4-C-3	all	all	Light1	4.996095	4.677531	1
wetland	NA	V4-C-3	all	all	Light2	4.996095	5.309336	7
wetland	NA	V8-C-1	all	all	Light1	4.709785	4.464570	1
wetland	NA	V8-C-1	all	all	Light2	4.709785	4.972819	7
wetland	NA	V8-C-2	all	all	Light1	4.586035	4.425726	1
wetland	NA	V8-C-2	all	all	Light2	4.586035	4.779277	7
wetland	NA	V8-C-3	all	all	Light1	4.635103	4.490343	1
wetland	NA	V8-C-3	all	all	Light2	4.635103	4.790331	7
wetland	NA	V3-C-1	all	all	Light1	5.594595	4.892675	1
wetland	NA	V3-C-1	all	all	Light2	5.594595	6.567474	7
wetland	NA	V1-C-1	all	all	Light1	5.575057	4.958721	1
wetland	NA	V1-C-1	all	all	Light2	5.575057	6.315900	7
wetland	NA	V1-C-2	all	all	Light1	5.522530	5.139566	1
wetland	NA	V1-C-2	all	all	Light2	5.522530	6.264477	7
wetland	NA	V1-C-3	all	all	Light1	5.450281	5.089262	1
wetland	NA	V1-C-3	all	all	Light2	5.450281	6.152334	7
wetland	NA	V1-C-4	all	all	Light1	5.647878	5.126794	1
wetland	NA	V1-C-4	all	all	Light2	5.647878	6.162793	7
wetland	NA	V1-C-6	all	all	Light1	4.816224	4.525242	1
wetland	NA	V1-C-6	all	all	Light2	4.816224	5.124243	7
wetland	NA	V1-C-7	all	all	Light1	4.804564	4.504224	1
wetland	NA	V1-C-7	all	all	Light2	4.804564	5.059866	7
wetland	NA	V1-C-8	all	all	Light1	4.915094	4.558245	1
wetland	NA	V1-C-8	all	all	Light2	4.915094	5.303735	7
wetland	NA	V2-C-1	all	all	Light1	4.173047	3.974819	1
wetland	NA	V2-C-1	all	all	Light2	4.173047	4.598703	7
wetland	NA	V2-C-2	all	all	Light1	4.256174	3.685444	1
wetland	NA	V2-C-2	all	all	Light2	4.256174	4.700616	7
wetland	NA	V2-C-3	all	all	Light1	4.125515	3.631843	1
wetland	NA	V2-C-3	all	all	Light2	4.125515	4.457249	7
wetland	NA	V3-C-2	all	all	Moist1	5.736585	5.201878	1
wetland	NA	V3-C-2	all	all	Moist2	5.736585	6.151852	12
wetland	NA	V1-C-5	all	all	Moist1	5.821132	5.391474	1
wetland	NA	V1-C-5	all	all	Moist2	5.821132	6.265910	12
wetland	NA	V4-C-2	all	all	Moist1	7.088321	6.708205	1
wetland	NA	V4-C-2	all	all	Moist2	7.088321	7.496453	12
wetland	NA	V4-C-3	all	all	Moist1	6.935913	6.599488	1
wetland	NA	V4-C-3	all	all	Moist2	6.935913	7.284055	12
wetland	NA	V8-C-1	all	all	Moist1	7.217742	6.803539	1
wetland	NA	V8-C-1	all	all	Moist2	7.217742	7.663696	12
wetland	NA	V8-C-2	all	all	Moist1	7.251258	6.925680	1
wetland	NA	V8-C-2	all	all	Moist2	7.251258	7.577229	12
wetland	NA	V8-C-3	all	all	Moist1	7.382590	7.024416	1
wetland	NA	V8-C-3	all	all	Moist2	7.382590	7.755543	12
wetland	NA	V3-C-1	all	all	Moist1	6.907166	5.528727	1
wetland	NA	V3-C-1	all	all	Moist2	6.907166	8.486166	12
wetland	NA	V1-C-1	all	all	Moist1	7.280132	5.735756	1
wetland	NA	V1-C-1	all	all	Moist2	7.280132	8.718040	12
wetland	NA	V1-C-2	all	all	Moist1	7.915167	6.426842	1
wetland	NA	V1-C-2	all	all	Moist2	7.915167	8.924438	12
wetland	NA	V1-C-3	all	all	Moist1	8.216718	7.113426	1
wetland	NA	V1-C-3	all	all	Moist2	8.216718	9.016390	12
wetland	NA	V1-C-4	all	all	Moist1	8.040734	7.065645	1
wetland	NA	V1-C-4	all	all	Moist2	8.040734	8.935004	12
wetland	NA	V1-C-6	all	all	Moist1	7.233970	6.016056	1
wetland	NA	V1-C-6	all	all	Moist2	7.233970	8.688171	12
wetland	NA	V1-C-7	all	all	Moist1	7.271581	6.326991	1
wetland	NA	V1-C-7	all	all	Moist2	7.271581	8.463120	12
wetland	NA	V1-C-8	all	all	Moist1	7.091729	6.258498	1
wetland	NA	V1-C-8	all	all	Moist2	7.091729	7.784715	12
wetland	NA	V2-C-1	all	all	Moist1	5.859056	5.112448	1
wetland	NA	V2-C-1	all	all	Moist2	5.859056	8.126484	12
wetland	NA	V2-C-2	all	all	Moist1	6.488092	5.513549	1
wetland	NA	V2-C-2	all	all	Moist2	6.488092	7.920221	12
wetland	NA	V2-C-3	all	all	Moist1	6.772113	6.004292	1
wetland	NA	V2-C-3	all	all	Moist2	6.772113	7.759210	12
wetland	NA	V3-C-2	all	all	pH1	2.243827	1.831766	1
wetland	NA	V3-C-2	all	all	pH2	2.243827	2.533659	8
wetland	NA	V1-C-5	all	all	pH1	2.516520	2.155809	1
wetland	NA	V1-C-5	all	all	pH2	2.516520	2.836431	8
wetland	NA	V4-C-2	all	all	pH1	3.787654	3.437916	1
wetland	NA	V4-C-2	all	all	pH2	3.787654	4.212773	8
wetland	NA	V4-C-3	all	all	pH1	4.388411	3.953983	1
wetland	NA	V4-C-3	all	all	pH2	4.388411	4.813829	8
wetland	NA	V8-C-1	all	all	pH1	4.163711	3.889345	1
wetland	NA	V8-C-1	all	all	pH2	4.163711	4.486426	8
wetland	NA	V8-C-2	all	all	pH1	4.575937	4.381745	1
wetland	NA	V8-C-2	all	all	pH2	4.575937	4.795269	8
wetland	NA	V8-C-3	all	all	pH1	4.876376	4.682191	1
wetland	NA	V8-C-3	all	all	pH2	4.876376	5.070794	8
wetland	NA	V3-C-1	all	all	pH1	2.013841	1.494071	1
wetland	NA	V3-C-1	all	all	pH2	2.013841	2.868778	8
wetland	NA	V1-C-1	all	all	pH1	2.099490	1.697084	1
wetland	NA	V1-C-1	all	all	pH2	2.099490	3.207349	8
wetland	NA	V1-C-2	all	all	pH1	2.837740	2.076494	1
wetland	NA	V1-C-2	all	all	pH2	2.837740	3.691609	8
wetland	NA	V1-C-3	all	all	pH1	3.604344	2.768193	1
wetland	NA	V1-C-3	all	all	pH2	3.604344	4.317388	8
wetland	NA	V1-C-4	all	all	pH1	4.508963	3.410753	1
wetland	NA	V1-C-4	all	all	pH2	4.508963	5.601644	8
wetland	NA	V1-C-6	all	all	pH1	3.091580	2.655389	1
wetland	NA	V1-C-6	all	all	pH2	3.091580	3.699116	8
wetland	NA	V1-C-7	all	all	pH1	3.547120	3.246374	1
wetland	NA	V1-C-7	all	all	pH2	3.547120	3.916825	8
wetland	NA	V1-C-8	all	all	pH1	4.212409	3.808178	1
wetland	NA	V1-C-8	all	all	pH2	4.212409	4.677554	8
wetland	NA	V2-C-1	all	all	pH1	2.942611	2.628539	1
wetland	NA	V2-C-1	all	all	pH2	2.942611	3.685835	8
wetland	NA	V2-C-2	all	all	pH1	3.618644	3.196779	1
wetland	NA	V2-C-2	all	all	pH2	3.618644	4.074977	8
wetland	NA	V2-C-3	all	all	pH1	4.097649	3.783875	1
wetland	NA	V2-C-3	all	all	pH2	4.097649	4.434308	8
wetland	NA	V3-C-2	all	all	Nitrogen1	1.952862	1.811728	1
wetland	NA	V3-C-2	all	all	Nitrogen2	1.952862	2.142723	9
wetland	NA	V1-C-5	all	all	Nitrogen1	2.167561	1.853110	1
wetland	NA	V1-C-5	all	all	Nitrogen2	2.167561	2.438285	9
wetland	NA	V4-C-2	all	all	Nitrogen1	3.323729	3.008427	1
wetland	NA	V4-C-2	all	all	Nitrogen2	3.323729	3.678983	9
wetland	NA	V4-C-3	all	all	Nitrogen1	3.807000	3.262604	1
wetland	NA	V4-C-3	all	all	Nitrogen2	3.807000	4.234730	9
wetland	NA	V8-C-1	all	all	Nitrogen1	4.462636	4.148721	1
wetland	NA	V8-C-1	all	all	Nitrogen2	4.462636	4.812027	9
wetland	NA	V8-C-2	all	all	Nitrogen1	5.459677	5.142488	1
wetland	NA	V8-C-2	all	all	Nitrogen2	5.459677	5.816705	9
wetland	NA	V8-C-3	all	all	Nitrogen1	5.900868	5.629957	1
wetland	NA	V8-C-3	all	all	Nitrogen2	5.900868	6.152190	9
wetland	NA	V3-C-1	all	all	Nitrogen1	1.494737	1.072464	1
wetland	NA	V3-C-1	all	all	Nitrogen2	1.494737	1.996790	9
wetland	NA	V1-C-1	all	all	Nitrogen1	1.667114	1.389752	1
wetland	NA	V1-C-1	all	all	Nitrogen2	1.667114	2.134328	9
wetland	NA	V1-C-2	all	all	Nitrogen1	1.974614	1.567785	1
wetland	NA	V1-C-2	all	all	Nitrogen2	1.974614	2.519928	9
wetland	NA	V1-C-3	all	all	Nitrogen1	2.396491	1.746633	1
wetland	NA	V1-C-3	all	all	Nitrogen2	2.396491	3.048414	9
wetland	NA	V1-C-4	all	all	Nitrogen1	2.685903	2.053753	1
wetland	NA	V1-C-4	all	all	Nitrogen2	2.685903	3.184157	9
wetland	NA	V1-C-6	all	all	Nitrogen1	2.761761	2.303992	1
wetland	NA	V1-C-6	all	all	Nitrogen2	2.761761	3.191950	9
wetland	NA	V1-C-7	all	all	Nitrogen1	3.186189	2.846378	1
wetland	NA	V1-C-7	all	all	Nitrogen2	3.186189	3.487493	9
wetland	NA	V1-C-8	all	all	Nitrogen1	3.517028	3.058802	1
wetland	NA	V1-C-8	all	all	Nitrogen2	3.517028	4.025923	9
wetland	NA	V2-C-1	all	all	Nitrogen1	3.211879	2.899263	1
wetland	NA	V2-C-1	all	all	Nitrogen2	3.211879	3.767787	9
wetland	NA	V2-C-2	all	all	Nitrogen1	4.032121	3.605832	1
wetland	NA	V2-C-2	all	all	Nitrogen2	4.032121	4.559971	9
wetland	NA	V2-C-3	all	all	Nitrogen1	4.900636	4.489434	1
wetland	NA	V2-C-3	all	all	Nitrogen2	4.900636	5.338214	9

Once test data (ANO) and the scaling values from the reference data (NiN) are in place, we can calculate CWMs of the selected indicators for the ANO community data and scale them against the scaling values from the reference distribution. Note that we scale each ANO plot’s CWM against either the lower threshold value and the min value OR the upper threshold value and the max value based on whether the CWM is smaller or higher than the reference value. Since the scaled values for both sides range between 0 and 1, we generate separate lower and upper indicators for each plant functional indicator type. An ANO plot can only have a scaled value in either the lower or the upper indicator (the other one will be ‘NA’), except for the unlikely event that the CWM exactly matches the reference value, in which case both lower and upper indicator will receive a scaled indicator value of 1.

Here is the scaling function:


#### scaled values ####
r.s <- 1    # reference value
l.s <- 0.6  # limit value
a.s <- 0    # abscence of indicator, or indicator at maximum

#### function for calculating scaled values for measured value ####

## scaling function including truncation
scal <- function() {
  # place to hold the result
   x <- numeric()
  if (maxmin < ref) {
    # values >= the reference value equal 1
    if (val >= ref) {x <- 1}
    # values < the reference value and >= the limit value can be deducted from the linear relationship between these two
    if (val < ref & val >= lim) {x <- (l.s + (val-lim) * ( (r.s-l.s) / (ref-lim) ) )}
    # values < the limit value and > maxmin can be deducted from the linear relationship between these two
    if (val < lim & val > maxmin) {x <- (a.s + (val-maxmin) * ( (l.s-a.s) / (lim-maxmin) ) )}
    # value equals or lower than maxmin
    if (val <= maxmin) {x <-0}
  } else {
    # values <= the reference value equal 1
    if (val <= ref) {x <- 1}
    # values > the reference value and <= the limit value can be deducted from the linear relationship between these two
    if (val > ref & val <= lim) {x <- ( r.s - ( (r.s - l.s) * (val - ref) / (lim - ref) ) )}
    # values > the limit value and < maxmin can be deducted from the linear relationship between these two
    if (val > lim) {x <- ( l.s - (l.s * (val - lim) / (maxmin - lim) ) )}
    # value equals or larger than maxmin
    if (val >= maxmin) {x <-0}
  }
  return(x)
  
}

We then can prepare a list of data frames to hold the results and perform the scaling according to the principles described in NINA report 1967 (Töpper and Jakobsson 2021).


unique(ANO.geo$hovedtype_rute) # NiN types in data
unique(substr(wetland.ref.cov.val$grunn,1,2)) # NiN types in reference
#### creating dataframe to hold the results for wetlands ####
# all ANO points
nrow(ANO.geo)
# all wetland ANO points
nrow(ANO.geo[ANO.geo$hovedtype_rute %in% list("V1","V2","V3","V4","V5","V6","V7","V8","V9","V10","V11","V12","V13"),])
# all wetland ANO points with a NiN-type represented in the reference
nrow(ANO.geo[ANO.geo$hovedtype_rute %in% unique(substr(wetland.ref.cov.val$grunn,1,2)),])
# ok, we'll be losing 70 (our of 1349) that are not covered by the reference
ANO.wetland <- ANO.geo[ANO.geo$hovedtype_rute %in% list("V1","V2","V3","V4","V5","V6","V7","V8","V9","V10","V11","V12","V13"),]

head(ANO.wetland)
# update row-numbers
row.names(ANO.wetland) <- 1:nrow(ANO.wetland)
head(ANO.wetland)
dim(ANO.wetland)
colnames(ANO.wetland)

length(levels(as.factor(ANO.wetland$ano_flate_id)))
length(levels(as.factor(ANO.wetland$ano_punkt_id)))
summary(as.factor(ANO.wetland$ano_punkt_id))
# four points that are double
ANO.wetland[ANO.wetland$ano_punkt_id=="ANO0159_55",] # double registration, said so in comment. -> choose row 207 over 206
ANO.wetland <- ANO.wetland[-206,]
row.names(ANO.wetland) <- 1:nrow(ANO.wetland) # update row-numbers
ANO.wetland[ANO.wetland$ano_punkt_id=="ANO0283_22",] # 2019 & 2021. Lot of NA's in 2019 -> omit 2019
ANO.wetland <- ANO.wetland[-156,]
row.names(ANO.wetland) <- 1:nrow(ANO.wetland) # update row-numbers
ANO.wetland[ANO.wetland$ano_punkt_id=="ANO0363_24",]
ANO.wetland[ANO.wetland$ano_flate_id=="ANO0363","ano_punkt_id"] # point-ID 15 is missing, but 24 is double. Likely that registrations are valid, but wrong point-ID.  -> keep both, call the second obs the one that's missing
ANO.wetland[311,"ano_punkt_id"] <- "ANO0363_15"
ANO.wetland[ANO.wetland$ano_punkt_id=="ANO1550_64",] # point-ID 66 is missing, but 64 is double. Likely that registrations are valid, but wrong point-ID.  -> keep both
ANO.wetland[ANO.wetland$ano_flate_id=="ANO1550","ano_punkt_id"] # point-ID 66 is missing, but 64 is double. Likely that registrations are valid, but wrong point-ID.  -> keep both
ANO.wetland[1273,"ano_punkt_id"] <- "ANO1550_66"

unique(ANO.wetland$hovedoekosystem_punkt)
unique(ANO.wetland$hovedtype_rute)
unique(ANO.wetland$kartleggingsenhet_1m2)
ANO.wetland$hovedtype_rute <- factor(ANO.wetland$hovedtype_rute)
ANO.wetland$kartleggingsenhet_1m2 <- factor(ANO.wetland$kartleggingsenhet_1m2)
summary(ANO.wetland$Hovedtype_rute)
summary(ANO.wetland$Kartleggingsenhet_rute)

results.wetland <- list()
ind <- unique(wetland.ref.cov.val$Ind)
# choose columns for site description
colnames(ANO.wetland)
results.wetland[['original']] <- ANO.wetland
# drop geometry
st_geometry(results.wetland[['original']]) <- NULL
results.wetland[['original']]

# add columns for indicators
nvar.site <- ncol(results.wetland[['original']])
for (i in 1:length(ind) ) {results.wetland[['original']][,i+nvar.site] <- NA}
colnames(results.wetland[['original']])[(nvar.site+1):(length(ind)+nvar.site)] <- paste(ind)
for (i in (nvar.site+1):(length(ind)+nvar.site) ) {results.wetland[['original']][,i] <- as.numeric(results.wetland[['original']][,i])}
summary(results.wetland[['original']])
#results.wetland[['original']]$Region <- as.factor(results.wetland[['original']]$Region)
results.wetland[['original']]$GlobalID <- as.factor(results.wetland[['original']]$GlobalID)
results.wetland[['original']]$ano_flate_id <- as.factor(results.wetland[['original']]$ano_flate_id)
results.wetland[['original']]$ano_punkt_id <- as.factor(results.wetland[['original']]$ano_punkt_id)
results.wetland[['original']]$hovedoekosystem_punkt <- as.factor(results.wetland[['original']]$hovedoekosystem_punkt)
#results.wetland[['original']]$Hovedoekosystem_rute  <- as.factor(results.wetland[['original']]$Hovedoekosystem_rute )
results.wetland[['original']]$kartleggingsenhet_1m2 <- as.factor(results.wetland[['original']]$kartleggingsenhet_1m2)
results.wetland[['original']]$hovedtype_rute    <- as.factor(results.wetland[['original']]$hovedtype_rute)


# roll out
results.wetland[['scaled']] <- results.wetland[['non-truncated']] <- results.wetland[['original']]


#### calculating scaled and non-truncated values for the indicators based on the dataset ####
for (i in 1:nrow(ANO.wetland) ) {  #
  tryCatch({
    print(i)
    print(paste(ANO.wetland$ano_flate_id[i]))
    print(paste(ANO.wetland$ano_punkt_id[i]))
#    ANO.wetland$Hovedoekosystem_sirkel[i]
#    ANO.wetland$Hovedoekosystem_rute[i]

    # if the ANO.hovedtype exists in the reference
    if (ANO.wetland$hovedtype_rute[i] %in% unique(substr(wetland.ref.cov.val$grunn,1,2)) ) {
      
      # if there is any species present in current ANO point  
      if ( length(ANO.sp.ind[ANO.sp.ind$ParentGlobalID==as.character(ANO.wetland$GlobalID[i]),'Species']) > 0 ) {
        

          
          # Light
          dat <- ANO.sp.ind[ANO.sp.ind$ParentGlobalID==as.character(ANO.wetland$GlobalID[i]),c('art_dekning','Light')]
          results.wetland[['original']][i,'richness'] <- nrow(dat)
          dat <- dat[!is.na(dat$Light),]
          
          if ( nrow(dat)>0 ) {
            
            val <- sum(dat[,'art_dekning'] * dat[,'Light'],na.rm=T) / sum(dat[,'art_dekning'],na.rm=T)
            # lower part of distribution
            ref <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='Light1' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'Rv']
            lim <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='Light1' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'Gv']
            maxmin <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='Light1' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'maxmin']
            # coercing x into results.wetland dataframe
            results.wetland[['scaled']][i,'Light1'] <- scal() 
            results.wetland[['non-truncated']][i,'Light1'] <- scal.2() 
            results.wetland[['original']][i,'Light1'] <- val 
            
            # upper part of distribution
            ref <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='Light2' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'Rv']
            lim <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='Light2' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'Gv']
            maxmin <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='Light2' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'maxmin']
            # coercing x into results.wetland dataframe
            results.wetland[['scaled']][i,'Light2'] <- scal() 
            results.wetland[['non-truncated']][i,'Light2'] <- scal.2() 
            results.wetland[['original']][i,'Light2'] <- val
          }
          
          
          # Moisture
          dat <- ANO.sp.ind[ANO.sp.ind$ParentGlobalID==as.character(ANO.wetland$GlobalID[i]),c('art_dekning','Moisture')]
          results.wetland[['original']][i,'richness'] <- nrow(dat)
          dat <- dat[!is.na(dat$Moisture),]
          
          if ( nrow(dat)>0 ) {
            
            val <- sum(dat[,'art_dekning'] * dat[,'Moisture'],na.rm=T) / sum(dat[,'art_dekning'],na.rm=T)
            # lower part of distribution
            ref <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='Moist1' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'Rv']
            lim <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='Moist1' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'Gv']
            maxmin <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='Moist1' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'maxmin']
            # coercing x into results.wetland dataframe
            results.wetland[['scaled']][i,'Moist1'] <- scal() 
            results.wetland[['non-truncated']][i,'Moist1'] <- scal.2() 
            results.wetland[['original']][i,'Moist1'] <- val 
            
            # upper part of distribution
            ref <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='Moist2' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'Rv']
            lim <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='Moist2' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'Gv']
            maxmin <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='Moist2' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'maxmin']
            # coercing x into results.wetland dataframe
            results.wetland[['scaled']][i,'Moist2'] <- scal() 
            results.wetland[['non-truncated']][i,'Moist2'] <- scal.2() 
            results.wetland[['original']][i,'Moist2'] <- val
          }
          
          
          # Soil_reaction_pH
          dat <- ANO.sp.ind[ANO.sp.ind$ParentGlobalID==as.character(ANO.wetland$GlobalID[i]),c('art_dekning','Soil_reaction_pH')]
          results.wetland[['original']][i,'richness'] <- nrow(dat)
          dat <- dat[!is.na(dat$Soil_reaction_pH),]
          
          if ( nrow(dat)>0 ) {
            
            val <- sum(dat[,'art_dekning'] * dat[,'Soil_reaction_pH'],na.rm=T) / sum(dat[,'art_dekning'],na.rm=T)
            # lower part of distribution
            ref <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='pH1' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'Rv']
            lim <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='pH1' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'Gv']
            maxmin <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='pH1' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'maxmin']
            # coercing x into results.wetland dataframe
            results.wetland[['scaled']][i,'pH1'] <- scal() 
            results.wetland[['non-truncated']][i,'pH1'] <- scal.2() 
            results.wetland[['original']][i,'pH1'] <- val 
            
            # upper part of distribution
            ref <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='pH2' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'Rv']
            lim <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='pH2' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'Gv']
            maxmin <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='pH2' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'maxmin']
            # coercing x into results.wetland dataframe
            results.wetland[['scaled']][i,'pH2'] <- scal() 
            results.wetland[['non-truncated']][i,'pH2'] <- scal.2() 
            results.wetland[['original']][i,'pH2'] <- val
          }
          
          
          # Nitrogen
          dat <- ANO.sp.ind[ANO.sp.ind$ParentGlobalID==as.character(ANO.wetland$GlobalID[i]),c('art_dekning','Nitrogen')]
          results.wetland[['original']][i,'richness'] <- nrow(dat)
          dat <- dat[!is.na(dat$Nitrogen),]
          
          if ( nrow(dat)>0 ) {
            
            val <- sum(dat[,'art_dekning'] * dat[,'Nitrogen'],na.rm=T) / sum(dat[,'art_dekning'],na.rm=T)
            # lower part of distribution
            ref <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='Nitrogen1' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'Rv']
            lim <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='Nitrogen1' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'Gv']
            maxmin <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='Nitrogen1' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'maxmin']
            # coercing x into results.wetland dataframe
            results.wetland[['scaled']][i,'Nitrogen1'] <- scal() 
            results.wetland[['non-truncated']][i,'Nitrogen1'] <- scal.2() 
            results.wetland[['original']][i,'Nitrogen1'] <- val 
            
            # upper part of distribution
            ref <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='Nitrogen2' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'Rv']
            lim <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='Nitrogen2' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'Gv']
            maxmin <- wetland.ref.cov.val[wetland.ref.cov.val$Ind=='Nitrogen2' & wetland.ref.cov.val$grunn==as.character(results.wetland[['original']][i,"kartleggingsenhet_1m2"]),'maxmin']
            # coercing x into results.wetland dataframe
            results.wetland[['scaled']][i,'Nitrogen2'] <- scal() 
            results.wetland[['non-truncated']][i,'Nitrogen2'] <- scal.2() 
            results.wetland[['original']][i,'Nitrogen2'] <- val
          }
          
        }
      }
      

    
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
}

# for using both sides of the plant indicators
results.wetland[['2-sided']] <- results.wetland[['non-truncated']]

# remove values >1 for 2-sided indicators
results.wetland[['2-sided']]$Light1[results.wetland[['2-sided']]$Light1>1] <- NA
results.wetland[['2-sided']]$Light2[results.wetland[['2-sided']]$Light2>1] <- NA

results.wetland[['2-sided']]$Moist1[results.wetland[['2-sided']]$Moist1>1] <- NA
results.wetland[['2-sided']]$Moist2[results.wetland[['2-sided']]$Moist2>1] <- NA

results.wetland[['2-sided']]$pH1[results.wetland[['2-sided']]$pH1>1] <- NA
results.wetland[['2-sided']]$pH2[results.wetland[['2-sided']]$pH2>1] <- NA

results.wetland[['2-sided']]$Nitrogen1[results.wetland[['2-sided']]$Nitrogen1>1] <- NA
results.wetland[['2-sided']]$Nitrogen2[results.wetland[['2-sided']]$Nitrogen2>1] <- NA

head(results.wetland[['2-sided']]) %>%
      kable("html", caption = "ANO plots for wetland ecosystems with scaled indicator values in the far-right columns.") %>% 
  kable_styling("striped") %>% scroll_box(width = "100%")

Table 6.5: ANO plots for wetland ecosystems with scaled indicator values in the far-right columns.
GlobalID	registeringsdato	klokkeslett_start	ano_flate_id	ano_punkt_id	ssb_id	program	instruks	aar	dataansvarlig_mdir	dataeier	vaer	hovedoekosystem_punkt	andel_hovedoekosystem_punkt	utilgjengelig_punkt	utilgjengelig_begrunnelse	gps	noeyaktighet	kommentar_posisjon	klokkeslett_karplanter_start	art_alle_registrert	karplanter_dekning	klokkeslett_karplanter_slutt	karplanter_feltsjikt	moser_dekning	torvmoser_dekning	stroe_dekning	jord_grus_stein_berg_dekning	stubber_kvister_dekning	alger_fjell_dekning	kommentar_ruteanalyse	fastmerker	kommentar_fastmerker	kartleggingsenhet_1m2	hovedtype_1m2	ke_beskrivelse_1m2	kartleggingsenhet_250m2	hovedtype_250m2	ke_beskrivelse_250m2	andel_kartleggingsenhet_250m2	groeftingsintensitet	bruksintensitet	beitetrykk	slatteintensitet	forekomst_ntyp	ntyp	kommentar_naturtyperegistering	side_5_note	ikke_krypende_vier_dekning	vedplanter_total_dekning	busker_dekning	tresjikt_dekning	treslag_registrert	roesslyng_dekning	roesslyngblad	pa_dekning	pa_note	pa_registrert	fa_registrert	kommentar_250m2_flate	klokkeslett_slutt	vedlegg_url	creator	creationdate	editor	editdate	hovedtype_rute	Light1	Light2	Moist1	Moist2	pH1	pH2	Nitrogen1	Nitrogen2
{36BAE83E-F412-4C90-8B56-DE6B4225BB17}	2019-08-16 11:59:59	13:19	ANO0003	ANO0003_35	20940006559500	ANO	https://nedlasting.miljodirektoratet.no/naturovervaking/ano_instruks.pdf	2019	Ole Einar Butli Hårstad	Miljødirektoratet, VAL/VMK	overskyet	vaatmark	NA	NA	NA	NA	NA	NA	NA	NA	85.1	NA	NA	20	20	60	NA	NA	NA	NA	NA	NA	V1-C-5	Åpen jordvannsmyr	svært og temmelig kalkfattige myrkanter	NA	NA	NA	NA	1	1	NA	1	nei	NA	en del gjengroing med bjørk og einer, noen store furuer	NA	1	30	12	15	NA	NA	NA	NA	NA	NA	NA	litt nakent berg	13:33	https://nin-felles-test.miljodirektoratet.no/api/Overvaking/list/36bae83e-f412-4c90-8b56-de6b4225bb17	chrpot_miljodirektoratet	2019-09-09 10:54:17	chrpot_miljodirektoratet	2019-09-09 10:54:17	V1	NA	0.9880821	NA	0.5073251	0.5130170	NA	NA	0.6919856
{662B047A-411B-42F1-8F2F-E24602CBF350}	2019-08-16 11:59:59	15:38	ANO0003	ANO0003_44	20940006559500	ANO	https://nedlasting.miljodirektoratet.no/naturovervaking/ano_instruks.pdf	2019	Ole Einar Butli Hårstad	Miljødirektoratet, VAL/VMK	sol	vaatmark	NA	NA	NA	NA	NA	NA	NA	NA	81.0	NA	NA	5	0	45	NA	NA	NA	NA	NA	NA	V1-C-1	Åpen jordvannsmyr	svært og temmelig kalkfattige myrflater	NA	NA	NA	NA	1	1	NA	1	nei	NA	NA	NA	0	7	1	3	NA	NA	NA	NA	NA	NA	NA	NA	15:45	https://nin-felles-test.miljodirektoratet.no/api/Overvaking/list/662b047a-411b-42f1-8f2f-e24602cbf350	chrpot_miljodirektoratet	2019-09-09 10:56:19	chrpot_miljodirektoratet	2019-09-09 10:56:19	V1	0.5557351	NA	NA	0.9783315	0.9992800	NA	NA	0.5718196
{F003E630-94CD-4FFF-839B-772299BF697A}	2019-08-16 11:59:59	15:07	ANO0003	ANO0003_64	20940006559500	ANO	https://nedlasting.miljodirektoratet.no/naturovervaking/ano_instruks.pdf	2019	Ole Einar Butli Hårstad	Miljødirektoratet, VAL/VMK	sol	vaatmark	NA	NA	NA	NA	NA	NA	NA	NA	79.0	NA	NA	85	75	15	NA	NA	NA	NA	NA	NA	V1-C-1	Åpen jordvannsmyr	svært og temmelig kalkfattige myrflater	NA	NA	NA	NA	1	1	NA	1	nei	NA	40% bærlyngskog T4C5, rute i overgang mellom myr og skog	NA	0	45	12	8	NA	NA	NA	NA	NA	NA	NA	NA	15:17	https://nin-felles-test.miljodirektoratet.no/api/Overvaking/list/f003e630-94cd-4fff-839b-772299bf697a	chrpot_miljodirektoratet	2019-09-09 10:57:21	chrpot_miljodirektoratet	2019-09-09 10:57:21	V1	0.5506190	NA	0.5484802	NA	NA	0.8622486	NA	0.5796588
{F8239C33-99BD-4EEE-85A9-04A10772858D}	2019-08-16 11:59:59	13:58	ANO0003	ANO0003_55	20940006559500	ANO	https://nedlasting.miljodirektoratet.no/naturovervaking/ano_instruks.pdf	2019	Ole Einar Butli Hårstad	Miljødirektoratet, VAL/VMK	sol	vaatmark	NA	NA	NA	NA	NA	NA	NA	NA	51.1	NA	NA	80	60	75	NA	NA	NA	har noen arter som er i V1C2 men bare veldig lokalt	NA	NA	V1-C-1	Åpen jordvannsmyr	svært og temmelig kalkfattige myrflater	NA	NA	NA	NA	1	1	NA	1	nei	NA	noen store furuer, litt bjørk	NA	3	17	8	12	NA	NA	NA	NA	NA	NA	NA	NA	14:12	https://nin-felles-test.miljodirektoratet.no/api/Overvaking/list/f8239c33-99bd-4eee-85a9-04a10772858d	chrpot_miljodirektoratet	2019-09-09 10:58:29	chrpot_miljodirektoratet	2019-09-09 10:58:29	V1	NA	0.9702611	NA	0.7785145	NA	0.4743074	NA	0.5484617
{95E184F4-9560-4ADC-A7A9-75AD9A97D654}	2019-08-16 11:59:59	12:06	ANO0003	ANO0003_15	20940006559500	ANO	https://nedlasting.miljodirektoratet.no/naturovervaking/ano_instruks.pdf	2019	Ole Einar Butli Hårstad	Miljødirektoratet, VAL/VMK	sol	vaatmark	NA	NA	NA	NA	NA	NA	NA	NA	71.0	NA	NA	80	80	23	NA	NA	NA	NA	NA	NA	V1-C-1	Åpen jordvannsmyr	svært og temmelig kalkfattige myrflater	NA	NA	NA	NA	1	1	NA	1	nei	NA	delvis myr, ganske mye lyng, virker som det gror igjen med bjørk, einer og furu	NA	0	53	13	8	NA	NA	NA	NA	NA	NA	NA	en bekk går gjennom flaten	12:19	https://nin-felles-test.miljodirektoratet.no/api/Overvaking/list/95e184f4-9560-4adc-a7a9-75ad9a97d654	chrpot_miljodirektoratet	2019-09-09 10:59:26	chrpot_miljodirektoratet	2019-09-09 10:59:26	V1	0.5998524	NA	0.5246274	NA	0.8311036	NA	NA	0.7270619
{6C11D423-9559-4C61-9036-C1309E0C6328}	2019-08-15 11:59:59	15:43	ANO0003	ANO0003_51	20940006559500	ANO	https://nedlasting.miljodirektoratet.no/naturovervaking/ano_instruks.pdf	2019	Ole Einar Butli Hårstad	Miljødirektoratet, VAL/VMK	nedbør	vaatmark	NA	NA	NA	NA	NA	NA	NA	NA	11.3	NA	NA	85	85	1	NA	NA	NA	NA	NA	NA	V1-C-1	Åpen jordvannsmyr	svært og temmelig kalkfattige myrflater	NA	NA	NA	NA	1	1	NA	1	nei	NA	rute ligger i en streng i fattigmyr	NA	0	5	2	0	NA	NA	NA	NA	NA	NA	NA	NA	15:50	https://nin-felles-test.miljodirektoratet.no/api/Overvaking/list/6c11d423-9559-4c61-9036-c1309e0c6328	chrpot_miljodirektoratet	2019-09-09 10:59:58	chrpot_miljodirektoratet	2019-09-09 10:59:58	V1	0.9771330	NA	NA	0.9056025	NA	0.7803061	NA	0.9422970

6.7.1.2 Scaled value analyses

In order to visualize the results we need to rearrange the results-objects from wide to long format (note that there is both a lower and an upper condition indicator for each of the functional plant indicators).

#### plotting scaled values by main ecosystem type ####
## continuing with 2-sided
res.wetland <- results.wetland[['2-sided']]

# make long version of the scaled value part
res.wetland <-
res.wetland %>%
pivot_longer(
cols = c("Light1","Light2","Moist1","Moist2","pH1","pH2","Nitrogen1","Nitrogen2"),
names_to = "fp_ind",
values_to = "scaled_value",
values_drop_na = FALSE
)

# add original values as well
res.wetland <-
res.wetland %>% add_column(original = results.wetland[['original']] %>%
pivot_longer(
cols = c("Light1","Light2","Moist1","Moist2","pH1","pH2","Nitrogen1","Nitrogen2"),
names_to = NULL,
values_to = "original",
values_drop_na = FALSE
) %>%
pull(original)
)

6.7.2 Ecosystem sub-types

And we can show the resulting scaled values as Violin plots for each indicator and main ecosystem type

Scaled vegetation indicator values for wetlands. 1 = cases where variable is lower than the reference value- 2 = cases where variablei is higher than the reference value.

Figure 6.1: Scaled vegetation indicator values for wetlands. 1 = cases where variable is lower than the reference value- 2 = cases where variablei is higher than the reference value.

This figure is a bit easier to read if we remove the sub-types with very few data points:

ggplot(res.wetland[res.wetland$hovedtype_rute %in% c("V1", "V2", "V3", "V4"),], aes(x=factor(hovedtype_rute), y=scaled_value, fill=fp_ind)) + 
  geom_hline(yintercept=0.6, linetype="dashed") + 
  geom_violin() +
#  geom_boxplot(width=0.2, color="grey") +
  geom_point(size=0.7, shape=16, color="grey") +
  facet_wrap(~factor(fp_ind,levels=c("Light1","Moist1","pH1","Nitrogen1","Light2","Moist2","pH2","Nitrogen2")), ncol = 4) + 
  xlab("Wetland ecosystem sub-type") + 
  ylab("Scaled indicator value") +
  guides(fill='none')

Scaled vegetation indicator values for the four most common sub-types of wetland. 1 = cases where variable is lower than the reference value- 2 = cases where variablei is higher than the reference value.

Figure 6.2: Scaled vegetation indicator values for the four most common sub-types of wetland. 1 = cases where variable is lower than the reference value- 2 = cases where variablei is higher than the reference value.

6.7.3 Indicator index maps

We can also show the results as a map, for instance for pH1 (the lower pH indicator), either by directly plotting the data onto the map…

#### scaled value maps ####
# keep wide format and add geometry again
res.wetland2 <- results.wetland[['2-sided']]
st_geometry(res.wetland2) <- st_geometry(ANO.wetland)

#nor <- readRDS('P:/41201785_okologisk_tilstand_2022_2023/data/rds/norway_outline.RDS')%>%
#  st_as_sf() %>%
#  st_transform(crs = crs(ANO.geo))

nor <- st_read("data/outlineOfNorway_EPSG25833.shp")%>%
  st_as_sf() %>%
  st_transform(crs = crs(ANO.geo))
#> Reading layer `outlineOfNorway_EPSG25833' from data source 
#>   `/data/Egenutvikling/41001581_egenutvikling_anders_kolstad/github/ecosystemCondition_v2/data/outlineOfNorway_EPSG25833.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -113472.7 ymin: 6448359 xmax: 1114618 ymax: 7939917
#> Projected CRS: ETRS89 / UTM zone 33N

#reg <- st_read("P:/41201785_okologisk_tilstand_2022_2023/data/regioner/regNorway_wgs84 - MERGED.shp")%>%
#  st_as_sf() %>%
#  st_transform(crs = crs(ANO.geo))

reg <- st_read("data/regions.shp")%>%
  st_as_sf() %>%
  st_transform(crs = crs(ANO.geo))
#> Reading layer `regions' from data source 
#>   `/data/Egenutvikling/41001581_egenutvikling_anders_kolstad/github/ecosystemCondition_v2/data/regions.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 5 features and 2 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -99551.21 ymin: 6426048 xmax: 1121941 ymax: 7962744
#> Projected CRS: ETRS89 / UTM zone 33N

# change region names to something R-friendly
reg$region
#> [1] "Nord-Norge"      "Midt-Norge"      "Ã\u0098stlandet"
#> [4] "Vestlandet"      "SÃ¸rlandet"
reg$region <- c("Northern Norway","Central Norway","Eastern Norway","Western Norway","Southern Norway")

regnor <- st_intersection(reg,nor)


## scaled value maps
# pH1 (lower indicator)
tm_shape(regnor) +
  tm_fill('GID_0', labels="", title="", legend.show = FALSE) + 
  tm_borders() +
  tm_shape(res.wetland2) +
  tm_dots('pH1',midpoint=NA, palette=tmaptools::get_brewer_pal("YlOrRd", 7, plot = FALSE), scale=1, legend.show = FALSE) + # 
  tm_layout(main.title = "pH (lower), wetland",legend.position = c("right", "bottom"), main.title.size=1.2) + 
  tm_add_legend(type = "fill", 
                col = c(tmaptools::get_brewer_pal("YlOrRd", 5, plot = FALSE),'grey'),
                labels = c("0.0 - 0.2", "0.2 - 0.4", "0.4 - 0.6", 
                           "0.6 - 0.8", "0.8 - 1.0", "NA"),
                title = "index values")

Figure 6.3: Functional plant indicator values for pH (lower) in wetlands.

…but here the colors and values of the data points on the map are hard to make out.

6.7.4 Regions - maps and statistics

Alternatively we can calculate and show the region-wise means and their related standard errors. But note that calculating a simple mean would be inappropriate for these data. This is because: (i) the scaled data are bound between 0 and 1, and thus follow a beta-distribution rather than a Gaussian one (ii) the ANO dataset has a nested structure

Therefore, we need to (i) use a beta-model, that (ii) can account for the nested structure of the data. Here, we apply the following function using either a glmmTMB null-model with a beta-distribution, logit link, and the nesting as a random intercept, or a simple betareg null-model with logit link if the nesting is not extensive enough for a mixed model.

expit <- function(L) exp(L) / (1+exp(L)) # since the beta-models use a logit link, we need to calculate the estimates back to the identity scale

# the function performs a glmmTMB if there's >= 5 random levels in the nesting structure
# if that is not the case, then the function performs a betareg if theres >= 2 observations
# if that is not the case either, then the function returns the value of the single observation with NA standard error

indmean.beta <- function(df) {

  st_geometry(df) <- NULL
  colnames(df) <- c("y","ran")
  
  if ( nrow(df[!is.na(df[,1]),]) >= 2 ) {
    
    if ( length(unique(df[!is.na(df[,1]),2])) >=5 ) {
      
      mod1 <- glmmTMB(y ~ 1 +(1|ran), family=beta_family(), data=df)
      
      return(c(
        expit(summary( mod1 )$coefficients$cond[1]),
        
        expit( summary( mod1 )$coefficients$cond[1] + 
                 summary( mod1 )$coefficients$cond[2] )-
          expit( summary( mod1 )$coefficients$cond[1] ),
        
        nrow(df[!is.na(df$y),]),
        summary( mod1 )$coefficients$cond[1],
        summary( mod1 )$coefficients$cond[2]
      ))
      
    } else {
      
      mod2 <- betareg(y ~ 1, data=df)
      
      return(c(
        expit(summary( mod2 )$coefficients$mean[1]),
        expit( summary( mod2 )$coefficients$mean[1] + 
                 summary( mod2 )$coefficients$mean[2] )-
          expit( summary( mod2 )$coefficients$mean[1] ),
        nrow(df[!is.na(df$y),]),
        summary( mod2 )$coefficients$mean[1],
        summary( mod2 )$coefficients$mean[2]
      ))
      
    }
    
  } else {
    
    return(c(df$y,NA,1,NA,NA))
    
  }

}


res.wetland2 = st_join(res.wetland2, regnor, left = TRUE)

regnor <- regnor %>%
  mutate(
    pH1.reg.mean = c(indmean.beta(df=res.wetland2[res.wetland2$region=="Northern Norway",c("pH1","ano_flate_id")])[1],
            indmean.beta(df=res.wetland2[res.wetland2$region=="Central Norway",c("pH1","ano_flate_id")])[1],
            indmean.beta(df=res.wetland2[res.wetland2$region=="Eastern Norway",c("pH1","ano_flate_id")])[1],
            indmean.beta(df=res.wetland2[res.wetland2$region=="Western Norway",c("pH1","ano_flate_id")])[1],
            indmean.beta(df=res.wetland2[res.wetland2$region=="Southern Norway",c("pH1","ano_flate_id")])[1]
    ),
    pH1.reg.se = c(indmean.beta(df=res.wetland2[res.wetland2$region=="Northern Norway",c("pH1","ano_flate_id")])[2],
            indmean.beta(df=res.wetland2[res.wetland2$region=="Central Norway",c("pH1","ano_flate_id")])[2],
            indmean.beta(df=res.wetland2[res.wetland2$region=="Eastern Norway",c("pH1","ano_flate_id")])[2],
            indmean.beta(df=res.wetland2[res.wetland2$region=="Western Norway",c("pH1","ano_flate_id")])[2],
            indmean.beta(df=res.wetland2[res.wetland2$region=="Southern Norway",c("pH1","ano_flate_id")])[2]
    ),
    pH1.reg.n = c(nrow(res.wetland2[res.wetland2$region=="Northern Norway" & !is.na(res.wetland2$pH1),]),
            nrow(res.wetland2[res.wetland2$region=="Central Norway" & !is.na(res.wetland2$pH1),]),
            nrow(res.wetland2[res.wetland2$region=="Eastern Norway" & !is.na(res.wetland2$pH1),]),
            nrow(res.wetland2[res.wetland2$region=="Western Norway" & !is.na(res.wetland2$pH1),]),
            nrow(res.wetland2[res.wetland2$region=="Southern Norway" & !is.na(res.wetland2$pH1),])
    )
  )

# Colour palettes
myCol <- "RdYlGn"

## scaled value maps for pH1 (lower indicator)
# mean
tm_shape(regnor) +
  tm_polygons(col="pH1.reg.mean", title="pH (lower), mean", style="fixed", 
              palette=myCol,
              breaks = seq(0,1,.2)) +
  tm_text("pH1.reg.n",col="black",bg.color="grey")

Mean index value by region for the lower pH indicator (i.e. index shows deviations towards a more acidic environment). Numbers in grey fields show the number of observations in the respective region.

Figure 6.4: Mean index value by region for the lower pH indicator (i.e. index shows deviations towards a more acidic environment). Numbers in grey fields show the number of observations in the respective region.

We can also plot the uncertainty, here represented as the spatial variation:


tm_shape(regnor) +
  tm_polygons(col="pH1.reg.se", title="pH (lower), sd", style="quantile", palette=get_brewer_pal(palette="OrRd", n=5, plot=FALSE)) +
  tm_text("pH1.reg.n",col="black",bg.color="grey")

Standard error to the mean index value by region for the lower pH indicator. Numbers in grey fields show the number of observations in the respective region.

Figure 6.5: Standard error to the mean index value by region for the lower pH indicator. Numbers in grey fields show the number of observations in the respective region.

6.7.5 Unscaled values vs. reference

We can also compare the unscaled values to the reference distribution in order to identify ecosystem types and functional plant indicators showing a deviation from the expectation. Since pH and nitrogen show some deviation in the lower indicator we exemplify this with these two indicators for unscaled values. (Only the code for the first figure is shown).

    #summary(res.wetland$kartleggingsenhet_1m2)
    #length(unique(res.wetland$kartleggingsenhet_1m2))
# 16 NiN-types to plot
    #colnames(wetland.ref.cov[['Soil_reaction_pH']])

### pH
par(mfrow=c(2,4))
## V1s
# V1-C-1
plot(density( as.matrix(wetland.ref.cov[['Soil_reaction_pH']][,c("V1-C1a","V1-C1b","V1-C1c","V1-C1d","V1-C1e")]) ,na.rm=T),
     xlim=c(1,7), type="l", main='V1-C1',xlab='pH value')
points(density(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-1",]$original,na.rm=T),
       type="l", col="red")
points(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-1",]$original,
       rep(0,length(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-1",]$original)),
       col="red")

# V1-C-2
plot(density( as.matrix(wetland.ref.cov[['Soil_reaction_pH']][,c("V1-C2a","V1-C2b","V1-C2c","V1-C2d")]) ,na.rm=T),
     xlim=c(1,7), type="l", main='V1-C2',xlab='pH value')
points(density(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-2",]$original,na.rm=T),
       type="l", col="red")
points(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-2",]$original,
       rep(0,length(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-2",]$original)),
       col="red")

# V1-C-3
plot(density( as.matrix(wetland.ref.cov[['Soil_reaction_pH']][,c("V1-C3a","V1-C3b","V1-C3c","V1-C3d")]) ,na.rm=T),
     xlim=c(1,7), type="l", main='V1-C3',xlab='pH value')
points(density(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-3",]$original,na.rm=T),
       type="l", col="red")
points(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-3",]$original,
       rep(0,length(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-3",]$original)),
       col="red")

# V1-C-4
plot(density( as.matrix(wetland.ref.cov[['Soil_reaction_pH']][,c("V1-C4a","V1-C4b","V1-C4c","V1-C4d","V1-C4e","V1-C4f","V1-C4g","V1-C4h")]) ,na.rm=T),
     xlim=c(1,7), type="l", main='V1-C4',xlab='pH value')
points(density(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-4",]$original,na.rm=T),
       type="l", col="red")
points(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-4",]$original,
       rep(0,length(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-4",]$original)),
       col="red")

# V1-C-5
plot(density( as.matrix(wetland.ref.cov[['Soil_reaction_pH']][,c("V1-C5")]) ,na.rm=T),
     xlim=c(1,7), type="l", main='V1-C5',xlab='pH value')
points(density(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-5",]$original,na.rm=T),
       type="l", col="red")
points(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-5",]$original,
       rep(0,length(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-5",]$original)),
       col="red")

# V1-C-6
plot(density( as.matrix(wetland.ref.cov[['Soil_reaction_pH']][,c("V1-C6a","V1-C6b")]) ,na.rm=T),
     xlim=c(1,7), type="l", main='V1-C6',xlab='pH value')
points(density(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-6",]$original,na.rm=T),
       type="l", col="red")
points(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-6",]$original,
       rep(0,length(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-6",]$original)),
       col="red")

# V1-C-7
plot(density( as.matrix(wetland.ref.cov[['Soil_reaction_pH']][,c("V1-C7a","V1-C7b")]) ,na.rm=T),
     xlim=c(1,7), type="l", main='V1-C7',xlab='pH value')
points(density(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-7",]$original,na.rm=T),
       type="l", col="red")
points(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-7",]$original,
       rep(0,length(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-7",]$original)),
       col="red")

# V1-C-8
plot(density( as.matrix(wetland.ref.cov[['Soil_reaction_pH']][,c("V1-C8a","V1-C8b")]) ,na.rm=T),
     xlim=c(1,7), type="l", main='V1-C8',xlab='pH value')
points(density(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-8",]$original,na.rm=T),
       type="l", col="red")
points(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-8",]$original,
       rep(0,length(res.wetland[res.wetland$fp_ind=="pH1" & res.wetland$kartleggingsenhet_1m2=="V1-C-8",]$original)),
       col="red")

Figure 6.6: Comparing the reference and variable distribution for the pH indicator in ecosystem sub-type V1.

Figure 6.7: Comparing the reference and variable distribution for the pH indicator in ecosystem sub-types other than V1.

Figure 6.8: Comparing the reference and variable distribution for the nitrogen affiliation indicator in ecosystem sub-type V1.

Figure 6.9: Comparing the reference and variable distribution for the nitrogen affiliation indicator in ecosystem sub-types other than V1.

Deviations from the reference occur mainly in V2-C1, which represents limestone-poor swamp forests. According to the functional signature from the plant community composition many of the occurrences of this ecosystem type are too acidic and may have too low availability of nitrogen.

6.7.6 Eksport file (final product)

Nothing exported yet.

5 Alien species

7 Functional Plant Indicators - Semi-Natural Ecosystems