In rivm-syso/KRWTrends: Functies voor trendrapportage KRW

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

# Laad benodigde packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, maps, sf, ggmap, 
               ggspatial)

# Laad KRWTrends pakket
# eerst installeren met devtools::install()
devtools::load_all()

Introductie

Deze handleiding beschrijft de werkwijze voor de trendanalyse algemene grondwaterkwaliteit voor de Kaderrichtlijn Water (KRW). Als input wordt gebruik gemaakt van de data uit het KRW Meetnet Grondwater (KMG), welke bestaat uit een selectie van putten uit het Landelijk- en Provinciaal Meetnet Grondwater (LMG en PMG). De werkwijze is verdeeld in zes stappen welke alle handelingen van dataselectie tot het komen tot een oordeel op grondwaterlichaam niveau beschrijven. Iedere stap bestaat uit meerdere onderdelen. Daarnaast wordt per stap aangegeven welke informatie gerapporteerd dient te worden in de achtergrondrapportage van de provincies. De volledige werkwijze en benodigde rapportage staat beschreven in het document 'Werkinstructie trendanalyse algemene grondwaterkwaliteit KRW'.

KRWTrends R-pakket

Voor de trendanalyse wordt gebruik gemaakt van het door het RIVM opgestelde R-pakket 'KRWTrends'. Dit pakket bestaat uit functies voor de analyse van trends en trendomkering in grondwater, specifiek voor het opstellen van de trendrapportages die provincies moeten opstellen voor de KRW. Daarnaast zijn er functies voor het ophalen van drempelwaardes en de grondwaterlichamen. In de stappen van de KRW trendanalyse worden de functies van KRWTrends toegelicht. Voorbeelden worden gegeven aan de hand van een dataset selectie uit het Landelijk Meetnet Grondwaterkwaliteit (LMG).

1: Samenstellen dataset voor trendanalyse KRW

Voor het uitvoeren van de trendanalyse grondwaterkwaliteit moeten de juiste gegevens uit de KMG database worden geselecteerd. Het uitgangspunt bij de dataselectie en trendanalyse is dat de kwaliteit van de data op orde is. De selectie vindt op de volgende manier plaats:

a) selecteer alle KMG putten in een grondwaterlichaam; b) selecteer alle ondiepe en diepe putfilters (volgens KMG aanduiding); c) selecteer de EU-relevante stoffen (NO3, Cd, Pb, Ni, P, Cl, As); d) selecteer alle waarnemingen over de periode 1998-2018; e) herdefinieer detectieteken "<" als 1 en geen detectieteken als 0.

In de achtergrondrapportage dient bovenstaande dataselectie te worden beschreven en eventuele aanvullende keuzes ( andere gebiedsindeling dan grondwaterlichamen, dataselectie op basis van infiltratiejaren, etc.) te worden gemotiveerd. Uiteindelijk moet bovenstaande leiden tot een data.frame, waarbij dit pakket ervan uitgaat dat de gegevens in de volgende velden staan:

veldnaam | inhoud | type ---------------|-----------------------------|-------- putfilter | identificatie put & filter | karakter xcoord | x-coordinaat, in meters volgens RD | numeriek ycoord | y-coordinaat, in meters volgens RD | numeriek meetjaar | jaar van meting | numeriek diepte | aanduiding diepte | karakter parameter | naam van stof of parameter | karakter detectielimiet | geeft aan of een waarde onder de detectielimiet ligt (1) of erboven (0) | karakter waarde | gemeten waarden | numeriek eenheid | eenheid van gemeten waarde | karakter grondwaterlichaam | grondwaterlichaam identificatiecode | karakter verwijder | markeer reeks voor verwijdering, 0 of >0 | numeriek

putfilter is de naam, of de code, van de put en filter waarop de metingen van toepassing zijn.

xcoord,ycoord zijn de coordinaten van de putten volgens het Rijksdriehoekstelsel. De coordinaten zijn gegeven in meters.

meetjaar is het jaar van de meting. Per put & filter mag er slechts 1 meting per meetjaar voorkomen.

diepte is de aanduiding van de diepte. Deze aanduiding kan zijn "diep" of "ondiep", in kleine letters. Andere aanduidingen zijn niet toegestaan.

parameter is de naam van de gemeten parameter of stof in kleine letters. Alleen de parameters cl, as, pb, no3, ni, cd en ptot zijn toegestaan.

detectielimiet geeft aan of een meting onder de rapportagegrens zit (1) of erboven (0).

waarde dit is de gemeten waarde. Alle waardes dienen positief te zijn.

eenheid is de eenheid van de gemeten parameter. Deze eenheid mag per parameter niet verschillen.

grondwaterlichaam geeft met een identificatiecode aan in welk grondwaterlichaam het putfilter zich bevindt.

verwijder is een veld om aan te geven welke metingen gemarkeerd worden. Dit veld is 0 als de stof niet gemarkeerd is, iedere waarde >0 geeft aan dat de meting is gemarkeerd.

Overzicht van data

str(grondwaterdata)
head(grondwaterdata)

2: Toets overschrijding 75% van de norm

Per grondwaterlichaam wordt er getoetst of er waarnemingen zijn die 75% van de norm overschrijden van de stoffen die geselecteerd zijn in stap 1. Hiervoor worden de normen gebruikt zoals opgenomen in het BKMW (2015). De stoffen met 75% normoverschrijdingen gaan in ieder geval door naar stap 4.

Lijst met grondwaterlichamen

Gebruik de functie getgwl() om een lijst van grondwaterlichamen te krijgen.

x <- getgwl()
print(x)

Lijst met drempelwaarden

Bepaal de drempelwaarden voor een grondwaterlichaam en stof met de funtie dw().

dw(param = "no3", gwl = "NLGW0001")

Met deze functie kan ook een kolom met drempelwaardes worden toegevoegd op basis van de parameter en het grondwaterlichaam:

grondwaterdata$norm <- mapply(dw, grondwaterdata$parameter, grondwaterdata$grondwaterlichaam)
head(grondwaterdata)

Help functies

Binnen utils.R staan enkele hulp functies voor de KRW trendanalyse. Deze functies worden aangeroepen binnen de andere functie voor het testen van de dataset en of de benodigde data aanwezig is. Bijvoorbeeld de functie testKolommen() gaat na of de verplichte kolommen aanwezig zijn. Deze hulp functies kunnen niet apart worden aangeroepen, maar worden dus binnen de andere functies gebruikt.

Detectielimieten

Gemeten concentraties beneden de detectielimiet worden in de kolom detectielimiet aangegeven met een 1. Binnen een reeks ( dat wil zeggen, binnen één filter met één parameter) worden de rapportagegrenzen welke lager dan de laagste waarneming (>RG) liggen, vervangen door 0,5 * de maximale RG onder de laagste waarneming. De rapportagegrenzen welke boven de laagste waarneming liggen houden de waarde van de rapportagegrens zelf. Dit gebeurt met de functie replaceDL op reeks niveau.

Deze functie wordt aangeroepen binnen de mktrends functie met de rpDL parameter en staat standaard aan, maar kan ook uitgezet worden. De dw.plot parameter wordt gebruikt om de drempelwaardes te plotten, en staat standaard ook aan. Zie onderstaand:

d <- grondwaterdata %>% 
    filter(parameter == "no3", diepte == "ondiep") 

# Pas RG niet aan
figuur <- mktrends("213f1", d, rpDL = FALSE, dw.plot = FALSE, make.plot = TRUE)
print(figuur)

# Pas RG wel aan
figuur <- mktrends("213f1", d, dw.plot = FALSE, make.plot = TRUE)
print(figuur)

Toets normoverschrijding

Voor het toetsten van de normoverschrijding wordt de gemeten waarde vergeleken met 75% van de drempelwaarde en met de drempelwaarde zelf. Hiervoor worden alleen de waarde >RG gebruikt. Dit wordt gedaan met de functie toetsNormoverschrijding per grondwaterlichaam. Binnen deze functie staat de norm standaard op 75% norm. Met het argument toetsnorm = 1 wordt deze op de drempelwaarde zelf gezet. De voorbeeld data is afkomstig van het grondwaterlichaam NLGW0003, toetsen van de data gaat als volgt:

# selecteer no3, vervang rapportagegrens met 0.5 * RG
d <- grondwaterdata %>% 
    filter(grondwaterlichaam == "NLGW0003", 
           parameter == "no3", diepte == "ondiep") 

# voorbeelddata is afkomstig van grondwaterlichaam NLGW0003
resultaat<-toetsNormoverschrijding(d)

knitr::kable(resultaat)

Met behulp van de digits statement kan de output van het aantal decimalen bepaald worden. Door middel van een for-loop kan er voor meerdere parameters tegelijk getoetst worden:

resultaat <- data.frame()
for (i in c("no3", "as")) {

    d <- grondwaterdata %>% 
        filter(grondwaterlichaam == "NLGW0003", parameter == i, diepte == "ondiep") 
    resultaat <- resultaat %>%
        bind_rows(data.frame(param = i, toetsNormoverschrijding(d),
                             stringsAsFactors = FALSE))
}

knitr::kable(resultaat, digits = 2)

Hiermee wordt een overzicht verkregen van stoffen waarvoor een trendtoets zal worden uitgevoerd. In de rapportage dient een tabel te worden opgenomen waarin per grondwaterlichaam is aangegeven voor welke stoffen een overschrijding is aangetroffen van 75% van de norm en van de norm, zie Tabel 1.

Tabel 1

| Grondwaterlichaam | Parameters | |-------------------+----------------------:|-------------------:| | |75% normoverschrijding | Normoverschrijding | |NLGW00XX | | |

3: Beoordeling geschiktheid dataset voor trendanalyse KRW

Voor de geselecteerde stoffen in stap 2 moet vervolgens de KMG-dataset worden beoordeeld op geschiktheid voor het uitvoeren van de KRW trendanalyse. De belangrijkste factoren hierbij zijn bemonsteringsfrequentie en rapportagegrenzen. De volgende checks dienen doorlopen te worden voor de ondiepe en diepe filters:

a) overzicht bemonsteringsfrequentie met een overzicht van het aantal putfilters met een x aantal bemonstering voor de ondiepe en diepe filters (zie document werkwijze, Tabel 2) kan de geschiktheid van de dataset voor trendanalyse KRW beoordeeld worden. Markeer de putfilters die minder dan 10x (arbitrair gekozen) zijn bemonsterd.

Onderstaand script maakt een overzicht wat kan helpen bij Tabel 2 uit de werkwijze:

# Selecteer één parameter om maar één meetronde mee te nemen, 
# voor meerdere parameters worden er extra bemonsteringen geteld

# maak overzicht van de meetfrequentie
d <- monsterFreq(grondwaterdata, param = "no3")
knitr::kable(d)

# Histogram 
d <- grondwaterdata %>% 
    filter(parameter == "no3", diepte == "ondiep") %>%
    droplevels()

a <- d %>% count(putfilter) %>%
  mutate(above = n > 10)

ggplot(a, aes(x = n, fill = above)) +
  geom_histogram(binwidth = 1,
    color = "darkblue") +
  geom_vline(aes(xintercept = mean(n), color = "gemiddelde"),
             linetype = "dashed", size = 1) +
  geom_vline(aes(xintercept = 10, color = "kritieke.grens"),
             linetype = "dashed", size = 1) +
  scale_y_continuous(breaks = 0:9) +
  ggtitle("Bemonsteringsfrequenties") + theme(plot.title = element_text(hjust = 0.5)) + 
  theme(legend.position = "right") +
  labs(x = "Aantal x bemonsterd", y = "Aantal filters") + 
  scale_color_manual(name = "lijn", values = c(gemiddelde = "blue", kritieke.grens = "red"))

b) verdeling waarnemingen 1998-2018 van de reeksen van waarnemingen per stof per putfilter dienen ten minste twee waarnemingen in de eerste planperiode KRW te zijn gelegen (2004-2009), ten minste twee waarnemingen in de tweede planperiode KRW (2010-2015) en ten minste één waarneming in de derde planperiode KRW (2016-2021). Markeer de reeks indien hier niet aan wordt voldaan.

Een overzicht per put kan verkregen worden met de meetverdeling() functie. Reeksen met onvoldoende verdeling in waarnemingen kunnen vervolgens worden gemarkeerd. Zie onderstaand script voor een voorbeeld:

d <- grondwaterdata %>%
  filter(parameter == "no3", diepte == "ondiep", putfilter == "5f1") 

resultaat <- meetverdeling(d)

knitr::kable(resultaat)

Door middel van een for-loop kan er voor meerdere putfilters een overzicht verkregen worden van de verdeling over de KRW planperiodes:

resultaat <- data.frame()
for (i in unique(grondwaterdata$putfilter)) {

  d <- grondwaterdata %>%
  filter(parameter == "no3", diepte == "ondiep", putfilter == i) 

  resultaat <- resultaat %>%
    bind_rows(data.frame(put = i, meetverdeling(d),
                         stringsAsFactors = FALSE)) 
}
# Voor een beter overzicht wordt de data van long naar wide format omgezet
resultaat <- resultaat %>% spread(key = tijdsperiode, value = aantal.metingen, fill = 0)
knitr::kable(head(resultaat))

Daarnaast kan ook een overzicht van het aantal metingen per put per meetjaar gemaakt worden met meetverdeling2():

d <- grondwaterdata

resultaat <- meetverdeling2(d, "no3", "ondiep")
knitr::kable(head(resultaat))

c) overzicht statistieken maak overzicht van de statistieken van de waarnemingen en de rapportagegrenzen in de dataset. Dit overzicht dient te worden gemaakt voor iedere stof per grondwaterlichaam, waarvoor de trendanalyse wordt uitgevoerd. Het gaat hierbij om (1) het aantal waarnemingen in de totale periode, (2) het 50e en 95e percentiel van de waarnmingen, (3) het percentage van de waarnemingen onder de rapportagegrens, (4) de laagste en hoogste rapportagegrens in de geselecteerde meetperiode en (5) de norm per stof. Zie onderstaande Tabel 3.

d <- grondwaterdata %>% 
    filter(grondwaterlichaam == "NLGW0003", parameter == "no3", diepte == "ondiep")

resultaat <- toetsStatistiek(d)

knitr::kable(resultaat, digits = 2)


# Voor meerdere stoffen
resultaat <- data.frame()
for (i in c("no3", "as")) {

    d <- grondwaterdata %>% 
        filter(grondwaterlichaam == "NLGW0003", 
               parameter == i, diepte == "ondiep")

    resultaat <- resultaat %>%
        bind_rows(data.frame(param = i, toetsStatistiek(d),
                             stringsAsFactors = FALSE))

}
knitr::kable(resultaat, digits = 2)

Tabel 3

| GWL | Parameter | Aantal waarnemingen | Waarnemingen| | \<RG (%) | Rapportage- grenzen| | Norm | |:-------------|:----------|---------------------|-------|-------------|----------|--------------------|---------|------| | | | | P50 | P95 | | Laagste | Hoogste | | | | NO3 | | | | | | | | | | Cl | | | | | | | | | | P | | | | | | | | | NLGW00XX | As | | | | | | | | | | Cd | | | | | | | | | | Ni | | | | | | | | | | Pb | | | | | | | |

Tabel 3 geeft een overzicht per grondwaterlichaam per stof m.b.t. de statistieken van de waarnemingen en rapportagegrenzen en in hoeverre die voldoende onderscheidend van elkaar zijn om de trendanalyse uit te kunnen voeren. Bij het uitvoeren van de trendanalyse in stap 4 moet er specifiek aandacht zijn voor de invloed van rapportagegrenzen. Dit wordt verder toegelicht in stap 4. De informatie in tabel 3 moet worden gebruikt bij de expert beoordeling van reeksen van waarnemingen in stap 4c.

d) Selecteer de reeksen met >= 4 waarnemingen Voor het uitvoeren van de statistische toets dienen ten minste 5 metingen in een reeks aanwezig te zijn. Dit betekent dat de trendanalyse alleen kan worden uitgevoerd als er ten minste 4 waarnemingen (<RG) en 1 rapportagegrens in een reeks zitten. Dit is een hard criterium. In deze stap worden de reeksen geselecteerd die aan dit criterium voldoen. Deze reeksen zijn goedgekeurd voor het uitvoeren van de statische toets (zie tabel 4). De overige reeksen worden niet meegenomen.

4: Statistische trendtoets

In stap 4 wordt voor de reeksen die in stap 3d zijn geselecteerd de trendtoets daadwerkelijk uitgevoerd. Dit is uitsluitend voor de grondwaterlichamen waarin een 75% normoverschrijding is aangetroffen (zie stap 2).

De trendtoets bestaat in stap 4 uit het uitvoeren van een statistische trendtoets op putfilterniveau (stap 4b), een correctie voor 'vals positieven' (stap 4c), het verzamelen van de statistieken (stap 4d) en, indien nodig, een visuele inspectie van de reeksen en een expert beoordeling (stap 4e). Voordat deze stappen worden toegelicht, wordt eerst ingegaan op de invloed van rapportagegrenzen en de invloed van gemarkeerde reeksen.

Toelichting invloed rapportagegrenzen

Binnen een reeks kunnen meerdere rapportagegrenzen met verschillende waarden aanwezig zijn die het trendresultaat beïnvloeden. In Bijlage 1 van de werkwijze zijn een aantal voorbeelden gegeven van de invloed van rapportagegrenzen op het trendresultaat, zoals kunstmatige trends en de invloed van uitbijters.

Voordat de statistische toets kan worden uitgevoerd dient de dataset te worden bewerkt. Binnen een reeks worden de rapportagegrenzen als volgt geherdefinieerd:

1) Rapportagegrenzen < alle waarnemingen (metingen >RG) binnen de reeks worden geherdefinieerd als de helft van de hoogste rapportagegrens onder de laagste waarneming; 2) Rapportagegrenzen > een of meerdere waarnemingen binnen de reeks worden geherdefinieerd als de rapportagegrens.

Hiervoor wordt gebruik gemaakt van de replaceDL functie, zie eerder.

Door de aanwezigheid van rapportagegrenzen die hoger zijn dan waarnemingen binnen een reeks kunnen kunstmatige trends worden berekend of kan het trendresultaat sterk worden beïnvloed door deze rapportagegrenzen. Daarom dienen de reeksen op filterniveau visueel te worden geïnspecteerd op de aanwezigheid van uitbijters in de rapportagegrens of op kunstmatige trends. Op basis van expert beoordeling (zie stap 4e) kunnen reeksen met kunstmatige trends worden benoemd als reeksen waarin geen trend aanwezig is. Dit zijn de trends die in tabel 4 afvallen na visuele inspectie.

Toelichting invloed van gemarkeerde reeksen

In stap 3a, 3b en 3d worden reeksen goedgekeurd en gemarkeerd op basis van het niet halen van de criteria. Het goedkeuren van reeksen wordt gedaan op basis van het minimaal benodigde aantal metingen voor een statistische toets. Dit is een hard criterium.

Het markeren van reeksen gebeurt op basis van arbitraire criteria. Met deze arbitraire criteria wordt onderscheid gemaakt in (1) reeksen waarvoor mogelijk trends kunnen worden benoemd op basis van het aantal metingen en de verdeling in de tijd en (2) reeksen waarvoor waarschijnlijk het aantal metingen niet voldoende zijn om trends te kunnen bepalen.

In tabel 4 wordt een overzicht gegeven van het aantal reeksen, het aantal goed- en afgekeurde reeksen en de benoemde trends in de goedgekeurde reeksen. Hierbij wordt onderscheid gemaakt in gemarkeerde en niet-gemarkeerde reeksen. Deze informatie dient in stap 4 te worden gebruikt bij de interpretatie van de trendresultaten en in stap 6 bij het komen tot een uitspraak op niveau van een grondwaterlichaam.

Hieronder worden de stappen 4a-e, die doorlopen moeten worden beschreven:

4.a Herdefinieer de rapportagegrenzen

Herdefinieer de rapportagegrenzen op basis van bovenstaande toelichting. Dit gebeurt op reeks niveau: per putfiler, per parameter. Dit kan gedaan worden met de replaceDL functie.

4.b Statistische analyse op putfilterniveau

De Mann-Kendall toets voor trends wordt toegepast op putfilterniveau om te testen of er een (stijgende en/of dalende) trend per filter is. Daarmee wordt een p-waarde per put en per stof berekend (de p-waarde is een getal dat aangeeft hoe waarschijnlijk het is dat een trend in een put aanwezig is. Hoe lager, hoe waarschijnlijker de trend).

De functie mktrends() berekent de Mann-Kendall statistiek per filter en per parameter. Uitbijters buiten 1,5x de interkwartielafstand kunnen binnen deze functie met rmoutlier() eruit gehaald volgens de Tukey methode. Dit wordt aangeroepen met parameter trim, welke standaard uit staat. De functie wordt als volgt gebruikt:

d <- grondwaterdata %>% 
    filter(parameter == "no3", diepte == "ondiep")

statistiek <- mktrends("5f1", d, make.plot = FALSE)
figuur <- mktrends("5f1", d, make.plot = TRUE)

print(statistiek)
print(figuur)

4.c Correctie voor 'vals positieven'

Om van een trendanalyse op putniveau te komen tot een uitspraak op grondwaterlichaamniveau moet rekening gehouden worden met de kans dat er een aantal putten kunnen zijn waar onterecht een trend is benoemd. Dit fenomeen is inherent aan het statistisch toetsen. Om voor deze zogenaamde 'vals positieven' te corrigeren wordt de Benjamini-Hochberg methode toegepast met een nominale waarde van FDR (false discovery rate) van 10%. Het resultaat van deze stap is een lijst van benoemde trends, ofwel van putfilters in een grondwaterlichaam waarvan we aannemen dat er een trend aanwezig is.

De functie bhfdr() berekent de Benjamini-Hochberg false discovery rate voor alle putfilters. Input is een vector met p-waardes van de Mann-Kendall statistiek. De laatste kolom van de output geeft aan of de p-waardes significant zijn volgens de BH-criteria. Met de volgende code wordt de functie gebruikt:

# selecteer data uit ondiepe filters voor NO3
d <- grondwaterdata %>% 
    filter(parameter == "no3", diepte == "ondiep") 

# Bereken de Mann-Kendall statistieken voor alle put locaties
i <- as.character(unique(d$putfilter))
res2  <- lapply(i, FUN = mktrends, d) %>%
    do.call("rbind", args = .) %>%
    na.omit()

# Bereken de false discovery rate voor alle locaties
fdr <- res2 %>% select(p) %>% bhfdr()

print(fdr)

4.d Verzamelen statistieken

Na het doorlopen van stappen 4a-4c kan een tabel worden samengesteld waarin per stof meerdere statistieken worden verzameld, zoals het percentage benoemde trends (stijgend of dalend), het percentage benoemde stijgende trends en het percentage benoemde dalende trends. Daarnaast wordt ook het percentage niet-benoemde trends (stijgend of dalend) gegeven. In deze tabel wordt tevens de gemiddelde helling voor de verschillende (benoemde en niet benoemde) trends gegeven.

De resultaten van de trendtoets wordt opgenomen in tabel 5 (zie 'werkinstructie'). Na het uitvoeren van de trendtoetsen vindt een visuele inspectie op alle grafieken plaats (4e) als gevolg van het vervangen van de iteraties met rapportagegrenzen.

De trendanalyse wordt gedaan met de functie trendAnalyse. Hierin zitten zowel stap 4a (aanpassen RG met replaceDL), stap 4b (Mann-Kendall trendanalyse met mktrends) als 4c (correctie vals-positieven met bhfd), om vervolgens de statistieken te verzamelen en benoemen. Om de trendanalyse voor NO3 uit te voeren gebruiken we de volgende code:

Trendanalyse tabel

d <- grondwaterdata %>% 
    filter(parameter == "no3", diepte == "ondiep") 

# Totaal aantal reeksen waarvoor trendanalyse wordt gedaan, nodig om 2 percentages te berekenen
aantal.reeks <- paste(d$putfilter, d$parameter, sep = " - ") %>% unique %>% length()

resultaat <- trendAnalyse(d, n.reeks = aantal.reeks)
knitr::kable(resultaat, digits = 2)

Om voor een aantal parameters of stoffen een tabel te maken kan een for-loop gebruikt worden:

resultaat <- data.frame()
for (i in c("no3", "as")) {

    d <- grondwaterdata %>% 
        filter(parameter == i, diepte == "ondiep") 

    aantal.reeks <- paste(d$putfilter, d$parameter, sep = " - ") %>% unique %>% length()

    resultaat <- resultaat %>%
        bind_rows(data.frame(param = i, trendAnalyse(d, n.reeks = aantal.reeks),
                             stringsAsFactors = FALSE))

}
knitr::kable(resultaat, digits = 2)

4.e Visuele inspectie van reeksen

Invloed van rapportagegrenzen Na het uitvoeren van de trendtoetsen worden de reeksen op putfilterniveau visueel geïnspecteerd. Hierdoor wordt duidelijk:

1) Welke metingen in de reeks rapportagegrenzen zijn; 2) In hoeverre een rapportagegrens een uitbijter is (zie hieronder); 3) In hoeverre een benoemde trend kunstmatig is (zie hieronder).

Uitbijters in rapportagegrenzen zijn hier gedefinieerd als rapportagegrenzen die afwijkend zijn van de andere rapportagegrenzen in de reeks en die het resultaat van de trendanalyse beïnvloeden. Met 'kunstmatige trends' bedoelen we hier trends die benoemd worden maar gebaseerd zijn op verschillen in de hoogte van de rapportagegrenzen in de dataset. De gevonden trend is dan eigenlijk een 'trend in rapportagegrenzen'.

Op basis van expert beoordeling kan besloten worden om uitbijters in rapportagegrenzen te verwijderen. In het geval van een kunstmatige trend kan besloten worden om de benoemde trend te verwijderen. Deze reeksen worden dan beschouwd als reeksen waarin geen trend aanwezig is. In 'Bijlage 1 van de Werkinstructie trendanalyse algemene grondwaterkwaliteit KRW' zijn een aantal voorbeelden gegeven van expert beoordeling van grafieken m.b.t. de invloed van rapportagegrenzen.

Voorbeeld

Middels onderstaand script kan voorbeeld 2 gevisualiseerd worden om de invloed van uitbijters in rapportagegrenzen te bekijken. De afwezigheid van een trendlijn laat zien dat de trend niet significant is (pas dan wordt deze geplot).

# Rapportagegrens aanpassen volgens methode
d <- grondwaterdata %>% 
    filter(parameter == "as", diepte == "ondiep") 

statistiek <- mktrends("17f1", d, rpDL = TRUE, make.plot = FALSE)
figuur <- mktrends("17f1", d, , rpDL = TRUE, dw.plot = FALSE, make.plot = TRUE)

print(statistiek)
print(figuur)

# Rapportagegrens gelijk aan 0.5 RG
d <- grondwaterdata %>% 
    filter(parameter == "as", diepte == "ondiep") %>%
  mutate(waarde = ifelse(detectielimiet == 1, 0.5 * waarde, waarde))

statistiek0.5 <- mktrends("17f1", d, rpDL = FALSE, make.plot = FALSE)
figuur0.5 <- mktrends("17f1", d, rpDL = FALSE, dw.plot = FALSE, make.plot = TRUE)

print(statistiek0.5)
print(figuur0.5)

# Rapportagegrens gelijk aan 0
d <- grondwaterdata %>% 
    filter(parameter == "as", diepte == "ondiep") %>%
  mutate(waarde = ifelse(detectielimiet == 1, 0, waarde))

statistiek0 <- mktrends("17f1", d, rpDL = FALSE, make.plot = FALSE)
figuur0 <- mktrends("17f1", d, rpDL = FALSE, dw.plot = FALSE, make.plot = TRUE)

print(statistiek0)
print(figuur0)

Uiteindelijk dienen de verzamelde statistieken voor de betreffende trendtoetsen verwerkt te worden in Tabel 5 van de 'Werkinstructie trendanalyse algemene grondwaterkwaliteit KRW'

4.f Verzamelen van kaartbeelden

Aangezien grondwaterlichamen een groter gebied bestrijken die niet homogeen zijn m.b.t. bodemtype, hydrologie en landgebruik, kan het lastig zijn om een eenduidige uitspraak te doen over stijgende/dalende trends in het hele grondwaterlichaam. Daarom is het zinvol te kijken naar de ruimtelijke spreiding van de benoemde trends. Hiervoor wordt een kaart samengesteld, waarbij de stijgende trends worden aangegeven met een rode stip, de dalende trends met een groene stip en de overige punten met een blauwe stip.

Met de basemap() functie kan een kaart worden gemaakt van de put locaties met onderliggend de grondwaterlichamen en provinciegrenzen. Dit wordt met het volgende stukje code gedaan:

# Data selecteren uit grondwaterlichaam
d <- grondwaterdata %>% 
    distinct(putfilter, .keep_all = TRUE)


pts <- st_as_sf(d, coords = c("xcoord", "ycoord"), crs = 28992) %>% select(waarde)

p <- basemap()
p <- p + geom_sf(data = pts) 
p

5: Statistische toets voor trendomkering

In deze stap worden statistische toetsen uitgevoerd voor trendomkering en de bijbehorende statistiek verzamelt voor de stoffen die geselecteerd zijn in stap 2. Voor het detecteren van trendomkering wordt dezelfde procedure gevolgd als voor het detecteren van trends (stap 4), met enkele uitzonderingen welke in de Werkwijze zijn toegelicht. Het statistisch toetsen van trendomkering is erg complex. Na het toepassen van de toetsen en de correctie voor vals-positieven, moeten de reeksen waarin een trendomkering wordt benoemd daarom altijd visueel geïnspecteerd worden.

5.a Trendomkering

Met de functie trendReversal() wordt gekeken of er binnen een reeks twee achtervolgende trends in verschillende richtingen zijn. Hiervoor wordt de robuuste Theil-Sen helling gebruikt. Vervolgens wordt het keerpunt bepaalt aan de hand van de kleinste discrepantie tussen de rechte lijnen.

Met de volgende code kan de trendomkering bepaalt worden:

d <- grondwaterdata %>% 
    filter(parameter == "no3", diepte == "ondiep") 

statistiek <- trendReversal("213f1", d, trim = FALSE, make.plot = FALSE)

knitr::kable(statistiek)

5.b Correctie voor 'vals-positieven'

De correctie voor 'vals-positieven' vindt wederom plaats volgens de Benjamini-Hochberg methode met de functie bhfdr(), zie stap 4b.

5.c Visuele inspectie trendomkering

Middels onderstaande code kan de trendomkering gevisualiseerd worden:

d <- grondwaterdata %>% 
    filter(parameter == "no3", diepte == "ondiep") 

statistiek <- trendReversal("213f1", d, trim = FALSE, make.plot = FALSE)
figuur <- trendReversal("213f1", d, dw.plot = FALSE, trim = FALSE, make.plot = TRUE)

print(statistiek)
print(figuur)

5.d Verzamelen van statistieken

De functie trendomkeringAnalyse() voert de trendomkeringanalyse uit met de functie trendReversal(). Vervolgens worden de trends gecorrigeerd voor 'vals-positieven' met de bhfdr() functie om tot slot de bijbehorende statistieken te genereren.

Voor het berekenen van de statistieken voor trendomkering voor NO3 gebruiken we de volgende code:

d <- grondwaterdata %>% 
    filter(parameter == "no3", diepte == "ondiep")

aantal.reeks <- paste(d$putfilter, d$parameter, sep = " - ") %>% unique %>% length()

resultaat <- trendomkeringAnalyse(d, n.reeks = aantal.reeks)

knitr::kable(resultaat, digits = 2)

Om voor een aantal parameters of stoffen een tabel te maken kan net als bij de trendanalyse een for-loop gebruikt worden:

resultaat <- data.frame()
for (i in c("no3", "as")) {

    d <- grondwaterdata %>% 
        filter(parameter == i, diepte == "ondiep") 
    aantal.reeks <- paste(d$putfilter, d$parameter, sep = " - ") %>% unique %>% length()

    resultaat <- resultaat %>%
        bind_rows(data.frame(param = i, trendomkeringAnalyse(d, n.reeks = aantal.reeks),
                             stringsAsFactors = FALSE))
}

knitr::kable(resultaat, digits = 2)

De resultaten van stap 5 geeft een overzicht van de statistieken betreffende trendomkering. Deze gegevens dienen in de rapportage te worden ingevuld in Tabel 6 uit de 'Werkinstructie'.

6: Oordeel op niveau van grondwaterlichaam

Het verzamelen van de informatie uit bovenstaande stappen geeft een overzicht van de statistieken die verzameld kunnen worden van een verzameling meetreeksen binnen een grondwaterlichaam.

Door het combineren van de verzamelde statistieken uit stap 4 met de kaarten waarop de geografische spreiding van de benoemde trends inzichtelijk is gemaakt, kan de afweging worden gemaakt of er wel of geen sprake is van een trend in een grondwaterlichaam.

rivm-syso/KRWTrends documentation built on Nov. 20, 2021, 9:52 a.m.