Difference between revisions of "Kopra"

From Testiwiki
Jump to: navigation, search
(first draft)
 
(started to document Kopra project's R code)
Line 1: Line 1:
 
[[Category:Classical air pollutants]]
 
[[Category:Classical air pollutants]]
−
{{encyclopedia}}
+
{{study|moderator=Jouni}}
−
'''Kopra''' is a research project about health impacts of fine particles in Finland.
+
'''Kopra''' was a research project about health impacts of fine particles in Finland.
  
 
* [http://www.ymparisto.fi/default.asp?contentid=69477&lan=en Project homepage]
 
* [http://www.ymparisto.fi/default.asp?contentid=69477&lan=en Project homepage]
 +
 +
== Research question ==
 +
 +
== Answer ==
 +
 +
== Rationale ==
 +
 +
=== Data ===
 +
 +
{{defend|# |See N:\YMAL\Projects\R83_piltti filelist_kopra.txt for data list.|--[[User:Jouni|Jouni]] ([[User talk:Jouni|talk]]) 14:12, 18 April 2017 (UTC)}}
 +
 +
=== Calculations ===
 +
 +
==== Marko's iF calculations ====
 +
 +
:''The codes are from U:\arkisto_kuopio\huippuyksikko\Tutkimus\R19_kopra\Mallit\IntakeFraction2. This section is for documentation only. The codes do not run.
 +
 +
'''Tainio_iF_emis_eur_pop_all_01.R
 +
 +
<rcode>
 +
#7.3.2008, Marko Tainio
 +
 +
# This model calculates intake fractions (iF) for European
 +
# primary PM2.5 emissions. Description of model:
 +
# - Emission source: European primary PM2.5 emissions from different countries
 +
# - Dispersion: Dispersion over the Europe
 +
# - Population: European population divided to different countries (38)
 +
# (Note: the population data does not include Finnish population!!)
 +
 +
#Work folders - the location of data and where the output is recorded
 +
wrk_input = "N:/Huippuyksikko/Tutkimus/R19_Kopra/Data 2/"
 +
wrk_output = "N:/Huippuyksikko/Tutkimus/R19_Kopra/Mallit/IntakeFraction2/Result/"
 +
 +
#Breathing rate, 20 m^3/day -> unit g/s
 +
BR = 20/(24 * 60 * 60)
 +
 +
# Calculation for different countries. The iF's are calculated separately
 +
# for all emission sources (countires). There are 48 sources in data.
 +
for (k in 1:48) {
 +
 +
# This reads from file the titles of countries to be used in calculations
 +
maa = read.table("N:/Huippuyksikko/Tutkimus/R19_Kopra/Mallit/IntakeFraction2/EUR_country.txt", sep="")
 +
maa = maa[,k]
 +
 +
#Emission volume - this part estimates emission volume in unit g/s
 +
Q = read.csv("N:/Huippuyksikko/Tutkimus/R19_Kopra/Mallit/IntakeFraction2/Eur_emission_05.csv", header=T)
 +
 +
#Concentration data
 +
#The order in list debends on the prepared emission data
 +
rm(conc)
 +
conc = "N:/Huippuyksikko/Tutkimus/R19_Kopra/Data 2/Concentration_primary_Europe_2000/EMEP_monthly_PM2_5_2000/"
 +
conc = paste(conc,"EMEP_2000_PM_2_5_hirlam_v3_8_1_",maa,"_200012.txt", sep="")
 +
conc = read.table(conc, header=T, sep="")
 +
 +
#Population data:
 +
rm(pop)
 +
pop = paste(wrk_input, "Population_Combined/Pop_eur00_all3.txt", sep="")
 +
pop = read.table(pop, header=T, sep="")
 +
 +
#Co-ordination check-in - this check that co-ordinate system is same
 +
#in population and concentration data.
 +
check = (conc["lon"] - pop["Lo"]) + (conc["lat"] - pop["La"])
 +
sum(check)
 +
 +
#This part slice the concentration data from the data file
 +
#(that is, LO and LA data is sliced oof)
 +
month = c("jan", "feb", "mar", "apr", "may", "jun", "jul", "aug", "sep", "oct", "nov", "dec", "annual")
 +
conc = data.frame(conc[month])
 +
 +
#Unit change in concentration data: kg/m^3 -> g/m^3
 +
conc = conc * 1000
 +
 +
#This part slice the population data from the data file
 +
#(that is, LO and LA data is sliced off)
 +
rm(country_eur)
 +
country_eur = c("Sum_Austri", "Sum_Belgiu", "Sum_Cyprus", "Sum_Czech_", "Sum_Denmar", "Sum_Estoni", "Sum_France", "Sum_German", "Sum_Greece", "Sum_Hungar", "Sum_Irish_", "Sum_Italy", "Sum_Latvia", "Sum_Lithua", "Sum_Luxemb", "Sum_Malta", "Sum_Nether", "Sum_Poland", "Sum_Portug", "Sum_Romani", "Sum_Slovak", "Sum_Sloven", "Sum_Spain", "Sum_UK", "Sum_Sweden", "Sum_Bulgar", "Sum_Norway", "Sum_Belaru", "Sum_Russia", "Sum_Switze", "Sum_Ukrain", "Sum_Moldov", "Sum_Croati", "Sum_Serbia", "Sum_Albani", "Sum_Macedo", "Sum_Turkey", "Sum_bosni", "Sum_Finl", "Sum_all")
 +
pop = data.frame(pop[country_eur])
 +
 +
#X must be defined before for -loop
 +
x = NULL
 +
 +
#Month/concentration - month 13 in annual mean
 +
for (i in 1:13) {
 +
 +
#Removal of huge values
 +
#There are some weird very high values in some concentration datasets.
 +
#This code replace those values with value 0.
 +
y = conc[,i]
 +
y = ifelse(y > 100, 0, y)
 +
 +
#Country/population - there are 40 different population data files
 +
for (j in 1:40) {
 +
 +
x[j] = sum(y * pop[,j])
 +
 +
}
 +
 +
#Calculation of iF
 +
x = (x * BR) / Q[k,i]
 +
 +
#Names for the data
 +
names(x) = c("Sum_Austri", "Sum_Belgiu", "Sum_Cyprus", "Sum_Czech_", "Sum_Denmar", "Sum_Estoni", "Sum_France", "Sum_German", "Sum_Greece", "Sum_Hungar", "Sum_Irish_", "Sum_Italy", "Sum_Latvia", "Sum_Lithua", "Sum_Luxemb", "Sum_Malta", "Sum_Nether", "Sum_Poland", "Sum_Portug", "Sum_Romani", "Sum_Slovak", "Sum_Sloven", "Sum_Spain", "Sum_UK", "Sum_Sweden", "Sum_Bulgar", "Sum_Norway", "Sum_Belaru", "Sum_Russia", "Sum_Switze", "Sum_Ukrain", "Sum_Moldov", "Sum_Croati", "Sum_Serbia", "Sum_Albani", "Sum_Macedo", "Sum_Turkey", "Sum_bosni", "Sum_Finl", "Sum_all")
 +
 +
#Recording of results
 +
write.table(x, paste(wrk_output, file = "European_emission/2000/if_emis_eur_year_2000_country_",maa,"_month_",i,"_pop_all.txt", sep=""), row.names = TRUE, col.names = FALSE)
 +
 +
}
 +
 +
}
 +
 +
#Remove of all data from memory
 +
rm(wrk_input, wrk_output, BR, k, maa, Q, pop, check, conc, month, country_eur, x, i, y, j)
 +
</rcode>
 +
 +
'''Tainio_iF_emis_eur_pop_all_data-conbination.R
 +
 +
<rcode>
 +
# 18.4.2008 Marko Tainio
 +
 +
# This model change the calculated 1-dimensional data to 2-dimensional so
 +
# that data can be downloaded e.g. to Analytica more easily.
 +
# Description of data:
 +
# - Intake fraction data: iF for European population due to European emissions
 +
 +
#The location of data files
 +
folder="N:/Huippuyksikko/Tutkimus/R19_Kopra/Mallit/IntakeFraction2/Result/European_emission/"
 +
 +
#Month/concentration - month 13 in annual mean
 +
for (i in 1:13) {
 +
 +
#Y must be defined before for -loop
 +
y = NULL
 +
 +
# Calculation for different countries. There are 48 sources in data.
 +
for (j in 1:48) {
 +
 +
# This reads from file the titles of countries to be used in calculations
 +
maa = read.table("N:/Huippuyksikko/Tutkimus/R19_Kopra/Mallit/IntakeFraction2/EUR_country.txt", sep="")
 +
maa = maa[,j]
 +
 +
# This read the intake fraction data files
 +
a = paste(folder,"2000/if_emis_eur_year_2000_country_",maa,"_month_",i,"_pop_all.txt", sep="")
 +
a = read.table(a, sep="")
 +
 +
#This combines all the data to 1-dimensional data
 +
y = c(y, a[,2])
 +
 +
}
 +
 +
#This change the 1-dimensional data to 2-dimensional
 +
# 40 is the number of populations used, 48 is number of emission sources (countries)
 +
dim(y) = c(40,48)
 +
 +
#Recording of results
 +
write.csv(y, paste(folder, file="comb_if_year_2000_country_all_month_",i,".csv", sep=""))
 +
 +
}
 +
 +
#Remove of all data from memory
 +
rm(folder, i, y, maa, a)
 +
</rcode>
 +
 +
'''Tainio_iF_emis_fin_pop_all_02.R
 +
 +
<rcode>
 +
#21.4.2008, Marko Tainio
 +
 +
# This model calculates intake fractions (iF) for Finnish
 +
# primary PM2.5 emissions. Description of model:
 +
# - Emission source: Finnish primary PM2.5 emissions from different sectors
 +
# - Dispersion: Dispersion over the Europe - small scale (11 countries)
 +
# - Population: European population (11 countries)
 +
 +
#Work folders - the location of data and where the output is recorded
 +
wrk_input = "N:/Huippuyksikko/Tutkimus/R19_Kopra/Data 2/"
 +
wrk_output = "N:/Huippuyksikko/Tutkimus/R19_Kopra/Mallit/IntakeFraction2/Result/"
 +
wrk_emis = "N:/Huippuyksikko/Tutkimus/R19_Kopra/Mallit/IntakeFraction2/"
 +
 +
#Breathing rate, 20 m^3/day -> unit g/s
 +
BR = 20/(24 * 60 * 60)
 +
 +
#Calculation of different years (2000, 2001, 2002)
 +
for (l in 0:2) {
 +
 +
#The calcultion of sectors - seven sectors are [DOM OTH TRA LPC LPP AGR ALL]
 +
for (k in 1:7) {
 +
 +
#Emission volume - this part copy emission data in unit g/s
 +
Q = paste(wrk_emis,"fin_emission_200",l,"_01.csv", sep="")
 +
Q = read.csv(Q, header=T)
 +
 +
#Concentration data
 +
#The data is read separately for different years because
 +
#year 2000 and 2001,2002 data has different dimensions
 +
rm(conc)
 +
if (l == 0) (conc = paste(wrk_input,"Concentration_primary_sector_200",l,"/monthly_PM_2_5/Fin_200",l,"_PM_hirlam_v3_8_1_", sep=""))
 +
if (l == 1) (conc = paste(wrk_input,"Concentration_primary_sector_200",l,"/monthly_200",l,"_PM_2_5/Fin_200",l,"_PM_EC_v3_8_1_", sep=""))
 +
if (l == 2) (conc = paste(wrk_input,"Concentration_primary_sector_200",l,"/monthly_200",l,"_PM_2_5/Fin_200",l,"_PM_EC_v3_8_1_", sep=""))
 +
 +
#The reading order for the data files
 +
if (k == 1) (conc = paste(conc,"_DOM_200",l,"12.txt", sep=""))
 +
if (k == 2) (conc = paste(conc,"_OTH_200",l,"12.txt", sep=""))
 +
if (k == 3) (conc = paste(conc,"_TRA_200",l,"12.txt", sep=""))
 +
if (k == 4) (conc = paste(conc,"_LPC_200",l,"12.txt", sep=""))
 +
if (k == 5) (conc = paste(conc,"_LPP_200",l,"12.txt", sep=""))
 +
if (k == 6) (conc = paste(conc,"_AGR_200",l,"12.txt", sep=""))
 +
if (k == 7) (conc = paste(conc,"all_srcs.txt", sep=""))
 +
conc = read.table(conc, header=T, sep="")
 +
 +
#Population data
 +
#The data is read separately for different years because
 +
#year 2000 and 2001,2002 data has different dimensions
 +
rm(pop)
 +
if (l == 0) (pop = paste(wrk_input, "Population_Combined/Pop_fin00_all.txt", sep=""))
 +
if (l == 1) (pop = paste(wrk_input, "Population_Combined/Pop_fin01_all.txt", sep=""))
 +
if (l == 2) (pop = paste(wrk_input, "Population_Combined/Pop_fin01_all.txt", sep=""))
 +
pop = read.table(pop, header=T, sep="")
 +
 +
#Co-ordination check-in - this check that co-ordinate system is same
 +
#in population and concentration data.
 +
check = (conc["lon"] - pop["Lo"]) + (conc["lat"] - pop["La"])
 +
sum(check)
 +
 +
#This part slice the concentration data from the data file
 +
#(that is, LO and LA data is sliced off)
 +
month = c("jan", "feb", "mar", "apr", "may", "jun", "jul", "aug", "sep", "oct", "nov", "dec", "annual")
 +
conc = data.frame(conc[month])
 +
 +
#Unit change in concentration data: kg/m^3 -> g/m^3
 +
conc = conc * 1000
 +
 +
#This part slice the population data from the data file
 +
#(that is, LO and LA data is sliced off)
 +
#The data is read separately for different years because
 +
#year 2000 and 2001,2002 data has different dimensions
 +
rm(country_fin)
 +
if (l == 0) (country_fin = c("Sum_ASUKKA", "Sum_Estoni", "Sum_Latvia", "Sum_Lithua", "Sum_Sweden", "Sum_Poland", "Sum_Denmar", "Sum_German", "Sum_Norway", "Sum_Belaru", "Sum_Russia", "Pop_all"))
 +
if (l == 1) (country_fin = c("Finland", "Sum_Estoni", "Sum_Latvia", "Sum_Lithua", "Sum_Poland", "Sum_Sweden", "Sum_Denmar", "Sum_German", "Sum_Norway", "Sum_Belaru", "Sum_Russia", "Pop_all"))
 +
if (l == 2) (country_fin = c("Finland", "Sum_Estoni", "Sum_Latvia", "Sum_Lithua", "Sum_Poland", "Sum_Sweden", "Sum_Denmar", "Sum_German", "Sum_Norway", "Sum_Belaru", "Sum_Russia", "Pop_all"))
 +
pop = data.frame(pop[country_fin])
 +
 +
#X must be defined before for -loop
 +
x = NULL
 +
 +
#Month/concentration - month 13 in annual mean
 +
for (i in 1:13) {
 +
 +
#Removal of huge values
 +
#There are some weird very high values in some concentration datasets.
 +
#This code replace those values with value 0.
 +
y = conc[,i]
 +
y = ifelse(y > 100, 0, y)
 +
 +
#Country/population - iF's are calculated separately for 12 different population
 +
for (j in 1:12) {
 +
 +
x[j] = sum(y * pop[,j])
 +
 +
}
 +
 +
#Calculation of iF
 +
x = (x * BR) / Q[k,i]
 +
 +
#Names for the data
 +
names(x) = c("FIN", "EST", "LAT", "LIT", "SWE", "POL", "DEN", "GER", "NOR", "BEL", "RUS", "All")
 +
 +
#Recording of results
 +
write.table(x, paste(wrk_output, file = "Finnish_emission/200",l,"/if_emis_fin_year_200",l,"_sector_",k,"_month_",i,".txt", sep=""), row.names = TRUE, col.names = FALSE)
 +
 +
}
 +
 +
}
 +
 +
}
 +
 +
rm(wrk_input, wrk_output, wrk_emis, BR, l, k, Q, conc, pop, check, month, country_fin, x, i, y, j)
 +
</rcode>
 +
 +
'''Tainio_iF_emis_fin_pop_all_data_conbination.R
 +
 +
<rcode>
 +
# 21.4.2008 Marko Tainio
 +
 +
# This model change the calculated 1-dimensional data to 2-dimensional so
 +
# that data can be downloaded e.g. to Analytica more easily.
 +
# Description of data:
 +
# - Intake fraction data: iF for European population due to Finnish emissions
 +
 +
#The location of data files
 +
folder="N:/Huippuyksikko/Tutkimus/R19_Kopra/Mallit/IntakeFraction2/Result/Finnish_emission/"
 +
 +
#Year (2000, 2001, 2002)
 +
for (k in 0:2) {
 +
 +
#The calcultion of sectors - seven sectors are [DOM OTH TRA LPC LPP AGR ALL]
 +
for (j in 1:7) {
 +
 +
#Y must be defined before for -loop
 +
y = NULL
 +
 +
#Month/concentration - month 13 in annual mean
 +
for (i in 1:13) {
 +
 +
# This read the intake fraction data files
 +
a = paste(folder,"200",k,"/if_emis_fin_year_200",k,"_sector_",j,"_month_",i,".txt", sep="")
 +
a = read.table(a, sep="")
 +
 +
#This combines all the data to 1-dimensional data
 +
y = c(y, a[,2])
 +
 +
}
 +
#This change the 1-dimensional data to 2-dimensional
 +
# 12 is the number of populations used, 13 is month (13 in annual mean)
 +
dim(y) = c(12,13)
 +
 +
#Recording of results
 +
write.csv(y, paste(folder, file="comb_if_year_200",k,"_sector_",j,"_month_all.csv", sep=""))
 +
 +
}
 +
 +
}
 +
 +
rm(folder, k, j, y, i, a)
 +
</rcode>
 +
 +
'''Tainio_population_01.R
 +
 +
<rcode>
 +
#7.3.2008 Marko Tainio
 +
 +
#This model summarise the population from different population datasets used
 +
#in intake fraction calculation
 +
 +
 +
#######
 +
#Data location for all datasets
 +
data_input = ("C:/Varsova/Kopra/Data 2/")
 +
data_output = ("C:/Varsova/Muutokset/Result/")
 +
 +
 +
#######
 +
#Finnish emissions, year 2000 population data
 +
pop_fin00_all = paste(data_input,"Population_Combined/Pop_fin00_all.txt", sep="")
 +
a_fin00 = read.table(pop_fin00_all, sep="", header=T)
 +
title_fin00 = c("Sum_ASUKKA", "Sum_Estoni", "Sum_Latvia", "Sum_Lithua", "Sum_Sweden", "Sum_Poland", "Sum_Denmar", "Sum_German", "Sum_Norway", "Sum_Belaru", "Sum_Russia", "Pop_all")
 +
a_fin00 = data.frame(a_fin00[title_fin00])
 +
 +
a = NULL
 +
 +
for (i in 1:12) {
 +
 +
a[i] = sum(a_fin00[,i])
 +
 +
}
 +
 +
names(a) = c("FIN", "EST", "LAT", "LIT", "SWE", "POL", "DEN", "GER", "NOR", "BEL", "RUS", "All")
 +
write.csv(a, paste(data_output, file="pop_fin00_sum.csv", sep=""))
 +
 +
rm(pop_fin00_all, a_fin00, title_fin00, a)
 +
 +
 +
#######
 +
#Finnish emissions, year 2001 population data
 +
pop_fin01_all = paste(data_input,"Population_Combined/Pop_fin01_all.txt", sep="")
 +
a_fin01 = read.table(pop_fin01_all, sep="", header=T)
 +
title_fin01 = c("Finland", "Sum_Estoni", "Sum_Latvia", "Sum_Lithua", "Sum_Sweden", "Sum_Poland", "Sum_Denmar", "Sum_German", "Sum_Norway", "Sum_Belaru", "Sum_Russia", "Pop_all")
 +
a_fin01 = data.frame(a_fin01[title_fin01])
 +
 +
a = NULL
 +
 +
for (i in 1:12) {
 +
 +
a[i] = sum(a_fin01[,i])
 +
 +
}
 +
 +
names(a) = c("FIN", "EST", "LAT", "LIT", "SWE", "POL", "DEN", "GER", "NOR", "BEL", "RUS", "All")
 +
write.csv(a, paste(data_output, file="pop_fin01_sum.csv", sep=""))
 +
 +
rm(pop_fin01_all, a_fin01, title_fin01, a)
 +
 +
 +
#######
 +
#European emissions, year 2000 population data
 +
pop_eur00_all = paste(data_input,"Population_Combined/Pop_Eur00_all.txt", sep="")
 +
a_eur00 = read.table(pop_eur00_all, sep="", header=T)
 +
title_eur00 = c("Sum_Austri", "Sum_Belgiu", "Sum_Cyprus", "Sum_Czech_", "Sum_Denmar", "Sum_Estoni", "Sum_France", "Sum_German", "Sum_Greece", "Sum_Hungar", "Sum_Irish_", "Sum_Italy", "Sum_Latvia", "Sum_Lithua", "Sum_Luxemb", "Sum_Malta", "Sum_Nether", "Sum_Poland", "Sum_Portug", "Sum_Romani", "Sum_Slovak", "Sum_Sloven", "Sum_Spain", "Sum_UK", "Sum_Sweden", "Sum_Bulgar", "Sum_Norway", "Sum_Belaru", "Sum_Russia", "Sum_Switze", "Sum_Ukrain", "Sum_Moldov", "Sum_Croati", "Sum_Serbia", "Sum_Albani", "Sum_Macedo", "Sum_Turkey", "pop_all")
 +
a_eur00 = data.frame(a_eur00[title_eur00])
 +
 +
a = NULL
 +
 +
for (i in 1:38) {
 +
 +
a[i] = sum(a_eur00[,i])
 +
 +
}
 +
 +
names(a) = c("Sum_Austri", "Sum_Belgiu", "Sum_Cyprus", "Sum_Czech_", "Sum_Denmar", "Sum_Estoni", "Sum_France", "Sum_German", "Sum_Greece", "Sum_Hungar", "Sum_Irish_", "Sum_Italy", "Sum_Latvia", "Sum_Lithua", "Sum_Luxemb", "Sum_Malta", "Sum_Nether", "Sum_Poland", "Sum_Portug", "Sum_Romani", "Sum_Slovak", "Sum_Sloven", "Sum_Spain", "Sum_UK", "Sum_Sweden", "Sum_Bulgar", "Sum_Norway", "Sum_Belaru", "Sum_Russia", "Sum_Switze", "Sum_Ukrain", "Sum_Moldov", "Sum_Croati", "Sum_Serbia", "Sum_Albani", "Sum_Macedo", "Sum_Turkey", "pop_all")
 +
write.csv(a, paste(data_output, file="pop_eur00_sum.csv", sep=""))
 +
 +
rm(pop_eur00_all, a_eur00, title_eur00, a)
 +
 +
 +
#European emissions, year 2000 population data, Finnish
 +
pop_eur00_fin = paste(data_input,"Population_Combined/Pop_Eur00_fin.txt", sep="")
 +
a_eur00 = read.table(pop_eur00_fin, sep="", header=T)
 +
a_eur00 = sum(a_eur00[,3])
 +
 +
write.csv(a_eur00, paste(data_output, file="pop_eur00_sum_fin.csv", sep=""))
 +
 +
rm(pop_eur00_fin, a_eur00)
 +
 +
#######
 +
 +
rm(data_input, data_output)
 +
</rcode>

Revision as of 14:12, 18 April 2017


Kopra was a research project about health impacts of fine particles in Finland.

Research question

Answer

Rationale

Data

←# : See N:\YMAL\Projects\R83_piltti filelist_kopra.txt for data list. --Jouni (talk) 14:12, 18 April 2017 (UTC)

Calculations

Marko's iF calculations

The codes are from U:\arkisto_kuopio\huippuyksikko\Tutkimus\R19_kopra\Mallit\IntakeFraction2. This section is for documentation only. The codes do not run.

Tainio_iF_emis_eur_pop_all_01.R

+ Show code

Tainio_iF_emis_eur_pop_all_data-conbination.R

+ Show code

Tainio_iF_emis_fin_pop_all_02.R

+ Show code

Tainio_iF_emis_fin_pop_all_data_conbination.R

+ Show code

Tainio_population_01.R

+ Show code

Retrieved from "http://www.opasnet.org/testiwiki/index.php?title=Kopra&oldid=40387"