Opbouw van UFP rond Schiphol

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code. Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ##

Code beschikbaar gesteld onder voorwaarden CC BY-NC-SA 4.0

[Attribution-NonCommercial-ShareAlike 4.0 International]

RIVM, 23 Februari2022

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

In dit document worden eerder los bepaalde UFP achtergronden, verkeersbijdragen en bijdragen van luchtvaart (grond en in de lucht) in de omgeving van Schiphol in 2017 en 2018 gefit aan gemeten totale UFP deeltjesaantallen in 36 windsectoren op 10 meetlocaties. De meetlocatie 6 (NH1, begin Polderbaan) wordt weggelaten, zie de eerdere RIVM rapportage voor de redenen daarvan.

Om op het laatst het verschil tussen de meetwaarden en de bijdragen vanuit de achtergrond en verkeer toe te wijzen aan de bijdragen van luchtvaart is het noodzakelijk dat de verschillende bijdragen zo onafhankelijk mogelijk worden bepaald. Daarom wordt de omrekening van NOx verkeersbijdragen naar UFP apart gecheckt.

De achtergronden worden in een aparte Excel file toegelicht. Bij de bepaling hiervan is deels”expert-judgement” gebruikt om te schatten vanuit welke windrichtingen (van buiten het studiegebied) er redelijkerwijs sprake is van achtergrond. Wegen buiten het studiegebied zijn namelijk wel voor een deel expliciet doorgerekend.

———————————————————————————

Stap 1, Invoer Datavelden:

Naam: naam meetlocatie, zie eerdere rapport.

windr: windrichting van de meting.

UFP: gemiddelde meetwaarde van ufp in aantal deeltjes per kubieke centimeter, zie eerder rapport.

NOx: berekende NOx bijdrage in microgram per kubieke meter.

Background: gemiddelde meetwaarde ufp in de aangenomen achtergrond in aantal deeltjes per kubieke centimeter.

BGR_P1: gemiddelde meetwaarde ufp in de aangenomen achtergrond in aantal deeltjes per kubieke centimeter gedurende enkel meetperode 1.

BGR_P2: gemiddelde meetwaarde ufp in de aangenomen achtergrond in aantal deeltjes per kubieke centimeter gedurende enkel meetperode 2.

Mazaheri: berekende ufp ten gevolge van luchtvaartactiviteiten in aantal deeltjes per kubieke centimeter, zie eerder rapport.

Uren: Aantal uren dat er bij deze combinatie van station en windrichting data beschikbaar is.

Periode: Meetperiode 1 of 2.

Eerst de data inlezen …

###############################################################################
library(simpleboot)

# setwd("HIER_PAD_INVOEGEN") 
UFPU = read.table("AlleData_18Feb22.txt", header=TRUE, sep = '\t')
head (UFPU)

———————————————————————————

Stap 2, Check UFP van verkeer

We verifiëren de omrekening van NOx verkeersbijdragen naar geschatte UFP bijdragen. Doe dit door een fit van de berekende NOx bijdragen van wegverkeer en de achtergronden aan de meetdata voor de locaties dat er naar verwachting wel wegverkeerbijdragen en weinig/geen luchtvaartbijdragen zijn. We gaan er van uit dat er weinig luchtvaartbijdragen en veel verkeersbijdragen zijn als in de data Mazaheri < 50 & Nox > 2. Deze waarden zijn een beetje arbitrair. Andere keuzes voor de grenzen geven iets andere resultaten. Voor de overall resultaten maakt dat niet wezenlijk uit. Let wel: deze fit is alleen om te verifiëren dat verkeersbijdragen goed los van de luchtvaartbijdragen kunnen worden bepaald.

##############################################################################
# ------- Fit verkeer met lm() met weging van het aantal uren!
##############################################################################
# Selecteer de data: veel verkeersbijdragen en weing van de luchtvaart.
ufp_uren <- subset(UFPU, select=c("naam", "UFP", "NOx", "Background", "Mazaheri", "Uren"),   
        NOx > 2 & Mazaheri < 50 )

# Definieer als weegfactoren de aantallen uren. 
wgt <- ufp_uren$Uren
ufps.fit <- lm( UFP ~ 0 + Background + NOx , data = ufp_uren, weights = wgt )
print(summary(ufps.fit))

Call:
lm(formula = UFP ~ 0 + Background + NOx, data = ufp_uren, weights = wgt)

Weighted Residuals:
   Min     1Q Median     3Q    Max 
-25539 -13655  -7287   3483  33838 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
Background   0.69242    0.05088   13.61 7.25e-14 ***
NOx        662.40101   67.24641    9.85 1.34e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15880 on 28 degrees of freedom
Multiple R-squared:  0.9784,    Adjusted R-squared:  0.9768 
F-statistic: 632.7 on 2 and 28 DF,  p-value: < 2.2e-16

predicted <- predict(ufps.fit,newdata=ufp_uren)
plot(ufp_uren$UFP, predicted, main = "Verkeer, 7 km, met gewogen fit",  
     cex = 0.95, xlim=c(0,30000), ylim=c(0,30000)) ; abline(0,1)

———————————————————————————

Stap 3, Gecombineerde fit

Vervolgens worden de afzonderlijke bijdragen (background, NOx en Mazaheri) integraal aan de gemeten UFP concentraties gefit. De achtergronden zijn apart bepaald, we schalen die alleen omdat we nu alles in een analyse integreren, de vorm en structuur blijft ongewijzigd. De bijdragen van wegverkeer worden ook alleen geschaald omdat we van hierboven weten dat de vorm van de verdeling in principe past bij de metingen. Voor de bijdragen van de luchtvaart verifiëren we of de verdeling past bij de metingen. De onzekerheden in de fit worden via een bootstrap procedure geschat.

##############################################################################
# ------- Fit alles gecombineerd met lm() met weging van het aantal uren.
# ------- Doe ook bootstrap voor onzekerheden.
##############################################################################
# Selecteer alle data.
ufp_uren <- subset(UFPU, select=c("naam", "UFP", "NOx", "Background", "Mazaheri", "Uren") )
# Definieer als weegfactoren de aantallen uren. 
wgt <- ufp_uren$Uren
ufps.fit <- lm( UFP ~ 0 + Background + NOx + Mazaheri , data = ufp_uren, weights = wgt )
print(summary(ufps.fit))

Call:
lm(formula = UFP ~ 0 + Background + NOx + Mazaheri, data = ufp_uren, 
    weights = wgt)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-220862  -22486   -3102   24868  318726 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
Background   1.04114    0.04123   25.25   <2e-16 ***
NOx        675.44635   52.15934   12.95   <2e-16 ***
Mazaheri     0.78259    0.02825   27.71   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 53960 on 393 degrees of freedom
Multiple R-squared:  0.9284,    Adjusted R-squared:  0.9279 
F-statistic:  1699 on 3 and 393 DF,  p-value: < 2.2e-16

predicted <- predict(ufps.fit,newdata=ufp_uren)
plot(ufp_uren$UFP, predicted, main = "Schiphol met gewogen fit, alle data",  cex = 0.65, 
     xlim=c(0,70000), ylim=c(0,70000)) ; abline(0,1)

##############################################################################
# Een simpele bootstrap ... 
# https://www.rdocumentation.org/packages/simpleboot/versions/1.1-7/topics/lm.boot 
lboot <- lm.boot(ufps.fit, rows = TRUE, R=10000)
summary(lboot)

BOOTSTRAP OF LINEAR MODEL  (method = rows)

Original Model Fit
------------------
Call:
lm(formula = UFP ~ 0 + Background + NOx + Mazaheri, data = ufp_uren, 
    weights = wgt)

Coefficients:
Background         NOx    Mazaheri  
    1.0411    675.4464      0.7826  

Bootstrap SD's:
 Background          NOx     Mazaheri  
 0.04206733  74.99985779   0.07070741  

———————————————————————————

Periode 1

Doe nu hetzelfde voor de data van periode 1:

• Achtergrond en luchtvaartbjdragen specifiek voor deze periode;

• Verkeer generiek voor 2018 omdat we die niet per periode kunnen bepalen.

##############################################################################
# ------- Fit PERIODE 1 gecombineerd met lm() met weging van het aantal uren.
# ------- Doe ook bootstrap voor onzekerheden.
##############################################################################
# Selecteer alle data.
ufp_uren <- subset(UFPU, select=c("naam", "UFP", "NOx", "Background", "Mazaheri", "Uren", "BGR_P1"), 
                   Periode == 1 )
# Definieer als weegfactoren de aantallen uren. 
wgt <- ufp_uren$Uren
ufps.fit <- lm( UFP ~ 0 + BGR_P1 + NOx + Mazaheri , data = ufp_uren, weights = wgt )
print(summary(ufps.fit))

Call:
lm(formula = UFP ~ 0 + BGR_P1 + NOx + Mazaheri, data = ufp_uren, 
    weights = wgt)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-298834  -19477   -2192   15721  263994 

Coefficients:
          Estimate Std. Error t value Pr(>|t|)    
BGR_P1   8.830e-01  6.538e-02  13.505  < 2e-16 ***
NOx      1.018e+03  1.239e+02   8.221 4.18e-14 ***
Mazaheri 9.509e-01  4.858e-02  19.574  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 51430 on 177 degrees of freedom
Multiple R-squared:  0.9237,    Adjusted R-squared:  0.9224 
F-statistic: 714.5 on 3 and 177 DF,  p-value: < 2.2e-16

predicted <- predict(ufps.fit,newdata=ufp_uren)
plot(ufp_uren$UFP, predicted, main = "Schiphol met gewogen fit, Periode 1",  cex = 0.65, 
     xlim=c(0,70000), ylim=c(0,70000)) ; abline(0,1)

lboot <- lm.boot(ufps.fit, rows = TRUE, R=10000)
summary(lboot)

BOOTSTRAP OF LINEAR MODEL  (method = rows)

Original Model Fit
------------------
Call:
lm(formula = UFP ~ 0 + BGR_P1 + NOx + Mazaheri, data = ufp_uren, 
    weights = wgt)

Coefficients:
   BGR_P1        NOx   Mazaheri  
   0.8830  1018.1519     0.9509  

Bootstrap SD's:
      BGR_P1           NOx      Mazaheri  
  0.08425056  267.55642282    0.18026863  

———————————————————————————

Periode 2

Doe nu hetzelfde voor de data van periode 2:

• Achtergrond en luchtvaartbjdragen specifiek voor deze periode;

• Vier windrichtingen kwamen in periode 2 in het meetgebied niet voor en er waren er geen gemeten data, hiervoor zijn de waarden van periode 1 gebruikt.

• Verkeer generiek voor 2018 omdat we die niet per periode kunnen bepalen.

##############################################################################
# ------- Fit PERIODE 2 gecombineerd met lm() met weging van het aantal uren.
# ------- Doe ook bootstrap voor onzekerheden.
##############################################################################
# Selecteer alle data.
ufp_uren <- subset(UFPU, select=c("naam", "UFP", "NOx", "Background", "Mazaheri", "Uren", "BGR_P2"), 
                   Periode == 2 )
# Definieer als weegfactoren de aantallen uren. 
wgt <- ufp_uren$Uren
ufps.fit <- lm( UFP ~ 0 + BGR_P2 + NOx + Mazaheri , data = ufp_uren, weights = wgt )
print(summary(ufps.fit))

Call:
lm(formula = UFP ~ 0 + BGR_P2 + NOx + Mazaheri, data = ufp_uren, 
    weights = wgt)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-139314  -25007   -3247   35956  173263 

Coefficients:
          Estimate Std. Error t value Pr(>|t|)    
BGR_P2     1.05561    0.04984   21.18   <2e-16 ***
NOx      619.67167   56.55073   10.96   <2e-16 ***
Mazaheri   0.68633    0.03296   20.83   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 51860 on 213 degrees of freedom
Multiple R-squared:  0.9418,    Adjusted R-squared:  0.941 
F-statistic:  1149 on 3 and 213 DF,  p-value: < 2.2e-16

predicted <- predict(ufps.fit,newdata=ufp_uren)
plot(ufp_uren$UFP, predicted, main = "Schiphol met gewogen fit, Periode 2",  cex = 0.65, 
     xlim=c(0,70000), ylim=c(0,70000)) ; abline(0,1)

lboot <- lm.boot(ufps.fit, rows = TRUE, R=10000)
summary(lboot)

BOOTSTRAP OF LINEAR MODEL  (method = rows)

Original Model Fit
------------------
Call:
lm(formula = UFP ~ 0 + BGR_P2 + NOx + Mazaheri, data = ufp_uren, 
    weights = wgt)

Coefficients:
  BGR_P2       NOx  Mazaheri  
  1.0556  619.6717    0.6863  

Bootstrap SD's:
     BGR_P2          NOx     Mazaheri  
 0.04246216  69.81639278   0.05656197  

———————————————————————————

Stap 4, Bijdragen luchtvaart

Nu gaan we de geschatte gemiddelde bijdragen van de luchtvaart vergelijken met de berekende gemiddelde bijdragen. We doen dit door de gefitte achtergrond en berekende verkeersbijdragen van de metingen af te trekken en het verschil dan te vergelijken met de gefitte berekende bijdragen van de luchtvaart. Deze laatste twee moeten goed bij elkaar passen, dat testen we met de correlatie tussen de twee sets. Net als in de eerdere rapportage nemen we de gemiddelden per meetlocatie.

ufp_uren <- subset(UFPU, select=c("naam", "UFP", "NOx", "Background", "Mazaheri", "Uren") )
wgt <- ufp_uren$Uren

SNAMES <- ufp_uren$naam
SUNAMES <- unique(SNAMES)
print(SUNAMES)

 [1] "Oude_Meer"    "Ookmeer"      "Amstelveen"   "Aalsmeer"     "Vijfhuizen"   "Nieuwe_Meer"  "Hoofddorp"    "Spaarnewoude" "Badhoevedorp"
[10] "Oude_Meer2"   "Ookmeer2"    

# uitgangsput is de fit aan alle data:
ufps.fit <- lm( UFP ~ 0 + Background + NOx + Mazaheri , data = ufp_uren, weights = wgt )
cf = coefficients(ufps.fit)
print(cf)

 Background         NOx    Mazaheri 
  1.0411405 675.4463519   0.7825934 

# Aggregeer de waarden per meetpunt in een nieuwe set tabellen.
Aircraft.naam <- array(1:15, dim=c(1))
Aircraft.estimate <- array(1:15, dim=c(1))
Aircraft.Mazaheri <- array(1:15, dim=c(1))
Aircraft.estimateUren <- array(1:15, dim=c(1))
Aircraft.MazaheriUren <- array(1:15, dim=c(1))

# Bepaal de gemiddelde waarden per meetlocatie.
for (i in 1:length(SUNAMES))  {
  single = ufp_uren[ufp_uren$naam == SUNAMES[i]  , ]
  single$mod = single$UFP - cf[1] * single$Background - cf[2] * single$NOx
  Aircraft.naam[i] = single$naam[1]
  
  # Gewogen met aantal uren:
  Aircraft.estimateUren[i] = mean(single$mod * single$Uren) / mean(single$Uren)
  Aircraft.MazaheriUren[i] = cf[3] * mean(single$Mazaheri * single$Uren)  / mean(single$Uren)
}

# Maak plots:
plot(Aircraft.estimateUren, Aircraft.MazaheriUren, main = "Gemiddelde luchtvaartbijdragen per locatie, gewogen met uren",  cex = 0.95, xlim=c(0,15000), ylim=c(0,15000)) ; abline(0,1)

# En de correlatie tussen de geaggregeerde geschatte luchtvaartbijdragen en Mazaheri:
print("Correlatie:")

[1] "Correlatie:"

print(cor(Aircraft.estimateUren, Aircraft.MazaheriUren))

[1] 0.8774251

———————————————————————————

Einde

———————————————————————————

@ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @ NOTE:

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I. When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed. @ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @