Kopra
From Opasnet
Jump to navigation
Jump to search
This page is a study.
The page identifier is Op_en3052 |
|---|
| Moderator:Jouni (see all) |
|
|
| Upload data
|
Kopra was a research project about health impacts of fine particles in Finland.
Contents
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) (type: truth; paradigms: science: defence)
Calculations
Mika's iF calculations
- This section is for documentation only. The codes do not run.
- The codes are from U:\arkisto_kuopio\huippuyksikko\Tutkimus\R19_kopra\Mallit\Mika\Concentration_Background (date 3.3.2006).
IF_kaikki.r
- Input: Used data from ..\R19_Kopra\Data\Model_Converted\ but this data was not used in the final calculations.
- Output: Two files IFConcentration_Backround[5|1]0km.txt
# Datat -> N:\Huippuyksikko\Tutkimus\R19_Kopra\Data\Model_Converted
# PM2.5
# Concentration_Background
# Kaikki
pm25kaikki <- function(koko){
# 50x50 ja 10x10 ruudut & PM2.5
AlkuAika <- Sys.time()
SummaPixCi <- NULL
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "c:\\omat\\dataa\\Concentration_Backround\\PM2.5\\"
if(koko==50){
PopulaatioData <- NULL
PopulaatioData <-read.table(paste(polku,"Population_Europe\\50km\\ciesin_center_63lon_50km_2",sep=""))
}
if(koko==10){
PopulaatioData <- NULL
PopulaatioData <-read.table(paste(polku,"Population_Europe\\10km\\ciesin_center_63lon_10km_2",sep=""))
}
for(i in 1:12){
if(koko==50){
tiedosto <- paste("Concentration_Background\\2.5\\joinfilelist_month",i,".joined.out",sep="")
}
if(koko==10){
tiedosto <- paste("Concentration_Background\\2.5\\joinfilelist_month",i,".joined.out.1",sep="")
}
PitoisuusData <- NULL
PitoisuusData <- read.table(paste(polku, tiedosto ,sep=""))
PitoisuusMatriisi <- NULL
for(j in 4:dim(PitoisuusData)[2])
PitoisuusMatriisi <- cbind(PitoisuusMatriisi, PitoisuusData[,j])
apu <- NULL
apu <- matrix(1,dim(PitoisuusMatriisi)[2],1)
PitoisuusVektori <- PitoisuusMatriisi%*%apu
SummaPixCi[i] <- (t(PitoisuusVektori)%*%PopulaatioData[,3])
}
# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
# Q <- 3479.5265 Gg/vuosi
# Q <- 3479.5265 * 10^6/(365*24*60*60) kg/s
BR <- 20/(24 * 60 * 60)
# Luetaan joka maalle oma Q_i tiedostosta
# C:\omat\dataa\MaaKohtainenQ.csv ja lasketaan yhteen.
MaaKohtainenQ <- read.csv2("C:\\omat\\dataa\\MaaKohtainenQ.csv", header=FALSE)
Q <- sum(MaaKohtainenQ[,5])
IFConcentrationBackround <- NULL
for(i in 1:12)
IFConcentrationBackround[i] <- SummaPixCi[i]/(Q * 10^6/(365*24*60*60))* BR
if(koko==50){
TalletusKohde <- paste(talletuspolku,"IFConcentration_Backround50km.txt",sep="")
}
if(koko==10){
TalletusKohde <- paste(talletuspolku,"IFConcentration_Backround10km.txt",sep="")
}
write.table(IFConcentrationBackround, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
LoppuAika <- Sys.time()
LoppuAika - AlkuAika
#rm(MaaKohtainenQ,Q, SummaPixCi, apu, polku, talletuspolku, PopulaatioData, AlkuAika, LoppuAika, tiedosto,
#PitoisuusData, PitoisuusMatriisi, BR, IFConcentrationBackround, TalletusKohde)
}
# Datat -> N:\Huippuyksikko\Tutkimus\R19_Kopra\Data\Model_Converted
# PM2.5-10
# Concentration_Background
# Kaikki
pm2510kaikki <- function(koko){
# 50x50 ja 10x10 ruudut & PM2.5-10
AlkuAika <- Sys.time()
SummaPixCi <- NULL
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "c:\\omat\\dataa\\Concentration_Backround\\PM2.5-10\\"
if(koko==50){
PopulaatioData <- NULL
PopulaatioData <-read.table(paste(polku,"Population_Europe\\50km\\ciesin_center_63lon_50km_2",sep=""))
}
if(koko==10){
PopulaatioData <- NULL
PopulaatioData <-read.table(paste(polku,"Population_Europe\\10km\\ciesin_center_63lon_10km_2",sep=""))
}
for(i in 1:12){
if(koko==50){
tiedosto <- paste("Concentration_Background\\2.5-10\\joinfilelist_month",i,".joined.out",sep="")
}
if(koko==10){
tiedosto <- paste("Concentration_Background\\2.5-10\\joinfilelist_month",i,".joined.out.1",sep="")
}
PitoisuusData <- NULL
PitoisuusData <- read.table(paste(polku, tiedosto ,sep=""))
PitoisuusMatriisi <- NULL
for(j in 4:dim(PitoisuusData)[2])
PitoisuusMatriisi <- cbind(PitoisuusMatriisi, PitoisuusData[,j])
apu <- NULL
apu <- matrix(1,dim(PitoisuusMatriisi)[2],1)
PitoisuusVektori <- PitoisuusMatriisi%*%apu
SummaPixCi[i] <- (t(PitoisuusVektori)%*%PopulaatioData[,3])
}
# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
# Q <- 3479.5265 Gg/vuosi
# Q <- 3479.5265 * 10^6/(365*24*60*60) kg/s
BR <- 20/(24 * 60 * 60)
# Luetaan joka maalle oma Q_i tiedostosta
# C:\omat\dataa\MaaKohtainenQ.csv ja lasketaan yhteen.
MaaKohtainenQ <- read.csv2("C:\\omat\\dataa\\MaaKohtainenQ.csv", header=FALSE)
Q <- sum(MaaKohtainenQ[,6])
IFConcentrationBackround <- NULL
for(i in 1:12)
IFConcentrationBackround[i] <- SummaPixCi[i]/(Q * 10^6/(365*24*60*60))* BR
if(koko==50){
TalletusKohde <- paste(talletuspolku,"IFConcentration_Backround50km.txt",sep="")
}
if(koko==10){
TalletusKohde <- paste(talletuspolku,"IFConcentration_Backround10km.txt",sep="")
}
write.table(IFConcentrationBackround, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
LoppuAika <- Sys.time()
LoppuAika - AlkuAika
#rm(MaaKohtainenQ,Q, SummaPixCi, apu, polku, talletuspolku, PopulaatioData, AlkuAika, LoppuAika, tiedosto,
#PitoisuusData, PitoisuusMatriisi, BR, IFConcentrationBackround, TalletusKohde)
}
|
IFmaille_erikseen.r
- Input: Used data from ..\R19_Kopra\Data\Model_Converted\ but this data was not used in the final calculations.
- Output: Two files maakohtainen_IFConcentration_Backround[5|1]0km.txt
# Datat -> N:\Huippuyksikko\Tutkimus\R19_Kopra\Data\Model_Converted
# Luetaan ensin joka maalle oma Q_i tiedostosta
# C:\omat\dataa\MaaKohtainenQ.txt
#PM2.5
MaaKohtainenQ <- NULL
MaaKohtainenQ <- read.csv2("C:\\omat\\dataa\\MaaKohtainenQ.csv", header=FALSE)
Vaesto <- "Eurooppa"
HilaKoko <- 10
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
if(Vaesto == "Eurooppa"){
#Vaesto = Eurooppa
PopulaatioData <- NULL
PopulaatioData<-read.table(paste(polku,"Population_Europe\\",HilaKoko,"km\\ciesin_center_63lon_", HilaKoko,"km_2",sep=""))
}
if(Vaesto == "Suomi"){
#Vaesto = Suomi
PopulaatioData <- NULL
PopulaatioData<-read.table(paste(polku,"Population_Finland\\","fin_v00_",HilaKoko,"km.txt",sep=""))
}
apuPD <- NULL
apuPD <- PopulaatioData[PopulaatioData[,3]!=0,]
if(HilaKoko == 10) paate <- "out.1"
if(HilaKoko == 50) paate <- "out"
#Maakohtainen_Concentration_Background
AlkuAika <- Sys.time()
SummaPixCi <- NULL
talletuspolku <- "c:\\omat\\dataa\\Concentration_Backround\\PM2.5\\"
apu<-NULL
for(i in 1:12){
tiedosto<-paste("Concentration_Background\\2.5\\joinfilelist_month",i,".joined.",paate,sep="")
apuPitData <- NULL
apuPitData <- read.table(paste(polku, tiedosto,sep=""))
PitoisuusData <- NULL
PitoisuusData <- apuPitData[PopulaatioData[,3]!=0,]
# Yhdistetään Saksat Saksaksi ja Venäjän alueet Venäjäksi. Sekä poistetaan LO, LA ja All
PitoisuusMatriisi <- cbind(PitoisuusData[,4:18],PitoisuusData[,19] + PitoisuusData[,20],
PitoisuusData[,21:40], PitoisuusData[,41]+ PitoisuusData[,42]+PitoisuusData[,43]+PitoisuusData[,44],
PitoisuusData[,45:49])
apu <- (t(PitoisuusMatriisi)%*%apuPD[,3])
SummaPixCi <- rbind(SummaPixCi,t(apu))
}
# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
# Q <- 3479.5265 Gg/vuosi
# Q <- 3479.5265 * 10^6/(365*24*60*60) kg/s
BR <- 20/(24 * 60 * 60)
maakohtainenIFConcentrationBackround <- NULL
#PM2.5 on Q taulussa sarakkeessa 5
for(i in 1:12)
maakohtainenIFConcentrationBackround <-
rbind(maakohtainenIFConcentrationBackround,(SummaPixCi[i,]/(t(MaaKohtainenQ[,5]) * 10^6/(365*24*60*60))) * BR)
TalletusKohde <- paste(talletuspolku,"maakohtainen_IFConcentration_Backround_",HilaKoko,"km_",Vaesto,".txt",sep="")
write.table(maakohtainenIFConcentrationBackround, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
LoppuAika <- Sys.time()
LoppuAika - AlkuAika
rm(MaaKohtainenQ, SummaPixCi50km, apu, polku, talletuspolku, PopulaatioData, AlkuAika, LoppuAika, tiedosto,
PitoisuusData, PitoisuusMatriisi, BR, maakohtainenIFConcentrationBackround50km, TalletusKohde)
#----------------------------------------------------------------------
# Piirto PM2.5 maat erikseen
Vaesto <- "Suomi"
yAkseliYlaraja <- 4*10^(-7)
IFmat10<-
read.table(paste("c:\\omat\\dataa\\concentration_backround\\PM2.5\\maakohtainen_IFConcentration_Backround_10km_",
Vaesto,".txt",sep=""))
IFmat50<-
read.table(paste("c:\\omat\\dataa\\concentration_backround\\PM2.5\\maakohtainen_IFConcentration_Backround_50km_",
Vaesto,".txt",sep=""))
maat <- read.table("c:\\omat\\dataa\\maat.txt")
postscript(file=
paste("c:\\omat\\dataa\\concentration_backround\\PM2.5\\maaterikseen_",Vaesto,".ps",sep=""),
width=7,height=15,horizontal=F)
# par(mfrow = c(2,1))
for(i in 1:42){
maa <- paste(maat[i+1,1], maat[i+1,2],sep=", ")
plot(c(1:12),IFmat10[,i],type="b", xlim=c(1,12), ylim=c(0, yAkseliYlaraja),xlab="kuukausi",
ylab="IF 10km (punainen) ja 50km (musta)", col=rgb(1,0,0), main=maa)
par(new=T)
plot(c(1:12),IFmat50[,i],type="b", xlim=c(1,12), ylim=c(0, yAkseliYlaraja),xlab="kuukausi",
ylab="IF 10km (punainen) ja 50km (musta)", col=rgb(0,0,0), main=maa)
}
dev.off()
#---------------------------------------------------------
#---------------------
pm2510 <- function(koko){
# 50x50 ja 10x10 ruudut & PM2.5-10
AlkuAika <- Sys.time()
# Luetaan ensin joka maalle oma Q_i tiedostosta
# C:\omat\dataa\MaaKohtainenQ.txt
MaaKohtainenQ <- read.csv2("C:\\omat\\dataa\\MaaKohtainenQ.csv", header=FALSE)
SummaPixCi <- NULL
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "c:\\omat\\dataa\\Concentration_Backround\\PM2.5-10\\"
if(koko==50){
PopulaatioData <- NULL
PopulaatioData <-read.table(paste(polku,"Population_Europe\\50km\\ciesin_center_63lon_50km_2",sep=""))
}
if(koko==10){
PopulaatioData <- NULL
PopulaatioData <-read.table(paste(polku,"Population_Europe\\10km\\ciesin_center_63lon_10km_2",sep=""))
}
apu<-NULL
for(i in 1:12){
if(koko==50){
tiedosto <- paste("Concentration_Background\\2.5-10\\joinfilelist_month",i,".joined.out",sep="")
}
if(koko==10){
tiedosto <- paste("Concentration_Background\\2.5-10\\joinfilelist_month",i,".joined.out.1",sep="")
}
PitoisuusData <- NULL
PitoisuusData <- read.table(paste(polku, tiedosto ,sep=""))
# Yhdistetään Saksat Saksaksi ja Venäjän alueet Venäjäksi. Sekä poistetaan LO, LA ja All
PitoisuusMatriisi <- cbind(PitoisuusData[,4:18],PitoisuusData[,19] + PitoisuusData[,20],
PitoisuusData[,21:40], PitoisuusData[,41]+ PitoisuusData[,42]+PitoisuusData[,43]+PitoisuusData[,44],
PitoisuusData[,45:49])
apu <- (t(PitoisuusMatriisi)%*%PopulaatioData[,3])
SummaPixCi <- rbind(SummaPixCi,t(apu))
}
# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
# Q <- 3479.5265 Gg/vuosi
# Q <- 3479.5265 * 10^6/(365*24*60*60) kg/s
BR <- 20/(24 * 60 * 60)
maakohtainenIFConcentrationBackround <- NULL
for(i in 1:12)
maakohtainenIFConcentrationBackround <- rbind(maakohtainenIFConcentrationBackround,(SummaPixCi[i,]/(t(MaaKohtainenQ[,6]) * 10^6/(365*24*60*60))) * BR)
if(koko==50){
TalletusKohde <- paste(talletuspolku,"maakohtainen_IFConcentration_Backround50km.txt",sep="")
}
if(koko==10){
TalletusKohde <- paste(talletuspolku,"maakohtainen_IFConcentration_Backround10km.txt",sep="")
}
write.table(maakohtainenIFConcentrationBackround, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
LoppuAika <- Sys.time()
LoppuAika - AlkuAika
#rm(MaaKohtainenQ, SummaPixCi, apu, polku, talletuspolku, PopulaatioData, AlkuAika, LoppuAika, tiedosto,
#PitoisuusData, PitoisuusMatriisi, BR, maakohtainenIFConcentrationBackround, TalletusKohde)
}
#---------------------
|
- Code from U:\arkisto_kuopio\huippuyksikko\Tutkimus\R19_kopra\Mallit\Mika\Concentration_Finland (date 1.3.2006)
PM2_5Erikseen.r
- Input: Used data from ..\R19_Kopra\Data\Model_Converted\ but this data was not used in the final calculations.
- Output: Two files IFConcentration_Backround[5|1]0km.txt
# Datat -> N:\Huippuyksikko\Tutkimus\R19_Kopra\Data\Model_Converted
#PM2.5 Erikseen DOM, OTH, TRA, AGR,...
#Concentration_Finland
AlkuAika <- Sys.time()
PaastoLahde <- "AGR"
HilaKoko <- "10"
Vaesto <- "Eurooppa"
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- paste("c:\\omat\\dataa\\Concentration_Finland\\",PaastoLahde,"\\", sep="")
Qtaulu <- NULL
Qtaulu <- read.table(paste("N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Origin\\Emission_Silam_massbudget\\",
"Suomi","PM25Q.txt",sep=""))
if(Vaesto == "Eurooppa"){
#Vaesto = Eurooppa
PopulaatioData <- NULL
PopulaatioData<-read.table(paste(polku,"Population_Europe\\",HilaKoko,"km\\ciesin_center_63lon_", HilaKoko,"km_2",sep=""))
}
if(Vaesto == "Suomi"){
#Vaesto = Suomi
PopulaatioData <- NULL
PopulaatioData<-read.table(paste(polku,"Population_Finland\\","fin_v00_",HilaKoko,"km.txt",sep=""))
}
SummaPixCi <- NULL
apu <- NULL
Pitoisuus <- NULL
#Paasto <- NULL
IFConcentration <-NULL
TalletusKohde <- NULL
Q <- NULL
#tiedosto <- "Emission_grads\\fmsyke_data_joined_original_hila"
#Paasto <- read.table(paste(polku, tiedosto, sep=""))
#AGR:ssä vain 2 saraketta
#Q <- Paasto[,14] + Paasto[,20]
#Q <- Paasto[,12] + Paasto[,18] + Paasto[,24]
tiedosto <- paste("Concentration_Finland\\", HilaKoko,"km\\pm_0.1\\joinfilelist_PM_", PaastoLahde,".joined.out", sep="")
pm01 <- read.table(paste(polku, tiedosto, sep=""))
tiedosto <- paste("Concentration_Finland\\", HilaKoko,"km\\pm_0.1_1\\joinfilelist_PM_", PaastoLahde,".joined.out", sep="")
pm011 <- read.table(paste(polku, tiedosto, sep=""))
tiedosto <- paste("Concentration_Finland\\", HilaKoko,"km\\pm_1_2.5\\joinfilelist_PM_", PaastoLahde,".joined.out", sep="")
pm125 <- read.table(paste(polku, tiedosto, sep=""))
apu <- pm01 + pm011 + pm125
Pitoisuus <- cbind(Pitoisuus,apu[,3])
for( i in 2:9) Pitoisuus <- cbind(Pitoisuus,apu[,(i + 2 + 3)])
for( i in 10:12) Pitoisuus <- cbind(Pitoisuus,apu[,(i + 2 - 8)])
for(i in 1:12){
SummaPixCi[i] <- sum(PopulaatioData[,3]*Pitoisuus[,i])
}
# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
BR <- 20/(24 * 60 * 60)
#Q Kiloa sekunnissa
Q <- sum(Qtaulu[PaastoLahde==Qtaulu[,1],2:4])*(365/372)/(365*24*60*60)
IFConcentration <- SummaPixCi * BR/Q
TalletusKohde <- paste(talletuspolku,Vaesto,"_",HilaKoko,"km.txt",sep="")
write.table(IFConcentration, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
LoppuAika <- Sys.time()
LoppuAika - AlkuAika
#rm(apu, SummaPixCi50km, polku, talletuspolku, PopulaatioData, AlkuAika, LoppuAika, tiedosto,
#PitoisuusData, Q, BR, IFConcentration_FinlandAGR50km, TalletusKohde,PitoisuusAGR)
#-----------------------------------
# Kuvien piirto
#PaastoLahde <- "TRA"
piirtopolku <- paste("c:\\omat\\dataa\\Concentration_Finland\\",PaastoLahde,"\\", sep="")
postscript(file=paste(piirtopolku,"IF_kuva_",Vaesto,".ps", sep=""), width=10,height=10)#,horizontal=F)
IFmat10<-read.table(paste(piirtopolku,Vaesto,"_10km.txt", sep=""))
IFmat50<-read.table(paste(piirtopolku,Vaesto,"_50km.txt", sep=""))
plot(c(1:12),IFmat10[,1],type="b", xlim=c(1,12), ylim=c(0,1.2*10^(-6)),xlab="kuukausi", ylab="IF 10km ja 50 km (50km on musta)", col=rgb(1,0,0))
par(new=T)
plot(c(1:12),IFmat50[,1],type="b", xlim=c(1,12), ylim=c(0,1.2*10^(-6)),xlab="kuukausi", ylab="IF 10km ja 50 km (50km on musta)", col=rgb(0,0,0))
# musta on 0,0,0
dev.off()
|
- Codes from U:\arkisto_kuopio\huippuyksikko\Tutkimus\R19_kopra\Mallit\Mika\Concentration_Point (dates 10.11.2005 - 16.6.2006)
R_Kopra.r (10.11.2005)
# Datat -> N:\Huippuyksikko\Tutkimus\R19_Kopra\Data\Model_Converted
#--------------
# Concentration_Point
# 50x50 ruudut
AlkuAika <- Sys.time()
SummaPixCi50km <- matrix(nrow=12,ncol=18)
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "c:\\omat\\dataa\\concentration_point\\"
PopulaatioData <- NULL
PopulaatioData<-read.table(paste(polku,"Population_Europe\\50km\\ciesin_center_63lon_50km_2",sep=""))
for(i in 1:12){
tiedosto1<-paste("joinfilelist_",i,"_PM_0_1.joined.out",sep="")
tiedosto2<-paste("joinfilelist_",i,"_PM_0_1_1.joined.out",sep="")
tiedosto3<-paste("joinfilelist_",i,"_PM_1_2_5.joined.out",sep="")
PitoisuusData1 <- NULL
PitoisuusData1 <-read.table(paste(polku,"Concentration_Point\\",tiedosto1,sep=""))
PitoisuusData2 <- NULL
PitoisuusData2 <-read.table(paste(polku,"Concentration_Point\\",tiedosto2,sep=""))
PitoisuusData3 <- NULL
PitoisuusData3 <-read.table(paste(polku,"Concentration_Point\\",tiedosto3,sep=""))
PM2.5PitoisuusData <- NULL
PM2.5PitoisuusData <- PitoisuusData1 + PitoisuusData2 + PitoisuusData3
for(j in 1:18)
SummaPixCi50km[i,j] <- sum( PM2.5PitoisuusData[,(j+2)] * PopulaatioData[,3] )
}
# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
BR <- 20/(24 * 60 * 60)
# Q <- 3 * 1/(365 * 24 * 60 * 60)kg/s
Q <- 3/(365 * 24 * 60 * 60)
IFConcentrationPoint50km <- SummaPixCi50km * BR/Q
TalletusKohde <- paste(talletuspolku,"IFConcentrationPoint50km.txt",sep="")
write.table(IFConcentrationPoint50km, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
LoppuAika <- Sys.time()
LoppuAika - AlkuAika
rm(SummaPixCi50km, polku, talletuspolku, PopulaatioData, AlkuAika, LoppuAika, tiedosto1, tiedosto2, tiedosto3, PitoisuusData1, PitoisuusData2, PitoisuusData3, PM2.5PitoisuusData, Q, BR, IFConcentrationPoint50km, TalletusKohde)
#-----------------------------
# 10x10 ruudut
AlkuAika <- Sys.time()
SummaPixCi10km <- matrix(nrow=12,ncol=18)
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "c:\\omat\\dataa\\concentration_point\\"
PopulaatioData <- NULL
PopulaatioData<-read.table(paste(polku,"Population_Europe\\10km\\ciesin_center_63lon_10km_2",sep=""))
for(i in 1:12){
tiedosto1<-paste("joinfilelist_",i,"_PM_0_1.joined.out.1",sep="")
tiedosto2<-paste("joinfilelist_",i,"_PM_0_1_1.joined.out.1",sep="")
tiedosto3<-paste("joinfilelist_",i,"_PM_1_2_5.joined.out.1",sep="")
PitoisuusData1 <- NULL
PitoisuusData1 <-read.table(paste(polku,"Concentration_Point\\",tiedosto1,sep=""))
PitoisuusData2 <- NULL
PitoisuusData2 <-read.table(paste(polku,"Concentration_Point\\",tiedosto2,sep=""))
PitoisuusData3 <- NULL
PitoisuusData3 <-read.table(paste(polku,"Concentration_Point\\",tiedosto3,sep=""))
PM2.5PitoisuusData <- NULL
PM2.5PitoisuusData <- PitoisuusData1 + PitoisuusData2 + PitoisuusData3
for(j in 1:18)
SummaPixCi10km[i,j] <- sum( PM2.5PitoisuusData[,(j+2)] * PopulaatioData[,3] )
}
# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
BR <- 20/(24 * 60 * 60)
# Q <- 3 * 1/(365 * 24 * 60 * 60)kg/s
Q <- 3/(365 * 24 * 60 * 60)
IFConcentrationPoint10km <- SummaPixCi10km * BR/Q
TalletusKohde <- paste(talletuspolku,"IFConcentrationPoint10km.txt",sep="")
write.table(IFConcentrationPoint10km, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
LoppuAika <- Sys.time()
LoppuAika - AlkuAika
rm(SummaPixCi10km, polku, talletuspolku, PopulaatioData, AlkuAika, LoppuAika, tiedosto1, tiedosto2, tiedosto3, PitoisuusData1, PitoisuusData2, PitoisuusData3, PM2.5PitoisuusData, Q, BR, IFConcentrationPoint10km, TalletusKohde)
|
Distance_piirto.r (25.1.2006)
#-----------------------------------------------------------------------
#Maan pinnalla sijaitsevien pisteiden etäisyyden laskeva funktio
#-----------------------------------------------------------------------
#maan säde n. 6371 km
Etaisyys <- function(sade, p1LO, p1LA, p2LO, p2LA)
{
piste1LA <- p1LA * (2*pi/360)
piste1LO <- p1LO * (2*pi/360)
piste2LA <- p2LA * (2*pi/360)
piste2LO <- p2LO * (2*pi/360)
#pallokoorinaatistoa käyttäen
SuoraValiMatka <- sade * sqrt( (cos(piste1LO)*cos(piste1LA) - cos(piste2LO)*cos(piste2LA))^2
+ (sin(piste1LO)*cos(piste1LA) - sin(piste2LO)*cos(piste2LA))^2
+ (sin(piste1LA) - sin(piste2LA))^2 )
#cosini lause
Kulma <- acos(1 - SuoraValiMatka^2/(2*sade^2))
#kaaren pituus on kulma(radiaaneina)*säde
sade * Kulma
}
#----------------------------------------------------------------------------
# PopulaatioData siirtävä funktio
#-------------------------------------------------------------------------
LohkoPituus <- 354
Lohkoja <- dim(PopulaatioData)[1]/LohkoPituus
# PD <- PDSiirto(0,0,2,0,PopulaatioData, Lohkoja, LohkoPituus)
PDSiirto <- function(poh, ita, ete, lan, PopulaatioData, Lohkoja, LohkoPituus){
UusiPD <- NULL
UusiPD <- PopulaatioData
#Siirto etelään
if(ete > 0){
#Siirretään dataa
for( i in 1:(Lohkoja-ete)){
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <-
PopulaatioData[((i + ete - 1)*LohkoPituus + 1):((i+ete)*LohkoPituus),3]
}
#Täytetään lopusta nollilla
for(i in 1:ete){
UusiPD[((Lohkoja-ete)*LohkoPituus+1):((Lohkoja - ete +1)*LohkoPituus),3] <- t(t(rep(0,LohkoPituus)))
}
}
if(lan > 0){
#Siirto länteen
#Siirretään dataa
for( i in 1:Lohkoja){
apu1 <- NULL
apu1 <- PopulaatioData[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3]
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <-
t(t(c(apu1[(1+lan):LohkoPituus],rep(0,lan))))
}
}
if(poh > 0){
#Siirto pohjoiseen
#Siirretään dataa
for( i in (poh + 1):(Lohkoja)){
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <-
PopulaatioData[((i-poh)*LohkoPituus + 1):((i-poh+1)*LohkoPituus),3]
}
#Täytetään lopusta nollilla
for(i in 1:poh){
UusiPD[((i-1)*LohkoPituus+1):((i)*LohkoPituus),3] <- t(t(rep(0,LohkoPituus)))
}
}
if(ita > 0){
#Siirto itään
#Siirretään dataa
for( i in 1:Lohkoja){
apu1 <- NULL
apu1 <- PopulaatioData[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3]
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <-
t(t(c(rep(0,ita),apu1[1:(LohkoPituus - ita)])))
}
}
UusiPD
}
#-----------------------------------------------------------
#Lähimpien Ruutujen etsintä
#-----------------------------------------------------------
# Annetaan pisteen Lo ja La koordinaatit, piste itse, Populaatio data ja Pienhiukkasdata.
# Viimeisin annetaan, jotta ruutujen määrää voidaan vähentää.
#-----------------------------------------------------------
LahimmatRuudut <- function(LO,LA, PD, PM, Piste){
LOLisays <- 0.1981
LOLisaysArvio <- 0.199
LALisays <- 0.08995
LALisaysArvio <- 0.09
apuPD <- NULL
apuPD <- PD[PD[,3]!=0 & (PM[,(1 + 2*Piste)]!=0 | PM[,(2 + 2*Piste)]!=0),]
apu1 <- NULL
#Sade 50 km
for(k in 1:dim(apuPD)[1]){
if( (apuPD[k,1] >= LO & apuPD[k,2] >= LA) | (apuPD[k,1] < LO & apuPD[k,2] < LA ) ){
apu1[k] <-
(Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2])) <= 50 &
Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2] + LALisays)) <= 50)
}
if((apuPD[k,1] < LO & apuPD[k,2] >= LA) | (apuPD[k,1] >= LO & apuPD[k,2] < LA ) ){
apu1[k] <-
(Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2])) <= 50 &
Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2] + LALisays)) <= 50)
}
}
apuPD <- apuPD[apu1,]
apu2 <- NULL
apu2 <- rep(FALSE,(dim(apuPD)[1]))
for(k in 1:dim(apuPD)[1]){
if( (apuPD[k,1] <= LO & apuPD[k,2] <= LA) &
(LO < (apuPD[k,1] + LOLisays) & LA < (apuPD[k,2] + LALisays)) )
{
apu2[k] <- (1==1)
Ruutu <- k
}
}
for(k in 1:dim(apuPD)[1]){
if
(
((apuPD[Ruutu,1]-LOLisaysArvio)<= apuPD[k,1]) & (apuPD[k,1] <= (apuPD[Ruutu,1]+LOLisaysArvio)) &
((apuPD[Ruutu,2]-LALisaysArvio) <= apuPD[k,2]) & (apuPD[k,2] <= (apuPD[Ruutu,2]+LALisaysArvio))
)
{
apu2[k] <- (1==1)
if(length(apu2[apu2]) == 9){ k <- dim(apuPD)[1] + 1}
}
}
apu2
}
#------------------------------------------------------
# IF-datan piirtävä koodi
#------------------------------------------------
#polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
polku <- "c:\\omat\\dataa\\concentration_point\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
SadeLisays <- 50
PieninSade <- 50
MaxSade <- 1500
for(i in 1:1){
Piste <- i
#data <- paste(polku,"IF_Piste_",Piste,"_Distance.txt",sep="")
data <- paste(polku,"IF_Piste_",Piste,"_Distance_Siirretty_kk3.txt",sep="")
#data <- paste(polku,"IF_Piste_",Piste,"_Distance_alle_50.txt",sep="")
#kohde <- paste(talletuspolku, "Piste_",Piste,".ps",sep="")
kohde <- paste(talletuspolku, "Piste_",Piste,"_Siirretty_kk3.ps",sep="")
#kohde <- paste(talletuspolku, "Piste_",Piste,"_alle_50.ps",sep="")
postscript(file=kohde, width=10,height=10,horizontal=F)
IFDist<-read.table(data)
xAkseli <- seq(PieninSade, MaxSade, by=SadeLisays)
for(kk in 1:1){
yAkseliHI <- NULL
yAkseliLO <- NULL
for(y in 1:(dim(IFDist)[2]/2)){
yAkseliHI[y] <- IFDist[kk,(2*y - 1)]
yAkseliLO[y] <- IFDist[kk,(2*y)]
}
plot(xAkseli,yAkseliHI, xlim=c(1,(MaxSade+ 2* SadeLisays)), ylim=c(0, 2.5*10^(-6)),
xlab="Säde", ylab="IF (high = punainen, low = musta)", col=rgb(1,0,0), main = paste("Piste ", Piste, " kuukausi ", kk, sep=""))
par(new=T)
plot(xAkseli,yAkseliLO, xlim=c(1,(MaxSade+ 2* SadeLisays)), ylim=c(0, 2.5*10^(-6)),
xlab="Säde", ylab="IF (high = punainen, low = musta)", col=rgb(0,0,0), main = paste("Piste ", Piste, " kuukausi ", kk, sep=""))
}
# musta on 0,0,0
dev.off()
}
#-------------------------------------------------------------------------------------------------------------
#Perustaso
DataRuudut <- 10
# HUOM! tiedostoissa pääte .out on 50km ja .out.1 on 10km
if (DataRuudut == 10) paate <- ".out.1"
if (DataRuudut == 50) paate <- ".out"
SadeLisays <- 50
PieninSade <- 50
MaxSade <- 50
LOLisays <- 0.1981
LOLisaysArvio <- 0.199
LALisays <- 0.08995
LALisaysArvio <- 0.09
PisteKK <- c(25,60,22,60.5,22,63,26,62,30,63,24.5,65,29.5,65,24,67.5,28,69)
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
PopulaatioData <- NULL
PopulaatioData <-
read.table(paste(polku,"Population_Europe\\",DataRuudut,"km\\ciesin_center_63lon_",DataRuudut,"km_2",sep=""))
for( Piste in 1:1){
LO <- PisteKK[(2*Piste-1)]
LA <- PisteKK[(2*Piste)]
SummaPixCi <- NULL
# Seuraavassa loopissa kk on kuukausi
for(kk in 1:12){
tiedosto1<-paste("joinfilelist_",kk,"_PM_0_1.joined", paate,sep="")
tiedosto2<-paste("joinfilelist_",kk,"_PM_0_1_1.joined", paate,sep="")
tiedosto3<-paste("joinfilelist_",kk,"_PM_1_2_5.joined", paate,sep="")
PitoisuusData1 <- NULL
PitoisuusData1 <-read.table(paste(polku,"Concentration_Point\\",tiedosto1,sep=""))
PitoisuusData2 <- NULL
PitoisuusData2 <-read.table(paste(polku,"Concentration_Point\\",tiedosto2,sep=""))
PitoisuusData3 <- NULL
PitoisuusData3 <-read.table(paste(polku,"Concentration_Point\\",tiedosto3,sep=""))
PM2.5PitoisuusData<- NULL
PM2.5PitoisuusData <- PitoisuusData1 + PitoisuusData2 + PitoisuusData3
apuPM <- NULL
apuPM <- PM2.5PitoisuusData[(PM2.5PitoisuusData[,(1 + 2*Piste)]!=0 | PM2.5PitoisuusData[,(2 + 2*Piste)]!=0)
& PopulaatioData[,3]!=0,]
apuPD <- NULL
apuPD <- PopulaatioData[PopulaatioData[,3]!=0 &
(PM2.5PitoisuusData[,(1 + 2*Piste)]!=0 | PM2.5PitoisuusData[,(2 + 2*Piste)]!=0),]
apuHI <- NULL
apuLO <- NULL
apuSumma <- NULL
#apu1 <- rep(FALSE, dim(apuPD)[1])
for(Sade in seq(MaxSade, PieninSade, by= -SadeLisays )){
apu1 <- NULL
for(k in 1:dim(apuPD)[1]){
if( (apuPD[k,1] >= LO & apuPD[k,2] >= LA) | (apuPD[k,1] < LO & apuPD[k,2] < LA ) ){
apu1[k] <-
(Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2])) < Sade &
Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2] + LALisays)) < Sade)
}
if((apuPD[k,1] < LO & apuPD[k,2] >= LA) | (apuPD[k,1] >= LO & apuPD[k,2] < LA ) ){
apu1[k] <-
(Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2])) < Sade &
Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2] + LALisays)) < Sade)
}
}
apuPM <- apuPM[apu1,]
apuPD <- apuPD[apu1,]
}
#-------------------------------------------------------
#Viimeisten ruutujen etsintä
apu2 <- NULL
apu2 <- LahimmatRuudut(LO,LA, apuPD, apuPM, Piste)
apuPM <- apuPM[apu2,]
apuPD <- apuPD[apu2,]
#------------------------------------------------------------------
apuHI[1] <- sum( apuPM[,(1 + 2*Piste)] * apuPD[,3] )
apuLO[1] <- sum( apuPM[,(2 + 2*Piste)] * apuPD[,3] )
apuSumma[1] <- apuHI[1]
apuSumma[2] <- apuLO[1]
SummaPixCi <- rbind(SummaPixCi, t(apuSumma))
}
# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
BR <- 20/(24 * 60 * 60)
# Q <- 3 * 1/(365 * 24 * 60 * 60)kg/s
Q <- 3/(365 * 24 * 60 * 60)
IF <- SummaPixCi * BR/Q
TalletusKohde <- paste(talletuspolku,"IF_Piste_",Piste,"_Perus_Distance.txt",sep="")
write.table(IF, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
}
|
Distance.r (17.2.2006) Codes varaDistance.r (31.1.2006) and varmuuskopioDistance.r (13.2.3006) omitted.
#-----------------------------------------------------------------------
#Maan pinnalla sijaitsevien pisteiden etäisyyden laskeva funktio
#-----------------------------------------------------------------------
#maan säde n. 6371 km
#jos laskutapa on 1, niin lasketaan 2 pisteen välimatka pallon pintaa pitkin
# jos laskutapa =2, niin lasketaan pisteiden (p1LO, p1LA) ja (p2LO, p2LA) euklidinen etäisyys,
# tällöin säteellä ei ole merkitystä
Etaisyys <- function(sade, p1LO, p1LA, p2LO, p2LA, laskutapa)
{
if(laskutapa == 1){
piste1LA <- p1LA * (2*pi/360)
piste1LO <- p1LO * (2*pi/360)
piste2LA <- p2LA * (2*pi/360)
piste2LO <- p2LO * (2*pi/360)
#pallokoorinaatistoa käyttäen
SuoraValiMatka <- sade * sqrt( (cos(piste1LO)*cos(piste1LA) - cos(piste2LO)*cos(piste2LA))^2
+ (sin(piste1LO)*cos(piste1LA) - sin(piste2LO)*cos(piste2LA))^2
+ (sin(piste1LA) - sin(piste2LA))^2 )
#cosini lause
Kulma <- acos(1 - SuoraValiMatka^2/(2*sade^2))
#kaaren pituus on kulma(radiaaneina)*säde
tulos <- sade * Kulma
}
if(laskutapa == 2)
{
tulos <- sqrt((p1LO - p2LO)^2 + (p1LA - p2LA)^2)
}
tulos
}
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
PopulaatioData <- NULL
PopulaatioData <-
read.table(paste(polku,"Population_Europe\\",DataRuudut,"km\\ciesin_center_63lon_",DataRuudut,"km_2",sep=""))
#----------------------------------------------------------------------------
# PopulaatioData siirtävä funktio
#-------------------------------------------------------------------------
LohkoPituus <- 354
Lohkoja <- dim(PopulaatioData)[1]/LohkoPituus
# sPD <- PDSiirto(0,0,2,0,PopulaatioData, Lohkoja, LohkoPituus)
PDSiirto <- function(poh, ita, ete, lan, PopulaatioData, Lohkoja, LohkoPituus){
UusiPD <- NULL
UusiPD <- PopulaatioData
#Siirto etelään
if(ete > 0){
#Siirretään dataa
for( i in 1:(Lohkoja-ete)){
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <-
PopulaatioData[((i + ete - 1)*LohkoPituus + 1):((i+ete)*LohkoPituus),3]
}
#Täytetään lopusta nollilla
for(i in 1:ete){
UusiPD[((Lohkoja-ete)*LohkoPituus+1):((Lohkoja - ete +1)*LohkoPituus),3] <- t(t(rep(0,LohkoPituus)))
}
}
if(lan > 0){
#Siirto länteen
#Siirretään dataa
for( i in 1:Lohkoja){
apu1 <- NULL
apu1 <- PopulaatioData[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3]
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <-
t(t(c(apu1[(1+lan):LohkoPituus],rep(0,lan))))
}
}
if(poh > 0){
#Siirto pohjoiseen
#Siirretään dataa
for( i in (poh + 1):(Lohkoja)){
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <-
PopulaatioData[((i-poh)*LohkoPituus + 1):((i-poh+1)*LohkoPituus),3]
}
#Täytetään lopusta nollilla
for(i in 1:poh){
UusiPD[((i-1)*LohkoPituus+1):((i)*LohkoPituus),3] <- t(t(rep(0,LohkoPituus)))
}
}
if(ita > 0){
#Siirto itään
#Siirretään dataa
for( i in 1:Lohkoja){
apu1 <- NULL
apu1 <- PopulaatioData[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3]
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <-
t(t(c(rep(0,ita),apu1[1:(LohkoPituus - ita)])))
}
}
UusiPD
}
PData <- LahimmatRuudut(25,60.08,PopulaatioData,1, 58, 65)
apuri <- PopulaatioData[58 <= PopulaatioData[,2] & PopulaatioData[,2] <= 65, ]
apuri[PData,]
PData
#-----------------------------------------------------------
#Lähimpien Ruutujen etsintä
#-----------------------------------------------------------
# Annetaan pisteen Lo ja La koordinaatit, piste itse, Populaatio data ja Pienhiukkasdata.
# Viimeisin annetaan, jotta ruutujen määrää voidaan vähentää.
# Voidaan myös rajata antamalla alaraja ja yläraja Latitudelle
#-----------------------------------------------------------
LahimmatRuudut <-function(LO,LA, PD, Sade,Sateella, LAalaraja, LAylaraja){
#function(LO,LA, PD, PM, Piste){
LOLisays <- 0.1981
LOLisaysArvio <- 0.199
LALisays <- 0.08995
LALisaysArvio <- 0.09
apuPD <- NULL
#apuPD <- PD[PD[,3]!=0 & (PM[,(1 + 2*Piste)]!=0 | PM[,(2 + 2*Piste)]!=0),]
#apuPD <- PD[PD[,3]!=0 & LAalaraja <= PD[,2] & PD[,2] <= LAylaraja,]
apuPD <- PD[LAalaraja <= PD[,2] & PD[,2] <= LAylaraja,]
apu2 <- NULL
if(Sateella==1){
for(k in 1:dim(apuPD)[1]){
if( (apuPD[k,1] >= LO & apuPD[k,2] >= LA) | (apuPD[k,1] < LO & apuPD[k,2] < LA ) ){
apu2[k] <-
(Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2]),1) <= Sade &
Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2] + LALisays),1) <= Sade)
}
if((apuPD[k,1] < LO & apuPD[k,2] >= LA) | (apuPD[k,1] >= LO & apuPD[k,2] < LA ) ){
apu2[k] <-
(Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2]),1) <= Sade &
Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2] + LALisays),1) <= Sade)
}
}
#apuPD <- apuPD[apu2,]
}
if(Sateella==0)
{
apu2 <- rep(FALSE,(dim(apuPD)[1]))
for(k in 1:dim(apuPD)[1]){
if( (apuPD[k,1] <= LO & apuPD[k,2] <= LA) &
(LO < (apuPD[k,1] + LOLisays) & LA < (apuPD[k,2] + LALisays)) )
{
apu2[k] <- (1==1)
Ruutu <- k
}
}
for(k in 1:dim(apuPD)[1]){
if
(
((apuPD[Ruutu,1]-LOLisaysArvio)<= apuPD[k,1]) & (apuPD[k,1] <= (apuPD[Ruutu,1]+LOLisaysArvio)) &
((apuPD[Ruutu,2]-LALisaysArvio) <= apuPD[k,2]) & (apuPD[k,2] <= (apuPD[Ruutu,2]+LALisaysArvio))
)
{
apu2[k] <- (1==1)
if(length(apu2[apu2]) == 9){ k <- dim(apuPD)[1] + 1}
}
}
}
apuPD[apu2,]
}
#------------------------------------------------------
# IF-datan piirtävä koodi
#------------------------------------------------
#polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
#polku <- "c:\\omat\\dataa\\concentration_point\\"
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
SadeLisays <- 50
PieninSade <- 50
MaxSade <- 1700
for(i in 1:1){
Piste <- i
data <- paste(polku,"IF_Piste_",Piste,"_Distance.txt",sep="")
#data <- paste(polku,"IF_Piste_",Piste,"_Distance_alle_50.txt",sep="")
#data <- paste(polku,"IF_Piste_",Piste,"_Distance_Siirretty.txt",sep="")
#kohde <- paste(talletuspolku, "Piste_",Piste,".ps",sep="")
kohde <- paste(talletuspolku, "Piste_",Piste,"_Siirretty.ps",sep="")
#kohde <- paste(talletuspolku, "Piste_",Piste,"_alle_50.ps",sep="")
postscript(file=kohde, width=10,height=10,horizontal=F)
IFDist<-read.table(data)
xAkseli <- seq(PieninSade, MaxSade, by=SadeLisays)
for(kk in 1:12){
yAkseliHI <- NULL
yAkseliLO <- NULL
for(y in 1:(dim(IFDist)[2]/2)){
yAkseliHI[y] <- IFDist[kk,(2*y - 1)]
yAkseliLO[y] <- IFDist[kk,(2*y)]
}
plot(xAkseli,yAkseliHI, xlim=c(1,(MaxSade+ 2* SadeLisays)), ylim=c(0, 5*10^(-6)),
xlab="Säde", ylab="IF (high = punainen, low = musta)", col=rgb(1,0,0), main = paste("Piste ", Piste, " kuukausi ", kk, sep=""))
par(new=T)
plot(xAkseli,yAkseliLO, xlim=c(1,(MaxSade+ 2* SadeLisays)), ylim=c(0, 5*10^(-6)),
xlab="Säde", ylab="IF (high = punainen, low = musta)", col=rgb(0,0,0), main = paste("Piste ", Piste, " kuukausi ", kk, sep=""))
}
# musta on 0,0,0
dev.off()
}
#-------------------------------------------------------------------------------------------------------------
DataRuudut <- 10
# HUOM! tiedostoissa pääte .out on 50km ja .out.1 on 10km
if (DataRuudut == 10) paate <- ".out.1"
if (DataRuudut == 50) paate <- ".out"
SadeLisays <- 50
PieninSade <- 50
MaxSade <- 1700
LOLisays <- 0.1981
LOLisaysArvio <- 0.199
LALisays <- 0.08995
LALisaysArvio <- 0.09
PisteKK <- c(25,60,22,60.5,22,63,26,62,30,63,24.5,65,29.5,65,24,67.5,28,69)
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
PopulaatioData <- NULL
PopulaatioData <-
read.table(paste(polku,"Population_Europe\\",DataRuudut,"km\\ciesin_center_63lon_",DataRuudut,"km_2",sep=""))
PopulaatioData <- sPD
for( Piste in 1:1){
LO <- PisteKK[(2*Piste-1)]
LA <- PisteKK[(2*Piste)]
SummaPixCi <- NULL
# Seuraavassa loopissa kk on kuukausi
for(kk in 1:12){
tiedosto1<-paste("joinfilelist_",kk,"_PM_0_1.joined", paate,sep="")
tiedosto2<-paste("joinfilelist_",kk,"_PM_0_1_1.joined", paate,sep="")
tiedosto3<-paste("joinfilelist_",kk,"_PM_1_2_5.joined", paate,sep="")
PitoisuusData1 <- NULL
PitoisuusData1 <-read.table(paste(polku,"Concentration_Point\\",tiedosto1,sep=""))
PitoisuusData2 <- NULL
PitoisuusData2 <-read.table(paste(polku,"Concentration_Point\\",tiedosto2,sep=""))
PitoisuusData3 <- NULL
PitoisuusData3 <-read.table(paste(polku,"Concentration_Point\\",tiedosto3,sep=""))
PM2.5PitoisuusData<- NULL
PM2.5PitoisuusData <- PitoisuusData1 + PitoisuusData2 + PitoisuusData3
apuPM <- NULL
apuPM <- PM2.5PitoisuusData[(PM2.5PitoisuusData[,(1 + 2*Piste)]!=0 | PM2.5PitoisuusData[,(2 + 2*Piste)]!=0)
& PopulaatioData[,3]!=0,]
apuPD <- NULL
apuPD <- PopulaatioData[PopulaatioData[,3]!=0 &
(PM2.5PitoisuusData[,(1 + 2*Piste)]!=0 | PM2.5PitoisuusData[,(2 + 2*Piste)]!=0),]
apuHI <- NULL
apuLO <- NULL
apuSumma <- NULL
#apu1 <- rep(FALSE, dim(apuPD)[1])
sadeindeksi <- 0
for(Sade in seq(MaxSade, PieninSade, by= -SadeLisays )){
apu1 <- NULL
sadeindeksi <- sadeindeksi + 1
for(k in 1:dim(apuPD)[1]){
if( (apuPD[k,1] >= LO & apuPD[k,2] >= LA) | (apuPD[k,1] < LO & apuPD[k,2] < LA ) ){
apu1[k] <-
(Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2]),1) < Sade &
Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2] + LALisays),1) < Sade)
}
if((apuPD[k,1] < LO & apuPD[k,2] >= LA) | (apuPD[k,1] >= LO & apuPD[k,2] < LA ) ){
apu1[k] <-
(Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2]),1) < Sade &
Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2] + LALisays),1) < Sade)
}
}
apuPM <- apuPM[apu1,]
apuPD <- apuPD[apu1,]
apuHI[sadeindeksi] <- sum( apuPM[,(1 + 2*Piste)] * apuPD[,3] )
apuLO[sadeindeksi] <- sum( apuPM[,(2 + 2*Piste)] * apuPD[,3] )
aputaulu <- c(Sade, dim(apuPD)[1])
TalletusKohde <- paste(talletuspolku,"Piste_",Piste,"\\IF_Piste_",Piste,"_kk_",kk,"_Sade_",Sade,"_Distance.txt",sep="")
write.table(t(aputaulu), file = TalletusKohde, row.names = FALSE, col.names = FALSE)
}
#-------------------------------------------------------
#Viimeisten ruutujen etsintä
#apu2 <- NULL
#apu2 <- LahimmatRuudut(LO,LA, apuPD, apuPM, Piste)
#apuPM <- apuPM[apu2,]
#apuPD <- apuPD[apu2,]
#apuHI[1] <- sum( apuPM[,(1 + 2*Piste)] * apuPD[,3] )
#apuLO[1] <- sum( apuPM[,(2 + 2*Piste)] * apuPD[,3] )
#apuSumma[1] <- apuHI[1]
#apuSumma[2] <- apuLO[1]
#------------------------------------------------------------------
ylaraja <- length(apuHI)
for(i in 1: ylaraja){
apuSumma[(2*i - 1)] <- apuHI[(ylaraja - i + 1)]
apuSumma[2*i] <- apuLO[(ylaraja - i + 1)]
}
SummaPixCi <- rbind(SummaPixCi, t(apuSumma))
}
# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
BR <- 20/(24 * 60 * 60)
# Q <- 3 * 1/(365 * 24 * 60 * 60)kg/s
Q <- 3/(365 * 24 * 60 * 60)
IF <- SummaPixCi * BR/Q
TalletusKohde <- paste(talletuspolku,"IF_Piste_",Piste,"_Distance.txt",sep="")
write.table(IF, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
}
#----------------------------------------------------------
#Kolmion pinta-ala, kun tiedetään sivujen pituudet
HeronAla <- function(a,b,c){
p <- (a+b+c)/2
sqrt(p*(p-a)*(p-b)*(p-c))
}
KantaKorkeus <- function(k, h){
1/2 * k * h
}
#--------------------------------------------------------------------------
# Voi tulla ongelmia, jos piste on kulmapisteiden välisellä janalla !
# Jos palautetun vektorin 1. koordinaatti on 1, niin tarkasteltavan
# kolmion tilavuus on nolla. Muutoin ensimmäinen koordinaatti on 1, 2. koord.
# on A pisteen x-koord, 3. koord. on A pisteen y-koord., ..., 8. koord. on 1 jos
# piste on kolmion sisällä, muutoin 8. koord on nolla.
KolmionSisalla <- function(x1,x2,D1,D2,E1,E2,F1,F2){
xvek <- c(D1,E1,F1)
yvek <- c(D2,E2,F2)
jarjvek <- NULL
paha <- 0
#Huonot tapaukset
#pisteiden D ja E kautta kulkeva suora -> testataan ovatko pisteet samalla suoralla
if(xvek[1] != xvek[2]){
k <- (yvek[2] - yvek[1])/(xvek[2] - xvek[1])
y <- k*(xvek[3] - xvek[1]) + yvek[1]
if(yvek[3] == y){ paha <- 1}
}
# Kaksi pistettä samat tai samalla vaakasuoralla tai samalla pystysuoralla
if( (xvek[1] == xvek[2] & yvek[1] == yvek[2]) |
(xvek[1] == xvek[3] & yvek[1] == yvek[3]) |
(xvek[2] == xvek[3] & yvek[2] == yvek[3]) |
(xvek[1] == xvek[2] & xvek[1] == xvek[3]) |
(yvek[1] == yvek[2] & yvek[1] == yvek[3])
)
{
paha <- 1
}
jarjvek[1] <- paha
if(paha == 0)
{
apu1 <- min(xvek)
if( length(xvek[xvek==apu1])==1)
{
jarjvek[2] <- xvek[xvek == apu1]
jarjvek[3] <- yvek[xvek == apu1]
apu2 <- max(yvek[xvek != apu1])
if( length(yvek[(xvek!= apu1 & yvek == apu2)]) ==1)
{
jarjvek[4] <- xvek[(xvek!= apu1 & yvek == apu2)]
jarjvek[5] <- yvek[(xvek!= apu1 & yvek == apu2)]
jarjvek[6] <- xvek[(xvek!= apu1 & yvek != apu2)]
jarjvek[7] <- yvek[(xvek!= apu1 & yvek != apu2)]
}
if( length(yvek[(xvek!= apu1 & yvek == apu2)]) ==2)
{
apu3 <- min(xvek[xvek!= apu1])
jarjvek[4] <- xvek[(xvek!= apu1 & xvek == apu3)]
jarjvek[5] <- yvek[(xvek!= apu1 & xvek == apu3)]
jarjvek[6] <- xvek[(xvek!= apu1 & xvek != apu3)]
jarjvek[7] <- yvek[(xvek!= apu1 & xvek != apu3)]
}
#Jos piste C on janan AB yläpuolella, niin vaihdetaan B ja C
janaT <- NULL
janaT[1] <- (jarjvek[6] - jarjvek[2])/(jarjvek[4] - jarjvek[2])
janaY <- NULL
janaY[1] <- jarjvek[5] * janaT[1] + (1 - janaT[1])* jarjvek[3]
if(jarjvek[7] > janaY[1]){
apu4 <- jarjvek
jarjvek[4:5] <- apu4[6:7]
jarjvek[6:7] <- apu4[4:5]
}
}
if( length(xvek[xvek==apu1])==2)
{
apu2 <- min(yvek[xvek == apu1])
jarjvek[2] <- xvek[(xvek == apu1 & yvek == apu2)]
jarjvek[3] <- yvek[(xvek == apu1 & yvek == apu2)]
jarjvek[4] <- xvek[(xvek == apu1 & yvek != apu2)]
jarjvek[5] <- yvek[(xvek == apu1 & yvek != apu2)]
jarjvek[6] <- xvek[xvek != apu1 ]
jarjvek[7] <- yvek[xvek != apu1 ]
}
}
if(jarjvek[1]==0)
{
#Janan koordinaatit aina järjestyksessä AB, BC, AC
janaT <- NULL
janaT[1] <- (x1 - jarjvek[2])/(jarjvek[4] - jarjvek[2])
janaT[2] <- (x1 - jarjvek[4])/(jarjvek[6] - jarjvek[4])
janaT[3] <- (x1 - jarjvek[2])/(jarjvek[6] - jarjvek[2])
janaY <- NULL
janaY[1] <- jarjvek[5] * janaT[1] + (1 - janaT[1])* jarjvek[3]
janaY[2] <- jarjvek[7] * janaT[2] + (1 - janaT[2])* jarjvek[5]
janaY[3] <- jarjvek[7] * janaT[3] + (1 - janaT[3])* jarjvek[3]
if(jarjvek[2]==jarjvek[4])
{
jarjvek[8] <- 0
if((jarjvek[2]<= x1) & (x2 <= janaY[2]) & (janaY[3]<= x2))
{
jarjvek[8] <- 1
}
}
if(jarjvek[2]!=jarjvek[4])
{
if(jarjvek[4]==jarjvek[6])
{
jarjvek[8] <- 0
if((x1 <= jarjvek[4]) & (x2 <= janaY[1]) & (janaY[3]<= x2))
{
jarjvek[8] <- 1
}
}
if(jarjvek[4]!=jarjvek[6])
{
jarjvek[8] <- 0
if( ((jarjvek[6] > jarjvek[4]) & (x2 <= janaY[1]) &
(x2 <= janaY[2]) & (janaY[3] <= x2)) |
((jarjvek[6] < jarjvek[4]) & (x2 <= janaY[1]) &
(janaY[2] <= x2) & (janaY[3] <= x2))
)
{
jarjvek[8] <- 1
}
}
}
}
jarjvek
}
#------------------------------------------------------------------------
#Pistevektoriin syötetään KolmionSisalla tuottama vektori!!!!!
# x1, x2 on ympyrän keskipiste
PisteenAlue <- function(x1,x2, pistevek){
alue <- c(0,0,0)
# C Kärkenä piste C
A <- pistevek[1:2]
B <- pistevek[3:4]
C <- pistevek[5:6]
t1 <-
((x2 - C[2])*(A[1] - C[1]) - (x1 - C[1])*(A[2] - C[2]))/
((x1 - C[1])*(B[2] - A[2]) - (x2 - C[2])*(B[1] - A[1]))
if(t1==0){
suunta <- -1
if((A[1] - C[1])!=0 & ((x1 - C[1])/(A[1] - C[1]))>=0)
suunta <- 1
if((A[1] - C[1])==0 & ((x2 - C[2])/(A[2] - C[2]))>=0)
suunta <- 1
t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, A[1],A[2], C[1],C[2],2)* suunta
}
if(t1==1){
suunta <- -1
if((B[1] - C[1])!=0 & ((x1 - C[1])/(B[1] - C[1]))>=0)
suunta <- 1
if((B[1] - C[1])==0 & ((x2 - C[2])/(B[2] - C[2]))>=0)
suunta <- 1
t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, B[1],B[2], C[1],C[2],2)*suunta
}
if(t1!=0 & t1!=1)
ifelse(x1==C[1],
{
suunta <- -1
if(sum(C - t1*B + (1 - t1)*A)/sum(C - c(x1,x2))>=0)
suunta <- 1
t2 <- Etaisyys( 6371, x1,x2,C[1],C[2],2)/
Etaisyys(6371,C[1],C[2],t1*B[1] + (1 - t1)*A[1],t1*B[2] + (1 - t1)*A[2],2) *suunta
},
t2 <-
(x1 - C[1])/(t1*B[1] + (1 - t1)*A[1] - C[1])
)
if(0<= t1 & t1 <= 1 & 1 < t2)
{
alue[1] <- 1
alue[2] <- t1
alue[3] <- t2
}
if(0< t1 & t1 < 1 & t2 < 0)
{
alue[1] <- 4
alue[2] <- t1
alue[3] <- t2
}
if(alue[1] == 0){
# C Kärkenä piste A
#---------------------------
A <- pistevek[3:4]
B <- pistevek[5:6]
C <- pistevek[1:2]
t1 <-
((x2 - C[2])*(A[1] - C[1]) - (x1 - C[1])*(A[2] - C[2]))/
((x1 - C[1])*(B[2] - A[2]) - (x2 - C[2])*(B[1] - A[1]))
if(t1==0){
suunta <- -1
if((A[1] - C[1])!=0 & ((x1 - C[1])/(A[1] - C[1]))>=0)
suunta <- 1
if((A[1] - C[1])==0 & ((x2 - C[2])/(A[2] - C[2]))>=0)
suunta <- 1
t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, A[1],A[2], C[1],C[2],2)* suunta
}
if(t1==1){
suunta <- -1
if((B[1] - C[1])!=0 & ((x1 - C[1])/(B[1] - C[1]))>=0)
suunta <- 1
if((B[1] - C[1])==0 & ((x2 - C[2])/(B[2] - C[2]))>=0)
suunta <- 1
t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, B[1],B[2], C[1],C[2],2)*suunta
}
if(t1!=0 & t1!=1)
ifelse(x1==C[1],
{
suunta <- -1
if(sum(C - t1*B + (1 - t1)*A)/sum(C - c(x1,x2))>=0)
suunta <- 1
t2 <- Etaisyys( 6371, x1,x2,C[1],C[2],2)/
Etaisyys(6371,C[1],C[2],t1*B[1] + (1 - t1)*A[1],t1*B[2] + (1 - t1)*A[2],2) *suunta
},
t2 <-
(x1 - C[1])/(t1*B[1] + (1 - t1)*A[1] - C[1])
)
#--------------------------------------
if(0<= t1 & t1 <= 1 & 1 < t2)
{
alue[1] <- 3
alue[2] <- t1
alue[3] <- t2
}
if(0< t1 & t1 < 1 & t2 < 0)
{
alue[1] <- 6
alue[2] <- t1
alue[3] <- t2
}
if(alue[1] == 0){
# C Kärkenä piste B
#-----------------------------------
A <- pistevek[5:6]
B <- pistevek[1:2]
C <- pistevek[3:4]
t1 <-
((x2 - C[2])*(A[1] - C[1]) - (x1 - C[1])*(A[2] - C[2]))/
((x1 - C[1])*(B[2] - A[2]) - (x2 - C[2])*(B[1] - A[1]))
if(t1==0){
suunta <- -1
if((A[1] - C[1])!=0 & ((x1 - C[1])/(A[1] - C[1]))>=0)
suunta <- 1
if((A[1] - C[1])==0 & ((x2 - C[2])/(A[2] - C[2]))>=0)
suunta <- 1
t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, A[1],A[2], C[1],C[2],2)* suunta
}
if(t1==1){
suunta <- -1
if((B[1] - C[1])!=0 & ((x1 - C[1])/(B[1] - C[1]))>=0)
suunta <- 1
if((B[1] - C[1])==0 & ((x2 - C[2])/(B[2] - C[2]))>=0)
suunta <- 1
t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, B[1],B[2], C[1],C[2],2)*suunta
}
if(t1!=0 & t1!=1)
ifelse(x1==C[1],
{
suunta <- -1
if(sum(C - t1*B + (1 - t1)*A)/sum(C - c(x1,x2))>=0)
suunta <- 1
t2 <- Etaisyys( 6371, x1,x2,C[1],C[2],2)/
Etaisyys(6371,C[1],C[2],t1*B[1] + (1 - t1)*A[1],t1*B[2] + (1 - t1)*A[2],2) *suunta
},
t2 <-
(x1 - C[1])/(t1*B[1] + (1 - t1)*A[1] - C[1])
)
#---------------------------------------------------
if(0<= t1 & t1 <= 1 & 1 < t2)
{
alue[1] <- 5
alue[2] <- t1
alue[3] <- t2
}
if(0< t1 & t1 < 1 & t2 < 0)
{
alue[1] <- 2
alue[2] <- t1
alue[3] <- t2
}
}
}
alue
}
#-----------------------------------------------------------------------------
#Laskee ympyrän ja kolmion leikkauksen alan, kun ympyrän keskipiste on yksi
#kolmion kulmista
#Palauttaa vektorin jonka 1. komponentti on leikkausen pinta-ala
# 2. komp. kolmion pinta-ala, 3. on näiden suhde ja 4. kolmion korkeus.
AlaYmpyranKpKolmionKulma <- function(kp1, kp2, r, kv1,kv2,kv3,kv4)
{
ala <- 0
vek <- c(kp1, kp2, r, kv1,kv2,kv3,kv4)
tilavek <- KolmionSisalla(vek[1],vek[2],vek[4],vek[5],vek[6], vek[7],vek[1],vek[2])
if(tilavek[1]==1 ){
tulos <- c(0,0,0,0)
}
if(tilavek[1]==0){
hNollaraja <- 10^(-12)
kp <- c(kp1,kp2)
A <- tilavek[2:3]
B <- tilavek[4:5]
C <- tilavek[6:7]
if(kp[1]==A[1] & kp[2]==A[2]){
A <- B
B <- C
}
if(kp[1]==B[1] & kp[2]==B[2]){
B <- A
A <- C
}
#oikea kolmioa
if(tilavek[1]==0)
{
AB <- Etaisyys(6371,A[1],A[2],B[1],B[2],1)
kpB <- Etaisyys(6371,kp[1],kp[2],B[1],B[2],1)
kpA <- Etaisyys(6371,kp[1],kp[2],A[1],A[2],1)
cosAlpha <- (kpA^2 + AB^2 - kpB^2)/(2*AB*kpA)
cosBeta <- (kpB^2 + AB^2 - kpA^2)/(2*AB*kpB)
cosGamma <- (kpA^2 + kpB^2 - AB^2)/(2*kpB*kpA)
if( cosAlpha >= 0)
{
h <- sqrt(kpA^2 - kpA^2*cosAlpha^2)
#A ja B eivät kuulu ympyrään ja ympyrä ei leikkaa niiden välistä janaa.
if(r <= h){
ala <- 1/2 * r^2 * acos(cosGamma)
}
# A ja B eivät kuulu ympyrään ja ympyrä leikkaa niiden välistä janaa.
if(h < r & kpA > r & kpB > r){
if(cosBeta <=0) ala <- 1/2*r^2*acos(cosGamma)
if(cosBeta > 0)
ala <- 1/2 * r^2*(acos(h/kpB) - acos(h/r)) +
KantaKorkeus(2*sqrt(r^2 - h^2),h) +
1/2 * r^2* (acos(h/kpA) - acos(h/r))
}
# B kuuluu ympyrään ja A ei kuulu
if(h < r & kpB <= r & kpA > r){
if(cosBeta > 0)
x <- sqrt(abs(kpB^2 - h^2)) + sqrt(r^2 - h^2)
if(cosBeta <= 0)
x <- sqrt(r^2 - h^2) - sqrt(abs(kpB^2 - h^2))
ala <- KantaKorkeus(x,h) + 1/2 * r^2 * (acos(h/kpA) - acos(h/r))
}
# B ei kuulu ympyrään ja A kuuluu
if(h < r & kpB > r & kpA <= r){
x <- sqrt(abs(kpA^2 - h^2)) + sqrt(r^2 - h^2)
ala <- KantaKorkeus(x,h) + 1/2 * r^2 * (acos(h/kpB) - acos(h/r))
}
# A ja B kuuluvat ympyrään
if(h < r & kpA <= r & kpB <= r){
ala <- KantaKorkeus(AB, h)
}
tulos <- c(ala, (1/2*h*AB), (ala/(1/2*h*AB)),h)
}
if( cosAlpha < 0)
{
h <- sqrt(kpB^2 - kpB^2*cosBeta^2)
# A ja B eivät kuulu ympyrään ja ympyrä ei leikkaa niiden välistä janaa.
if(r <= h){
ala <- 1/2 * r^2 * acos(cosGamma)
}
# A ja B eivät kuulu ympyrään ja ympyrä leikkaa niiden välistä janaa.
# Tämä ei ole mahdollista, koska Alpha > 90 astetta
# B kuuluu ympyrään ja A ei kuulu
# Ei mahdollista, koska kpA < kpB
# B ei kuulu ympyrään ja A kuuluu
if(h < r & kpB > r & kpA <= r){
x <- sqrt(r^2 - h^2) - sqrt(abs(kpA^2 - h^2))
ala <- KantaKorkeus(x,h) + 1/2 * r^2 * (acos(h/kpB) - acos(h/r))
}
# A ja B kuuluvat ympyrään
if(h < r & kpA <= r & kpB <= r){
ala <- KantaKorkeus(AB, h)
}
tulos <- c(ala, (1/2*h*AB), (ala/(1/2*h*AB)),h)
}
}
}
tulos
}
#-----------------------------------------------------------------------------
#Laskee ympyrän ja kolmion leikkauksen alan, kun ympyrän keskipiste on kolmion sisällä
AlaYmpyranKpKolmionSisalla <- function(kp1, kp2, r, kv1,kv2,kv3,kv4,kv5,kv6)
{
tilavek <- KolmionSisalla(kp1,kp2,kv1,kv2,kv3,kv4,kv5,kv6)
hNollaraja <- 10^(-6)
kp <- c(kp1,kp2)
A <- tilavek[2:3]
B <- tilavek[4:5]
C <- tilavek[6:7]
KolmionKulmat <- rbind(A,B,C,A)
#oikea kolmioa ja piste kolmion sisalla
if(tilavek[1]==0 & tilavek[8]==1)
{
AB <- Etaisyys(6371,tilavek[2],tilavek[3],tilavek[4],tilavek[5],1)
BC <- Etaisyys(6371,tilavek[4],tilavek[5],tilavek[6],tilavek[7],1)
AC <- Etaisyys(6371,tilavek[2],tilavek[3],tilavek[6],tilavek[7],1)
etA <- Etaisyys(6371, kp1, kp2, tilavek[2], tilavek[3],1)
etB <- Etaisyys(6371, kp1, kp2, tilavek[4], tilavek[5],1)
etC <- Etaisyys(6371, kp1, kp2, tilavek[6], tilavek[7],1)
etv <- c(etA, etB, etC, etA)
valit <- c(AB, BC, AC)
#h <- c(-100,-100,-100)
h <- c(1,1,1)
# Tapaukset, joissa ymp. kp on kolmion kulmapiste
if( (kp == A)[1] & (kp == A)[2]){ h[1] <- 0
h[3] <- 0
}
if( (kp == B)[1] & (kp == B)[2]){ h[1] <- 0
h[2] <- 0
}
if( (kp == C)[1] & (kp == C)[2]){ h[2] <- 0
h[3] <- 0
}
# Tapaukset, joissa ymp. kp on kolmion kulmapisteiden välisellä suoralla
if(A[1]!=B[1])
{
k <- (B[2] - A[2])/(B[1] - A[1])
y <- k* (kp[1] - A[1]) + A[2]
if( abs(kp[2] - y) < hNollaraja) h[1] <- 0
}
if( (A[1] == B[1] & kp[1]==A[1]))
{
h[1] <- 0
}
if(B[1]!=C[1])
{
k <- (C[2] - B[2])/(C[1] - B[1])
y <- k* (kp[1] - B[1]) + B[2]
if( abs(kp[2] - y) < hNollaraja) h[2] <- 0
}
if( (B[1] == C[1] & kp[1]==B[1]))
{
h[2] <- 0
}
if(C[1]!=A[1])
{
k <- (A[2] - C[2])/(A[1] - C[1])
y <- k* (kp[1] - C[1]) + C[2]
if( abs(kp[2] - y) < hNollaraja) h[3] <- 0
}
if( (C[1] == A[1] & kp[1]==C[1]))
{
h[3] <- 0
}
ala <- NULL
for(i in 1:3){
if(h[i] >= hNollaraja){
ala[i] <- AlaYmpyranKpKolmionKulma(kp[1], kp[2], r, KolmionKulmat[i,1],
KolmionKulmat[i,2],KolmionKulmat[(i+1),1],KolmionKulmat[(i+1),2])[1]
}
if(h[i] < hNollaraja){
ala[i] <- 0
}
}
tulos <- NULL
tulos[1] <- sum(ala)
tulos[2] <- KantaKorkeus(valit[1],
sqrt(valit[3]^2 - ((valit[1]^2 + valit[3]^2 - valit[2]^2)/(2*valit[1]))^2))
tulos[3] <- tulos[1]/tulos[2]
}
if(tilavek[1]==1 & length(tilavek)==1)
{
tulos <- c(0,0,0)
}
if(length(tilavek) > 1 & tilavek[8]== 0)
{
tulos <- c(-1,-1,-1)
}
tulos
}
#-------------------------------------------------------------------------------------------
#Pisteiden koordinaatit syötettävä kiertosuunnassa vastapäivään!!!!
YmpyraNelikulmioKpSisalla <- function(kpLO, kpLA, r, ALO, ALA, BLO, BLA, CLO, CLA, DLO, DLA)
{
tyovek <- c(kpLO, kpLA, r, ALO, ALA, BLO, BLA, CLO, CLA, DLO, DLA)
kp <- tyovek[1:2]
KolmionKulmat <- rbind(c(tyovek[4:5],tyovek[6:7]),
c(tyovek[6:7],tyovek[8:9]),
c(tyovek[8:9],tyovek[10:11]),
c(tyovek[10:11],tyovek[4:5]))
ala <- NULL
for(i in 1:4){
ala[i] <- AlaYmpyranKpKolmionKulma(kp[1], kp[2], r,
KolmionKulmat[i,1], KolmionKulmat[i,2],
KolmionKulmat[i,3],KolmionKulmat[i,4])[1]
}
AB <- Etaisyys(6371, tyovek[4],tyovek[5], tyovek[6], tyovek[7],1)
BC <- Etaisyys(6371, tyovek[6], tyovek[7], tyovek[8], tyovek[9],1)
CD <- Etaisyys(6371, tyovek[8], tyovek[9], tyovek[10], tyovek[11],1)
DA <- Etaisyys(6371, tyovek[10], tyovek[11], tyovek[4], tyovek[5],1)
AC <- Etaisyys(6371, tyovek[4], tyovek[5], tyovek[8], tyovek[9],1)
tulos <- NULL
tulos[1] <- sum(ala)
tulos[2] <- HeronAla(AB, BC, AC) + HeronAla(AC, CD, DA)
tulos[3] <- tulos[1]/tulos[2]
tulos
}
#----------------------------------------------------------------------------
#kv on vektori, joka sisaltaa kolmion kulmat
KolmioYmpyranAla <- function(kp1, kp2, r, kv1,kv2,kv3,kv4,kv5,kv6)
{
ala <- NULL
tilavek <- KolmionSisalla(kp1,kp2,kv1,kv2,kv3,kv4,kv5,kv6)
#oikea kolmioa ja piste kolmion sisalla
if(tilavek[1]==0 & tilavek[8]==1)
{
ala <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r, kv1,kv2,kv3,kv4,kv5,kv6)
}
#Oikea kolmio ja piste kolmion ulkopuolella
if(tilavek[1]==0 & tilavek[8]==0)
{
paikkavek <- PisteenAlue(kp1,kp2,tilavek[2:7])
A <- tilavek[2:3]
B <- tilavek[4:5]
C <- tilavek[6:7]
AB <- Etaisyys(6371,A[1],A[2],B[1],B[2],1)
BC <- Etaisyys(6371,B[1],B[2],C[1],C[2],1)
AC <- Etaisyys(6371,A[1],A[2],C[1],C[2],1)
if(paikkavek[1] == 1)
{
# kp = (1 - paikkavek[3])*C + paikkavek[3]*(paikkavek[2]*B + (1-paikkavek[2])*A)
uusiA <- (1 - (paikkavek[3] + 1))*C + (paikkavek[3] + 1)*A
uusiB <- (1 - (paikkavek[3] + 1))*C + (paikkavek[3] + 1)*B
uusiC <- C
#Pisteiden koordinaatit syötettävä kiertosuunnassa vastapäivään!!!!!!!!!!
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r,
uusiB[1], uusiB[2], uusiA[1], uusiA[2], uusiC[1], uusiC[2])[1]-
YmpyraNelikulmioKpSisalla(kp1, kp2, r,
uusiB[1], uusiB[2], uusiA[1], uusiA[2], A[1], A[2], B[1], B[2])[1]
}
if(paikkavek[1] == 2)
{
#iso kolmio - pikkukolmiot
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], C[1],C[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, B[1], B[2], C[1],C[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], B[1],B[2], kp1, kp2)[1]
}
if(paikkavek[1] == 3)
{
uusiA <- A
uusiB <- (1 - (paikkavek[3] + 1))*A + (paikkavek[3] + 1)*B
uusiC <- (1 - (paikkavek[3] + 1))*A + (paikkavek[3] + 1)*C
#Pisteiden koordinaatit syötettävä kiertosuunnassa vastapäivään!!!!!!!!!!
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r,
uusiB[1], uusiB[2], uusiA[1], uusiA[2], uusiC[1], uusiC[2])[1]-
YmpyraNelikulmioKpSisalla(kp1, kp2, r,
uusiB[1], uusiB[2], B[1], B[2], C[1], C[2], uusiC[1], uusiC[2])[1]
}
if(paikkavek[1] == 4)
{
#iso kolmio - pikkukolmiot
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], B[1],B[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, B[1], B[2], C[1],C[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], C[1],C[2], kp1, kp2)[1]
}
if(paikkavek[1] == 5)
{
uusiA <- (1 - (paikkavek[3] + 1))*B + (paikkavek[3] + 1)*A
uusiB <- B
uusiC <- (1 - (paikkavek[3] + 1))*B + (paikkavek[3] + 1)*C
#Pisteiden koordinaatit syötettävä kiertosuunnassa vastapäivään!!!!!!!!!!
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r,
uusiB[1], uusiB[2], uusiA[1], uusiA[2], uusiC[1], uusiC[2])[1]-
YmpyraNelikulmioKpSisalla(kp1, kp2, r,
uusiA[1], uusiA[2], uusiC[1], uusiC[2], C[1], C[2], A[1], A[2])[1]
}
if(paikkavek[1] == 6)
{
#iso kolmio - pikkukolmiot
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r, C[1], C[2], B[1],B[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, B[1], B[2], A[1],A[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], C[1],C[2], kp1, kp2)[1]
}
ala[2] <- HeronAla(AB,BC,AC)
ala[3] <- ala[1]/ala[2]
}
ala
}
#----------------------------------------------------------------------------
# Syötä kulmat kiertojärjestyksessä vastapäivään!!!!!!
YmpyraNelikulmioAla <-function(kp,r,kulmat)
{
tulos <- NULL
apu1 <- KolmioYmpyranAla(kp[1], kp[2], r, kulmat[1],kulmat[2],kulmat[3],kulmat[4],kulmat[7],kulmat[8])
apu2 <- KolmioYmpyranAla(kp[1], kp[2], r, kulmat[5],kulmat[6],kulmat[3],kulmat[4],kulmat[7],kulmat[8])
tulos[1:2] <- apu1[1:2] + apu2[1:2]
tulos[3] <- tulos[1]/tulos[2]
tulos
}
pisteLO <- kojetaulu[4,1]
pisteLA <- kojetaulu[4,2]
leik <- NULL
for(sade in 0:30){
leik[sade] <- YmpyraNelikulmioAla(c(25,60),sade, c(pisteLO,pisteLA,pisteLO+LOLisays,pisteLA,
pisteLO+LOLisays,pisteLA + LALisays,pisteLO,pisteLA+ LALisays))[3]
}
KolmioYmpyranAla(25,61,1,24,59,25,61,26,59)
for(km in seq(0,20,by=0.2)){
cat(km)
cat(" ")
cat(KolmioYmpyranAla(1,1,0.5*km,5,0,2,3,-2,0))
cat("\n")
}
kp1 <- -5
kp2 <- -1500
r <- 1503
tilavek <- KolmionSisalla(kp1,kp2,10,1,-4,1,-4,8)
tilavek
PisteenAlue(kp1,kp2,tilavek[2:7])
KolmioYmpyranAla(kp1,kp2,r,10,1,-4,1,-4,8)
kp1 <- 5
kp2 <- -4
#Ongelma säteellä 9!!!!!!!!!!!!!!!!!!!!!
r <- 9
KolmioYmpyranAla(kp1,kp2,9,5,2,12,2,5,6)
tilavek <- KolmionSisalla(kp1,kp2,5,2,12,2,5,6)
paikkavek <- PisteenAlue(kp1,kp2,tilavek[2:7])
x1 <- kp1
x2 <- kp2
pistevek <- tilavek[2:7]
kp <- c(1,500)
kp <- c(-4,4)
PisteenAlue(kp[1],kp[2],KolmionSisalla(kp[1],kp[2],A[1],A[2],B[1],B[2],C[1],C[2])[2:7])
x1 <- kp[1]
x2 <- kp[2]
pistevek <- KolmionSisalla(kp[1],kp[2],A[1],A[2],B[1],B[2],C[1],C[2])[2:7]
|
datansiirtoYMS.r (24.2.2006)
#-----------------------------------------------------------------------
#Maan pinnalla sijaitsevien pisteiden etäisyyden laskeva funktio
#-----------------------------------------------------------------------
#maan säde n. 6371 km
#jos laskutapa on 1, niin lasketaan 2 pisteen välimatka pallon pintaa pitkin
# jos laskutapa =2, niin lasketaan pisteiden (p1LO, p1LA) ja (p2LO, p2LA) euklidinen etäisyys,
# tällöin säteellä ei ole merkitystä
Etaisyys <- function(sade, p1LO, p1LA, p2LO, p2LA, laskutapa)
{
if(laskutapa == 1){
piste1LA <- p1LA * (2*pi/360)
piste1LO <- p1LO * (2*pi/360)
piste2LA <- p2LA * (2*pi/360)
piste2LO <- p2LO * (2*pi/360)
#pallokoorinaatistoa käyttäen
SuoraValiMatka <- sade * sqrt( (cos(piste1LO)*cos(piste1LA) - cos(piste2LO)*cos(piste2LA))^2
+ (sin(piste1LO)*cos(piste1LA) - sin(piste2LO)*cos(piste2LA))^2
+ (sin(piste1LA) - sin(piste2LA))^2 )
#cosini lause
Kulma <- acos(1 - SuoraValiMatka^2/(2*sade^2))
#kaaren pituus on kulma(radiaaneina)*säde
tulos <- sade * Kulma
}
if(laskutapa == 2)
{
tulos <- sqrt((p1LO - p2LO)^2 + (p1LA - p2LA)^2)
}
tulos
}
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
DataRuudut <- 10
# HUOM! tiedostoissa pääte .out on 50km ja .out.1 on 10km
if (DataRuudut == 10) paate <- ".out.1"
if (DataRuudut == 50) paate <- ".out"
PopulaatioData <- NULL
PopulaatioData <-
read.table(paste(polku,"Population_Europe\\",DataRuudut,"km\\ciesin_center_63lon_",DataRuudut,"km_2",sep=""))
#----------------------------------------------------------------------------
# PopulaatioData siirtävä funktio
#-------------------------------------------------------------------------
LohkoPituus <- 354
Lohkoja <- dim(PopulaatioData)[1]/LohkoPituus
# sPD <- PDSiirto(0,0,2,0,PopulaatioData, Lohkoja, LohkoPituus)
PDSiirto <- function(poh, ita, ete, lan, PopulaatioData, Lohkoja, LohkoPituus){
UusiPD <- NULL
UusiPD <- PopulaatioData
#Siirto etelään
if(ete > 0){
#Siirretään dataa
for( i in 1:(Lohkoja-ete)){
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <-
PopulaatioData[((i + ete - 1)*LohkoPituus + 1):((i+ete)*LohkoPituus),3]
}
#Täytetään lopusta nollilla
for(i in 1:ete){
UusiPD[((Lohkoja-ete)*LohkoPituus+1):((Lohkoja - ete +1)*LohkoPituus),3] <- t(t(rep(0,LohkoPituus)))
}
}
if(lan > 0){
#Siirto länteen
#Siirretään dataa
for( i in 1:Lohkoja){
apu1 <- NULL
apu1 <- PopulaatioData[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3]
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <-
t(t(c(apu1[(1+lan):LohkoPituus],rep(0,lan))))
}
}
if(poh > 0){
#Siirto pohjoiseen
#Siirretään dataa
for( i in (poh + 1):(Lohkoja)){
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <-
PopulaatioData[((i-poh)*LohkoPituus + 1):((i-poh+1)*LohkoPituus),3]
}
#Täytetään lopusta nollilla
for(i in 1:poh){
UusiPD[((i-1)*LohkoPituus+1):((i)*LohkoPituus),3] <- t(t(rep(0,LohkoPituus)))
}
}
if(ita > 0){
#Siirto itään
#Siirretään dataa
for( i in 1:Lohkoja){
apu1 <- NULL
apu1 <- PopulaatioData[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3]
UusiPD[((i-1)*LohkoPituus + 1):(i*LohkoPituus),3] <-
t(t(c(rep(0,ita),apu1[1:(LohkoPituus - ita)])))
}
}
UusiPD
}
kojetaulu2 <- LahimmatRuudut(1,25,60.08,sPD,60,1, 58, 65)
kojetaulu3 <- LahimmatRuudut(1,25,60,sPD,500,1, 52.2,64.8)
#Pietari lähellä pistettä 30.2, 60
PietariTot <- LahimmatRuudut(2.5,30.1,60,PopulaatioData,PopulaatioData,100,1, 0, 80, 1)
PData <- LahimmatRuudut(25,60.08,PopulaatioData,1, 58, 64.8)
apuri <- PopulaatioData[58 <= PopulaatioData[,2] & PopulaatioData[,2] <= 65, ]
apuri[PData,]
PData
#-----------------------------------------------------------
#Lähimpien Ruutujen etsintä
#-----------------------------------------------------------
# Annetaan pisteen Lo ja La koordinaatit, piste itse, Populaatio data ja Pienhiukkasdata.
# Viimeisin annetaan, jotta ruutujen määrää voidaan vähentää.
# Voidaan myös rajata antamalla alaraja ja yläraja Latitudelle
#-----------------------------------------------------------
LahimmatRuudut <-function(Piste,LO,LA, PD, PM, Sade,Sateella, LAalaraja, LAylaraja, totuustaulu){
LOLisays <- 0.1981
LOLisaysArvio <- 0.199
LALisays <- 0.08995
LALisaysArvio <- 0.09
apuPD <- NULL
apuPD <- PD
#apuPD <- PD[(PD[,3]!=0 & LAalaraja <= PD[,2] & PD[,2] <= LAylaraja),]
#apuPD <- PD[(PD[,3]!=0 & (PM[,]!=0 | PM[,]!=0)),]
#apuPD <- PD[PD[,3]!=0 & (PM[,(1 + 2*Piste)]!=0 | PM[,(2 + 2*Piste)]!=0),]
#apuPD <- PD[PD[,3]!=0 & LAalaraja <= PD[,2] & PD[,2] <= LAylaraja,]
#apuPD <- PD[LAalaraja <= PD[,2] & PD[,2] <= LAylaraja,]
apu2 <- NULL
if(Sateella==1){
for(k in 1:dim(apuPD)[1]){
if( (apuPD[k,1] >= LO & apuPD[k,2] >= LA) | (apuPD[k,1] < LO & apuPD[k,2] < LA ) ){
apu2[k] <-
(Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2]),1) <= Sade &
Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2] + LALisays),1) <= Sade)
}
if((apuPD[k,1] < LO & apuPD[k,2] >= LA) | (apuPD[k,1] >= LO & apuPD[k,2] < LA ) ){
apu2[k] <-
(Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2]),1) <= Sade &
Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2] + LALisays),1) <= Sade)
}
}
}
if(Sateella==0)
{
apu2 <- rep(FALSE,(dim(apuPD)[1]))
for(k in 1:dim(apuPD)[1]){
if( (apuPD[k,1] <= LO & apuPD[k,2] <= LA) &
(LO < (apuPD[k,1] + LOLisays) & LA < (apuPD[k,2] + LALisays)) )
{
apu2[k] <- (1==1)
Ruutu <- k
}
}
for(k in 1:dim(apuPD)[1]){
if
(
((apuPD[Ruutu,1]-LOLisaysArvio)<= apuPD[k,1]) & (apuPD[k,1] <= (apuPD[Ruutu,1]+LOLisaysArvio)) &
((apuPD[Ruutu,2]-LALisaysArvio) <= apuPD[k,2]) & (apuPD[k,2] <= (apuPD[Ruutu,2]+LALisaysArvio))
)
{
apu2[k] <- (1==1)
if(length(apu2[apu2]) == 9){ k <- dim(apuPD)[1] + 1}
}
}
}
if(totuustaulu == 1)
{
tulos <- apu2
}
if(totuustaulu == 0)
{
tulos <- apuPD[apu2,]
}
tulos
}
#------------------------------------------------------
# IF-datan piirtävä koodi
#------------------------------------------------
#polku <- "c:\\omat\\dataa\\concentration_point\\"
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
SadeLisays <- 1
PieninSade <- 0
MaxSade <- 100
PiirYlaRaja <- 3.5*10^(-6)
for(i in 4:4){
Piste <- i
data <- paste(polku,"IF_Piste_",Piste,"_Distance_Osuus_Kokeilu_SadeVali_1_Suurin_100.txt",sep="")
#data <- paste(polku,"IF_Piste_",Piste,"_Distance.txt",sep="")
#data <- paste(polku,"IF_Piste_",Piste,"_Distance_alle_50.txt",sep="")
#data <- paste(polku,"IF_Piste_",Piste,"_Distance_Siirretty.txt",sep="")
#kohde <- paste(talletuspolku, "Piste_",Piste,".ps",sep="")
#kohde <- paste(talletuspolku, "Piste_",Piste,"_Siirretty.ps",sep="")
#kohde <- paste(talletuspolku, "Piste_",Piste,"_Siirretty_Osuus_Kokeilu_SadeVali_1_Suurin_100.ps",sep="")
#kohde <- paste(talletuspolku, "Piste_",Piste,"_alle_50.ps",sep="")
kohde <- paste(talletuspolku, "Piste_",Piste,"_Siirretty_Osuus_Kokeilu_SadeVali_1_Suurin_100_Helsingin_Skaala.ps",sep="")
postscript(file=kohde, width=10,height=10,horizontal=F)
IFDist<-read.table(data)
xAkseli <- seq(PieninSade, MaxSade, by=SadeLisays)
for(kk in 1:12){
yAkseliHI <- NULL
yAkseliLO <- NULL
for(y in 1:(dim(IFDist)[2]/2)){
yAkseliHI[y] <- IFDist[kk,(2*y - 1)]
yAkseliLO[y] <- IFDist[kk,(2*y)]
}
plot(xAkseli,yAkseliHI, xlim=c(1,(MaxSade+ 2* SadeLisays)), ylim=c(0, PiirYlaRaja),
xlab="Säde", ylab="IF (high = punainen, low = musta)", col=rgb(1,0,0), main = paste("Piste ", Piste, " kuukausi ", kk, sep=""))
par(new=T)
plot(xAkseli,yAkseliLO, xlim=c(1,(MaxSade+ 2* SadeLisays)), ylim=c(0, PiirYlaRaja),
xlab="Säde", ylab="IF (high = punainen, low = musta)", col=rgb(0,0,0), main = paste("Piste ", Piste, " kuukausi ", kk, sep=""))
}
# musta on 0,0,0
dev.off()
}
#-------------------------------------------------------------------------------------------------------------
DataRuudut <- 10
# HUOM! tiedostoissa pääte .out on 50km ja .out.1 on 10km
if (DataRuudut == 10) paate <- ".out.1"
if (DataRuudut == 50) paate <- ".out"
SadeLisays <- 50
PieninSade <- 50
MaxSade <- 1700
LOLisays <- 0.1981
LOLisaysArvio <- 0.199
LALisays <- 0.08995
LALisaysArvio <- 0.09
PisteKK <- c(25,60,22,60.5,22,63,26,62,30,63,24.5,65,29.5,65,24,67.5,28,69)
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
PopulaatioData <- NULL
PopulaatioData <-
read.table(paste(polku,"Population_Europe\\",DataRuudut,"km\\ciesin_center_63lon_",DataRuudut,"km_2",sep=""))
PopulaatioData <- sPD
for( Piste in 1:1){
LO <- PisteKK[(2*Piste-1)]
LA <- PisteKK[(2*Piste)]
SummaPixCi <- NULL
# Seuraavassa loopissa kk on kuukausi
for(kk in 1:12){
tiedosto1<-paste("joinfilelist_",kk,"_PM_0_1.joined", paate,sep="")
tiedosto2<-paste("joinfilelist_",kk,"_PM_0_1_1.joined", paate,sep="")
tiedosto3<-paste("joinfilelist_",kk,"_PM_1_2_5.joined", paate,sep="")
PitoisuusData1 <- NULL
PitoisuusData1 <-read.table(paste(polku,"Concentration_Point\\",tiedosto1,sep=""))
PitoisuusData2 <- NULL
PitoisuusData2 <-read.table(paste(polku,"Concentration_Point\\",tiedosto2,sep=""))
PitoisuusData3 <- NULL
PitoisuusData3 <-read.table(paste(polku,"Concentration_Point\\",tiedosto3,sep=""))
PM2.5PitoisuusData<- NULL
PM2.5PitoisuusData <- PitoisuusData1 + PitoisuusData2 + PitoisuusData3
apuPM <- NULL
apuPM <- PM2.5PitoisuusData[(PM2.5PitoisuusData[,(1 + 2*Piste)]!=0 | PM2.5PitoisuusData[,(2 + 2*Piste)]!=0)
& PopulaatioData[,3]!=0,]
apuPD <- NULL
apuPD <- PopulaatioData[PopulaatioData[,3]!=0 &
(PM2.5PitoisuusData[,(1 + 2*Piste)]!=0 | PM2.5PitoisuusData[,(2 + 2*Piste)]!=0),]
apuHI <- NULL
apuLO <- NULL
apuSumma <- NULL
#apu1 <- rep(FALSE, dim(apuPD)[1])
sadeindeksi <- 0
for(Sade in seq(MaxSade, PieninSade, by= -SadeLisays )){
apu1 <- NULL
sadeindeksi <- sadeindeksi + 1
for(k in 1:dim(apuPD)[1]){
if( (apuPD[k,1] >= LO & apuPD[k,2] >= LA) | (apuPD[k,1] < LO & apuPD[k,2] < LA ) ){
apu1[k] <-
(Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2]),1) < Sade &
Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2] + LALisays),1) < Sade)
}
if((apuPD[k,1] < LO & apuPD[k,2] >= LA) | (apuPD[k,1] >= LO & apuPD[k,2] < LA ) ){
apu1[k] <-
(Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2]),1) < Sade &
Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2] + LALisays),1) < Sade)
}
}
apuPM <- apuPM[apu1,]
apuPD <- apuPD[apu1,]
apuHI[sadeindeksi] <- sum( apuPM[,(1 + 2*Piste)] * apuPD[,3] )
apuLO[sadeindeksi] <- sum( apuPM[,(2 + 2*Piste)] * apuPD[,3] )
aputaulu <- c(Sade, dim(apuPD)[1])
TalletusKohde <- paste(talletuspolku,"Piste_",Piste,"\\IF_Piste_",Piste,"_kk_",kk,"_Sade_",Sade,"_Distance.txt",sep="")
write.table(t(aputaulu), file = TalletusKohde, row.names = FALSE, col.names = FALSE)
}
ylaraja <- length(apuHI)
for(i in 1: ylaraja){
apuSumma[(2*i - 1)] <- apuHI[(ylaraja - i + 1)]
apuSumma[2*i] <- apuLO[(ylaraja - i + 1)]
}
SummaPixCi <- rbind(SummaPixCi, t(apuSumma))
}
# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
BR <- 20/(24 * 60 * 60)
# Q <- 3 * 1/(365 * 24 * 60 * 60)kg/s
Q <- 3/(365 * 24 * 60 * 60)
IF <- SummaPixCi * BR/Q
TalletusKohde <- paste(talletuspolku,"IF_Piste_",Piste,"_Distance.txt",sep="")
write.table(IF, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
}
write.table(t(t(PietariTot)), file = paste(talletuspolku,"Pietari_100km_totuus.txt",sep=""), row.names = FALSE, col.names = FALSE)
leikkausData <- NULL
osuus <- NULL
for(i in 1:dim(kojetaulu)[1]){
pisteLO <- NULL
pisteLA <- NULL
pisteLO <- kojetaulu[i,1]
pisteLA <- kojetaulu[i,2]
osuus[i] <- YmpyraNelikulmioAla(c(25,60),150, c(pisteLO,pisteLA,pisteLO+LOLisays,pisteLA,
pisteLO+LOLisays,pisteLA + LALisays,pisteLO,pisteLA+ LALisays))[3]
}
leikkausData <- osuus*kojetaulu[,3]
leik <- NULL
for(sade in 0:30){
leik[sade] <- YmpyraNelikulmioAla(c(25,60),sade, c(24.9272,59.9944,24.9272+LOLisays,59.9944,
24.9272+LOLisays,59.9944 + LALisays,24.9272,59.9944+ LALisays))[3]
}
pisteLO <- kojetaulu[4,1]
pisteLA <- kojetaulu[4,2]
leik <- NULL
for(sade in 0:30){
leik[sade] <- YmpyraNelikulmioAla(c(25,60),sade, c(pisteLO,pisteLA,pisteLO+LOLisays,pisteLA,
pisteLO+LOLisays,pisteLA + LALisays,pisteLO,pisteLA+ LALisays))[3]
}
KolmioYmpyranA25,61,1,24,59,25,61,26,59)
for(km in seq(0,20,by=0.2)){
cat(km)
cat(" ")
cat(KolmioYmpyranAla(1,1,0.5*km,5,0,2,3,-2,0))
cat("\n")
}
kp1 <- -5
kp2 <- -1500
r <- 1503
tilavek <- KolmionSisalla(kp1,kp2,10,1,-4,1,-4,8)
tilavek
PisteenAlue(kp1,kp2,tilavek[2:7])
KolmioYmpyranAla(kp1,kp2,r,10,1,-4,1,-4,8)
kp1 <- 5
kp2 <- -4
#Ongelma säteellä 9!!!!!!!!!!!!!!!!!!!!!
r <- 9
KolmioYmpyranAla(kp1,kp2,9,5,2,12,2,5,6)
tilavek <- KolmionSisalla(kp1,kp2,5,2,12,2,5,6)
paikkavek <- PisteenAlue(kp1,kp2,tilavek[2:7])
x1 <- kp1
x2 <- kp2
pistevek <- tilavek[2:7]
kp <- c(1,500)
kp <- c(-4,4)
PisteenAlue(kp[1],kp[2],KolmionSisalla(kp[1],kp[2],A[1],A[2],B[1],B[2],C[1],C[2])[2:7])
x1 <- kp[1]
x2 <- kp[2]
pistevek <- KolmionSisalla(kp[1],kp[2],A[1],A[2],B[1],B[2],C[1],C[2])[2:7]
|
ympyran_ja-monikulmion_leikkaus.r (16.6.2006)
#-----------------------------------------------------------------------
#Maan pinnalla sijaitsevien pisteiden etäisyyden laskeva funktio
#-----------------------------------------------------------------------
# maan säteen arviona voidaan käyttää 6371 km
# jos laskutapa on 1, niin lasketaan 2 pisteen välimatka pallon pintaa pitkin
# jos laskutapa on 2, niin lasketaan pisteiden (p1LO, p1LA) ja (p2LO, p2LA) euklidinen etäisyys,
# tällöin säteellä ei ole merkitystä
Etaisyys <- function(sade, p1LO, p1LA, p2LO, p2LA, laskutapa)
{
if(laskutapa == 1){
piste1LA <- p1LA * (2*pi/360)
piste1LO <- p1LO * (2*pi/360)
piste2LA <- p2LA * (2*pi/360)
piste2LO <- p2LO * (2*pi/360)
#pallokoorinaatistoa käyttäen
SuoraValiMatka <- sade * sqrt( (cos(piste1LO)*cos(piste1LA) - cos(piste2LO)*cos(piste2LA))^2
+ (sin(piste1LO)*cos(piste1LA) - sin(piste2LO)*cos(piste2LA))^2
+ (sin(piste1LA) - sin(piste2LA))^2 )
#cosini lause
Kulma <- acos(1 - SuoraValiMatka^2/(2*sade^2))
#kaaren pituus on kulma(radiaaneina)*säde
tulos <- sade * Kulma
}
if(laskutapa == 2)
{
tulos <- sqrt((p1LO - p2LO)^2 + (p1LA - p2LA)^2)
}
tulos
}
#----------------------------------------------------------------------------
#----------------------------------------------------------
#Kolmion pinta-ala, kun tiedetään sivujen pituudet
HeronAla <- function(a,b,c){
p <- (a+b+c)/2
sqrt(p*(p-a)*(p-b)*(p-c))
}
#Tavallisella tavalla laskettava kolmion pinta-ala
KantaKorkeus <- function(k, h){
1/2 * k * h
}
#--------------------------------------------------------------------------
#----------------------------------------------------------------------------
# Voi tulla ongelmia, jos piste on kulmapisteiden välisellä janalla !
# Jos palautetun vektorin 1. koordinaatti on 1, niin tarkasteltavan
# kolmion tilavuus on nolla. Muutoin ensimmäinen koordinaatti on 0, 2. koord.
# on A pisteen x-koord, 3. koord. on A pisteen y-koord., ..., 8. koord. on 1 jos
# piste on kolmion sisällä, muutoin 8. koord on nolla.
KolmionSisalla <- function(x1,x2,D1,D2,E1,E2,F1,F2){
xvek <- c(D1,E1,F1)
yvek <- c(D2,E2,F2)
jarjvek <- NULL
paha <- 0
#Huonot tapaukset
#pisteiden D ja E kautta kulkeva suora -> testataan ovatko pisteet samalla suoralla
if(xvek[1] != xvek[2]){
k <- (yvek[2] - yvek[1])/(xvek[2] - xvek[1])
y <- k*(xvek[3] - xvek[1]) + yvek[1]
if(yvek[3] == y){ paha <- 1}
}
# Kaksi pistettä samat tai samalla vaakasuoralla tai samalla pystysuoralla
if( (xvek[1] == xvek[2] & yvek[1] == yvek[2]) |
(xvek[1] == xvek[3] & yvek[1] == yvek[3]) |
(xvek[2] == xvek[3] & yvek[2] == yvek[3]) |
(xvek[1] == xvek[2] & xvek[1] == xvek[3]) |
(yvek[1] == yvek[2] & yvek[1] == yvek[3])
)
{
paha <- 1
}
jarjvek[1] <- paha
if(paha == 0)
{
apu1 <- min(xvek)
if( length(xvek[xvek==apu1])==1)
{
jarjvek[2] <- xvek[xvek == apu1]
jarjvek[3] <- yvek[xvek == apu1]
apu2 <- max(yvek[xvek != apu1])
if( length(yvek[(xvek!= apu1 & yvek == apu2)]) ==1)
{
jarjvek[4] <- xvek[(xvek!= apu1 & yvek == apu2)]
jarjvek[5] <- yvek[(xvek!= apu1 & yvek == apu2)]
jarjvek[6] <- xvek[(xvek!= apu1 & yvek != apu2)]
jarjvek[7] <- yvek[(xvek!= apu1 & yvek != apu2)]
}
if( length(yvek[(xvek!= apu1 & yvek == apu2)]) ==2)
{
apu3 <- min(xvek[xvek!= apu1])
jarjvek[4] <- xvek[(xvek!= apu1 & xvek == apu3)]
jarjvek[5] <- yvek[(xvek!= apu1 & xvek == apu3)]
jarjvek[6] <- xvek[(xvek!= apu1 & xvek != apu3)]
jarjvek[7] <- yvek[(xvek!= apu1 & xvek != apu3)]
}
#Jos piste C on janan AB yläpuolella, niin vaihdetaan B ja C
janaT <- NULL
janaT[1] <- (jarjvek[6] - jarjvek[2])/(jarjvek[4] - jarjvek[2])
janaY <- NULL
janaY[1] <- jarjvek[5] * janaT[1] + (1 - janaT[1])* jarjvek[3]
if(jarjvek[7] > janaY[1]){
apu4 <- jarjvek
jarjvek[4:5] <- apu4[6:7]
jarjvek[6:7] <- apu4[4:5]
}
}
if( length(xvek[xvek==apu1])==2)
{
apu2 <- min(yvek[xvek == apu1])
jarjvek[2] <- xvek[(xvek == apu1 & yvek == apu2)]
jarjvek[3] <- yvek[(xvek == apu1 & yvek == apu2)]
jarjvek[4] <- xvek[(xvek == apu1 & yvek != apu2)]
jarjvek[5] <- yvek[(xvek == apu1 & yvek != apu2)]
jarjvek[6] <- xvek[xvek != apu1 ]
jarjvek[7] <- yvek[xvek != apu1 ]
}
}
if(jarjvek[1]==0)
{
#Janan koordinaatit aina järjestyksessä AB, BC, AC
janaT <- NULL
janaT[1] <- (x1 - jarjvek[2])/(jarjvek[4] - jarjvek[2])
janaT[2] <- (x1 - jarjvek[4])/(jarjvek[6] - jarjvek[4])
janaT[3] <- (x1 - jarjvek[2])/(jarjvek[6] - jarjvek[2])
janaY <- NULL
janaY[1] <- jarjvek[5] * janaT[1] + (1 - janaT[1])* jarjvek[3]
janaY[2] <- jarjvek[7] * janaT[2] + (1 - janaT[2])* jarjvek[5]
janaY[3] <- jarjvek[7] * janaT[3] + (1 - janaT[3])* jarjvek[3]
if(jarjvek[2]==jarjvek[4])
{
jarjvek[8] <- 0
if((jarjvek[2]<= x1) & (x2 <= janaY[2]) & (janaY[3]<= x2))
{
jarjvek[8] <- 1
}
}
if(jarjvek[2]!=jarjvek[4])
{
if(jarjvek[4]==jarjvek[6])
{
jarjvek[8] <- 0
if((x1 <= jarjvek[4]) & (x2 <= janaY[1]) & (janaY[3]<= x2))
{
jarjvek[8] <- 1
}
}
if(jarjvek[4]!=jarjvek[6])
{
jarjvek[8] <- 0
if( ((jarjvek[6] > jarjvek[4]) & (x2 <= janaY[1]) &
(x2 <= janaY[2]) & (janaY[3] <= x2)) |
((jarjvek[6] < jarjvek[4]) & (x2 <= janaY[1]) &
(janaY[2] <= x2) & (janaY[3] <= x2))
)
{
jarjvek[8] <- 1
}
}
}
}
jarjvek
}
#----------------------------------------------------------------------------
#------------------------------------------------------------------------
#Pistevektoriin syötetään KolmionSisalla tuottama vektori!!!!!
# x1, x2 on ympyrän keskipiste
PisteenAlue <- function(x1,x2, pistevek){
alue <- c(0,0,0)
# C Kärkenä piste C
A <- pistevek[1:2]
B <- pistevek[3:4]
C <- pistevek[5:6]
t1 <-
((x2 - C[2])*(A[1] - C[1]) - (x1 - C[1])*(A[2] - C[2]))/
((x1 - C[1])*(B[2] - A[2]) - (x2 - C[2])*(B[1] - A[1]))
if(t1==0){
suunta <- -1
if((A[1] - C[1])!=0 & ((x1 - C[1])/(A[1] - C[1]))>=0)
suunta <- 1
if((A[1] - C[1])==0 & ((x2 - C[2])/(A[2] - C[2]))>=0)
suunta <- 1
t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, A[1],A[2], C[1],C[2],2)* suunta
}
if(t1==1){
suunta <- -1
if((B[1] - C[1])!=0 & ((x1 - C[1])/(B[1] - C[1]))>=0)
suunta <- 1
if((B[1] - C[1])==0 & ((x2 - C[2])/(B[2] - C[2]))>=0)
suunta <- 1
t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, B[1],B[2], C[1],C[2],2)*suunta
}
if(t1!=0 & t1!=1)
ifelse(x1==C[1],
{
suunta <- -1
if(sum(C - t1*B + (1 - t1)*A)/sum(C - c(x1,x2))>=0)
suunta <- 1
t2 <- Etaisyys( 6371, x1,x2,C[1],C[2],2)/
Etaisyys(6371,C[1],C[2],t1*B[1] + (1 - t1)*A[1],t1*B[2] + (1 - t1)*A[2],2) *suunta
},
t2 <-
(x1 - C[1])/(t1*B[1] + (1 - t1)*A[1] - C[1])
)
if(0<= t1 & t1 <= 1 & 1 < t2)
{
alue[1] <- 1
alue[2] <- t1
alue[3] <- t2
}
if(0< t1 & t1 < 1 & t2 < 0)
{
alue[1] <- 4
alue[2] <- t1
alue[3] <- t2
}
if(alue[1] == 0){
# C Kärkenä piste A
#---------------------------
A <- pistevek[3:4]
B <- pistevek[5:6]
C <- pistevek[1:2]
t1 <-
((x2 - C[2])*(A[1] - C[1]) - (x1 - C[1])*(A[2] - C[2]))/
((x1 - C[1])*(B[2] - A[2]) - (x2 - C[2])*(B[1] - A[1]))
if(t1==0){
suunta <- -1
if((A[1] - C[1])!=0 & ((x1 - C[1])/(A[1] - C[1]))>=0)
suunta <- 1
if((A[1] - C[1])==0 & ((x2 - C[2])/(A[2] - C[2]))>=0)
suunta <- 1
t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, A[1],A[2], C[1],C[2],2)* suunta
}
if(t1==1){
suunta <- -1
if((B[1] - C[1])!=0 & ((x1 - C[1])/(B[1] - C[1]))>=0)
suunta <- 1
if((B[1] - C[1])==0 & ((x2 - C[2])/(B[2] - C[2]))>=0)
suunta <- 1
t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, B[1],B[2], C[1],C[2],2)*suunta
}
if(t1!=0 & t1!=1)
ifelse(x1==C[1],
{
suunta <- -1
if(sum(C - t1*B + (1 - t1)*A)/sum(C - c(x1,x2))>=0)
suunta <- 1
t2 <- Etaisyys( 6371, x1,x2,C[1],C[2],2)/
Etaisyys(6371,C[1],C[2],t1*B[1] + (1 - t1)*A[1],t1*B[2] + (1 - t1)*A[2],2) *suunta
},
t2 <-
(x1 - C[1])/(t1*B[1] + (1 - t1)*A[1] - C[1])
)
#--------------------------------------
if(0<= t1 & t1 <= 1 & 1 < t2)
{
alue[1] <- 3
alue[2] <- t1
alue[3] <- t2
}
if(0< t1 & t1 < 1 & t2 < 0)
{
alue[1] <- 6
alue[2] <- t1
alue[3] <- t2
}
if(alue[1] == 0){
# C Kärkenä piste B
#-----------------------------------
A <- pistevek[5:6]
B <- pistevek[1:2]
C <- pistevek[3:4]
t1 <-
((x2 - C[2])*(A[1] - C[1]) - (x1 - C[1])*(A[2] - C[2]))/
((x1 - C[1])*(B[2] - A[2]) - (x2 - C[2])*(B[1] - A[1]))
if(t1==0){
suunta <- -1
if((A[1] - C[1])!=0 & ((x1 - C[1])/(A[1] - C[1]))>=0)
suunta <- 1
if((A[1] - C[1])==0 & ((x2 - C[2])/(A[2] - C[2]))>=0)
suunta <- 1
t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, A[1],A[2], C[1],C[2],2)* suunta
}
if(t1==1){
suunta <- -1
if((B[1] - C[1])!=0 & ((x1 - C[1])/(B[1] - C[1]))>=0)
suunta <- 1
if((B[1] - C[1])==0 & ((x2 - C[2])/(B[2] - C[2]))>=0)
suunta <- 1
t2 <- Etaisyys(6371, x1,x2, C[1],C[2],2)/
Etaisyys(6371, B[1],B[2], C[1],C[2],2)*suunta
}
if(t1!=0 & t1!=1)
ifelse(x1==C[1],
{
suunta <- -1
if(sum(C - t1*B + (1 - t1)*A)/sum(C - c(x1,x2))>=0)
suunta <- 1
t2 <- Etaisyys( 6371, x1,x2,C[1],C[2],2)/
Etaisyys(6371,C[1],C[2],t1*B[1] + (1 - t1)*A[1],t1*B[2] + (1 - t1)*A[2],2) *suunta
},
t2 <-
(x1 - C[1])/(t1*B[1] + (1 - t1)*A[1] - C[1])
)
#---------------------------------------------------
if(0<= t1 & t1 <= 1 & 1 < t2)
{
alue[1] <- 5
alue[2] <- t1
alue[3] <- t2
}
if(0< t1 & t1 < 1 & t2 < 0)
{
alue[1] <- 2
alue[2] <- t1
alue[3] <- t2
}
}
}
alue
}
#----------------------------------------------------------------------------
#kp1, kp2, r, A[1],A[2],C[1],C[2]
#kulm <- c(B[1],B[2],C[1],C[2])
#kv1 <- kulm[1]
#kv2 <- kulm[2]
#kv3 <- kulm[3]
#kv4 <- kulm[4]
#----------------------------------------------------------------------------
#-----------------------------------------------------------------------------
#Laskee ympyrän ja kolmion leikkauksen alan, kun ympyrän keskipiste on yksi
#kolmion kulmista
#Palauttaa vektorin jonka 1. komponentti on leikkausen pinta-ala
# 2. komp. kolmion pinta-ala, 3. on näiden suhde ja 4. kolmion korkeus.
AlaYmpyranKpKolmionKulma <- function(kp1, kp2, r, kv1,kv2,kv3,kv4)
{
ala <- 0
vek <- c(kp1, kp2, r, kv1,kv2,kv3,kv4)
tilavek <- KolmionSisalla(vek[1],vek[2],vek[4],vek[5],vek[6], vek[7],vek[1],vek[2])
if(tilavek[1]==1 ){
tulos <- c(0,0,0,0)
}
# Oikea kolmio
if(tilavek[1]==0){
hNollaraja <- 10^(-12)
kp <- c(kp1,kp2)
A <- tilavek[2:3]
B <- tilavek[4:5]
C <- tilavek[6:7]
if(kp[1]==A[1] & kp[2]==A[2]){
A <- B
B <- C
}
if(kp[1]==B[1] & kp[2]==B[2]){
B <- A
A <- C
}
AB <- Etaisyys(6371,A[1],A[2],B[1],B[2],1)
kpB <- Etaisyys(6371,kp[1],kp[2],B[1],B[2],1)
kpA <- Etaisyys(6371,kp[1],kp[2],A[1],A[2],1)
cosAlpha <- (kpA^2 + AB^2 - kpB^2)/(2*AB*kpA)
cosBeta <- (kpB^2 + AB^2 - kpA^2)/(2*AB*kpB)
cosGamma <- (kpA^2 + kpB^2 - AB^2)/(2*kpB*kpA)
if( cosAlpha >= 0)
{
h <- sqrt(kpA^2 - kpA^2*cosAlpha^2)
#A ja B eivät kuulu ympyrään ja ympyrä ei leikkaa niiden välistä janaa.
if(r <= h | (h < r & r <= kpB & cosBeta < 0)){
ala <- 1/2 * r^2 * acos(cosGamma)
}
# A ja B eivät kuulu ympyrään ja ympyrä leikkaa niiden välistä janaa.
if(h < r & kpA > r & kpB > r){
if(cosBeta <=0) ala <- 1/2*r^2*acos(cosGamma)
if(cosBeta > 0)
ala <- 1/2 * r^2*(acos(h/kpB) - acos(h/r)) +
KantaKorkeus(2*sqrt(r^2 - h^2),h) +
1/2 * r^2* (acos(h/kpA) - acos(h/r))
}
# B kuuluu ympyrään ja A ei kuulu
if(h < r & kpB <= r & kpA > r){
if(cosBeta > 0)
x <- sqrt(abs(kpB^2 - h^2)) + sqrt(r^2 - h^2)
if(cosBeta <= 0)
x <- sqrt(r^2 - h^2) - sqrt(abs(kpB^2 - h^2))
ala <- KantaKorkeus(x,h) + 1/2 * r^2 * (acos(h/kpA) - acos(h/r))
}
# B ei kuulu ympyrään ja A kuuluu
if(h < r & kpB > r & kpA <= r){
x <- sqrt(abs(kpA^2 - h^2)) + sqrt(r^2 - h^2)
ala <- KantaKorkeus(x,h) + 1/2 * r^2 * (acos(h/kpB) - acos(h/r))
}
# A ja B kuuluvat ympyrään
if(h < r & kpA <= r & kpB <= r){
ala <- KantaKorkeus(AB, h)
}
tulos <- c(ala, (1/2*h*AB), (ala/(1/2*h*AB)),h)
}
if( cosAlpha < 0)
{
h <- sqrt(kpB^2 - kpB^2*cosBeta^2)
# A ja B eivät kuulu ympyrään ja ympyrä ei leikkaa niiden välistä janaa.
if(r <= h | (h < r & r <= kpA)){
ala <- 1/2 * r^2 * acos(cosGamma)
}
# A ja B eivät kuulu ympyrään ja ympyrä leikkaa niiden välistä janaa.
# Tämä ei ole mahdollista, koska Alpha > 90 astetta
# B kuuluu ympyrään ja A ei kuulu
# Ei mahdollista, koska kpA < kpB
# B ei kuulu ympyrään ja A kuuluu
if(h < r & kpB > r & kpA <= r){
x <- sqrt(r^2 - h^2) - sqrt(abs(kpA^2 - h^2))
ala <- KantaKorkeus(x,h) + 1/2 * r^2 * (acos(h/kpB) - acos(h/r))
}
# A ja B kuuluvat ympyrään
if(h < r & kpA <= r & kpB <= r){
ala <- KantaKorkeus(AB, h)
}
tulos <- c(ala, (1/2*h*AB), (ala/(1/2*h*AB)),h)
}
}
tulos
}
#----------------------------------------------------------------------------
#leik <- NULL
# for(sade in 0:30){
# leik[sade+1] <- YmpyraNelikulmioAla(c(25,60),sade, c(24.9272,59.9944,24.9272+LOLisays,59.9944,
# 24.9272+LOLisays,59.9944 + LALisays,24.9272,59.9944+ LALisays))[3]
# }
#leik
#pisteLO <- kojetaulu[4,1]
#pisteLA <- kojetaulu[4,2]
#leik <- NULL
#for(sade in 0:30){
#leik[sade+1] <- YmpyraNelikulmioAla(c(25,60),sade, c(pisteLO,pisteLA,pisteLO+LOLisays,pisteLA,
#pisteLO+LOLisays,pisteLA + LALisays,pisteLO,pisteLA+ LALisays))[3]
#}
#leik
#kp <- c(25,60)
#kulmat <- c(pisteLO,pisteLA,pisteLO+LOLisays,pisteLA,
#pisteLO+LOLisays,pisteLA + LALisays,pisteLO,pisteLA+ LALisays)
#r <- 1
#KolmioYmpyranAla(kp[1], kp[2], r, kulmat[1],kulmat[2],kulmat[3],kulmat[4],kulmat[7],kulmat[8])
#KolmioYmpyranAla(kp[1], kp[2], r, kulmat[5],kulmat[6],kulmat[3],kulmat[4],kulmat[7],kulmat[8])
#kulmavek <- c(kulmat[1],kulmat[2],kulmat[3],kulmat[4],kulmat[7],kulmat[8])
#r <- 1
#kp1 <- kp[1]
#kp2 <- kp[2]
#kv1 <- kulmavek[1]
#kv2 <- kulmavek[2]
#kv3 <- kulmavek[3]
#kv4 <- kulmavek[4]
#kv5 <- kulmavek[5]
#kv6 <- kulmavek[6]
#----------------------------------------------------------------------------
#----------------------------------------------------------------------------
#kv on vektori, joka sisaltaa kolmion kulmat
KolmioYmpyranAla <- function(kp1, kp2, r, kv1,kv2,kv3,kv4,kv5,kv6)
{
ala <- NULL
tilavek <- KolmionSisalla(kp1,kp2,kv1,kv2,kv3,kv4,kv5,kv6)
A <- tilavek[2:3]
B <- tilavek[4:5]
C <- tilavek[6:7]
AB <- Etaisyys(6371,A[1],A[2],B[1],B[2],1)
BC <- Etaisyys(6371,B[1],B[2],C[1],C[2],1)
AC <- Etaisyys(6371,A[1],A[2],C[1],C[2],1)
#Lasketaan tarvittavien kolmioiden leikkaus
alaKPAB <- AlaYmpyranKpKolmionKulma(kp1, kp2, r, A[1],A[2],B[1],B[2])[1]
alaKPBC <- AlaYmpyranKpKolmionKulma(kp1, kp2, r, B[1],B[2],C[1],C[2])[1]
alaKPAC <- AlaYmpyranKpKolmionKulma(kp1, kp2, r, A[1],A[2],C[1],C[2])[1]
#oikea kolmio ja piste kolmion sisalla
if(tilavek[1]==0 & tilavek[8]==1)
{
ala[1] <- alaKPAB + alaKPAC + alaKPBC
ala[2] <- HeronAla(AB,BC,AC)
if(ala[1] < 0) ala[1] <- 0
if(ala[1] > ala[2]) ala[1] <- ala[2]
ala[3] <- ala[1]/ala[2]
}
#Oikea kolmio ja piste kolmion ulkopuolella
if(tilavek[1]==0 & tilavek[8]==0)
{
paikkavek <- PisteenAlue(kp1,kp2,tilavek[2:7])
if(paikkavek[1] == 1)
{
ala[1] <- alaKPAC + alaKPBC - alaKPAB
}
if(paikkavek[1] == 2)
{
ala[1] <- alaKPAC - alaKPAB - alaKPBC
}
if(paikkavek[1] == 3)
{
ala[1] <- alaKPAB + alaKPAC - alaKPBC
}
if(paikkavek[1] == 4)
{
ala[1] <- alaKPAB - alaKPAC - alaKPBC
}
if(paikkavek[1] == 5)
{
ala[1] <- alaKPAB + alaKPBC - alaKPAC
}
if(paikkavek[1] == 6)
{
ala[1] <- alaKPBC - alaKPAC - alaKPAB
}
ala[2] <- HeronAla(AB,BC,AC)
if(ala[1] < 0) ala[1] <- 0
if(ala[1] > ala[2]) ala[1] <- ala[2]
ala[3] <- ala[1]/ala[2]
}
ala
}
#----------------------------------------------------------------------------
#----------------------------------------------------------------------------
# Syötä kulmat kiertojärjestyksessä vastapäivään!!!!!!
YmpyraNelikulmioAla <-function(kp,r,kulmat)
{
tulos <- NULL
apu1 <- KolmioYmpyranAla(kp[1], kp[2], r, kulmat[1],kulmat[2],kulmat[3],kulmat[4],kulmat[7],kulmat[8])
apu2 <- KolmioYmpyranAla(kp[1], kp[2], r, kulmat[5],kulmat[6],kulmat[3],kulmat[4],kulmat[7],kulmat[8])
tulos[1:2] <- apu1[1:2] + apu2[1:2]
tulos[3] <- tulos[1]/tulos[2]
tulos
}
#-----------------------------------------------------------------------------
#leikkausData <- NULL
#osuus <- NULL
#for(i in 1:dim(kojetaulu3)[1]){
#for(i in 1:100){
#pisteLO <- NULL
#pisteLA <- NULL
#pisteLO <- kojetaulu[i,1]
#pisteLA <- kojetaulu[i,2]
#osuus[i] <- YmpyraNelikulmioAla(c(25,60),25, c(pisteLO,pisteLA,pisteLO+LOLisays,pisteLA,
#pisteLO+LOLisays,pisteLA + LALisays,pisteLO,pisteLA+ LALisays))[3]
#}
#leikkausData <- osuus*kojetaulu[,3]
#cbind(kojetaulu[,3],t(t(leikkausData)))
#summary(cbind(kojetaulu3[,3],t(t(leikkausData)))[,2])
#-----------------------------------------------------------------------------
#Laskeetaan pienempiä ympyröitä osuuksien avulla.
DataRuudut <- 10
# HUOM! tiedostoissa pääte .out on 50km ja .out.1 on 10km
if (DataRuudut == 10) paate <- ".out.1"
if (DataRuudut == 50) paate <- ".out"
SadeLisays <- 50
PieninSade <- 0
MaxSade <- 1700
LOLisays <- 0.1981
LOLisaysArvio <- 0.199
LALisays <- 0.08995
LALisaysArvio <- 0.09
PisteKK <- c(25,60,22,60.5,22,63,26,62,30,63,24.5,65,29.5,65,24,67.5,28,69)
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
PopulaatioData <- NULL
PopulaatioData <-
read.table(paste(polku,"Population_Europe\\",DataRuudut,"km\\ciesin_center_63lon_",DataRuudut,"km_2",sep=""))
PopulaatioData <- sPD
for( Piste in 4:4){
LO <- PisteKK[(2*Piste-1)]
LA <- PisteKK[(2*Piste)]
SummaPixCi <- NULL
# Seuraavassa loopissa kk on kuukausi
for(kk in 1:12){
tiedosto1<-paste("joinfilelist_",kk,"_PM_0_1.joined", paate,sep="")
tiedosto2<-paste("joinfilelist_",kk,"_PM_0_1_1.joined", paate,sep="")
tiedosto3<-paste("joinfilelist_",kk,"_PM_1_2_5.joined", paate,sep="")
PitoisuusData1 <- NULL
PitoisuusData1 <-read.table(paste(polku,"Concentration_Point\\",tiedosto1,sep=""))
PitoisuusData2 <- NULL
PitoisuusData2 <-read.table(paste(polku,"Concentration_Point\\",tiedosto2,sep=""))
PitoisuusData3 <- NULL
PitoisuusData3 <-read.table(paste(polku,"Concentration_Point\\",tiedosto3,sep=""))
PM2.5PitoisuusData<- NULL
PM2.5PitoisuusData <- PitoisuusData1 + PitoisuusData2 + PitoisuusData3
apuPM <- NULL
apuPM <- PM2.5PitoisuusData[(PM2.5PitoisuusData[,(1 + 2*Piste)]!=0 | PM2.5PitoisuusData[,(2 + 2*Piste)]!=0)
& PopulaatioData[,3]!=0,]
apuPD <- NULL
apuPD <- PopulaatioData[PopulaatioData[,3]!=0 &
(PM2.5PitoisuusData[,(1 + 2*Piste)]!=0 | PM2.5PitoisuusData[,(2 + 2*Piste)]!=0),]
apuHI <- NULL
apuLO <- NULL
apuSumma <- NULL
sadeindeksi <- 0
for(Sade in seq(MaxSade, PieninSade, by= -SadeLisays )){
apu1 <- NULL
sadeindeksi <- sadeindeksi + 1
#-------------------------------------------------------------
#Etäisyyden laskentatapa
for(k in 1:dim(apuPD)[1]){
if( (apuPD[k,1] >= LO & apuPD[k,2] >= LA) | (apuPD[k,1] < LO & apuPD[k,2] < LA ) ){
apu1[k] <-
(Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2]),1) < Sade &
Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2] + LALisays),1) < Sade)
}
if((apuPD[k,1] < LO & apuPD[k,2] >= LA) | (apuPD[k,1] >= LO & apuPD[k,2] < LA ) ){
apu1[k] <-
(Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2]),1) < Sade &
Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2] + LALisays),1) < Sade )
}
}
# for(k in 1:dim(apuPD)[1]){
# if( (apuPD[k,1] >= LO & apuPD[k,2] >= LA) | (apuPD[k,1] < LO & apuPD[k,2] < LA ) ){
# apu1[k] <-
# (Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2]),1) < Sade+25 &
# Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2] + LALisays),1) < Sade+25)
# }
# if((apuPD[k,1] < LO & apuPD[k,2] >= LA) | (apuPD[k,1] >= LO & apuPD[k,2] < LA ) ){
# apu1[k] <-
# (Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2]),1) < Sade+25 &
# Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2] + LALisays),1) < Sade+25)
# }
# }
apuPM <- apuPM[apu1,]
apuPD <- apuPD[apu1,]
# osuus <- NULL
# apuPDOsuus <- NULL
# for(i in 1:dim(apuPD)[1]){
# pisteLO <- NULL
# pisteLA <- NULL
# pisteLO <- apuPD[i,1]
# pisteLA <- apuPD[i,2]
# osuus[i] <- YmpyraNelikulmioAla(c(LO,LA),Sade, c(pisteLO,pisteLA,pisteLO+LOLisays,pisteLA,
# pisteLO+LOLisays,pisteLA + LALisays,pisteLO,pisteLA+ LALisays))[3]
# }
# apuPDOsuus <- apuPD[,3]*t(t(osuus))
#
# apuHI[sadeindeksi] <- sum( apuPM[,(1 + 2*Piste)] * apuPDOsuus )
# apuLO[sadeindeksi] <- sum( apuPM[,(2 + 2*Piste)] * apuPDOsuus )
apuHI[sadeindeksi] <- sum( apuPM[,(1 + 2*Piste)] * apuPD[,3] )
apuLO[sadeindeksi] <- sum( apuPM[,(2 + 2*Piste)] * apuPD[,3] )
#Etäisyyden laskentatapa
#--------------------------------------------------------------------------
aputaulu <- c(Sade, dim(apuPD)[1])
# TalletusKohde <- paste(talletuspolku,"Piste_",Piste,"\\IF_Piste_",Piste,"_kk_",kk,"_Sade_",Sade,"_Distance_Osuus_kokeilu.txt",sep="")
TalletusKohde <- paste(talletuspolku,"Piste_",Piste,"\\IF_Piste_",Piste,"_kk_",kk,"_Sade_",Sade,"_Distance.txt",sep="")
write.table(t(aputaulu), file = TalletusKohde, row.names = FALSE, col.names = FALSE)
}
ylaraja <- length(apuHI)
for(i in 1: ylaraja){
apuSumma[(2*i - 1)] <- apuHI[(ylaraja - i + 1)]
apuSumma[2*i] <- apuLO[(ylaraja - i + 1)]
}
SummaPixCi <- rbind(SummaPixCi, t(apuSumma))
}
# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
BR <- 20/(24 * 60 * 60)
# Q <- 3 * 1/(365 * 24 * 60 * 60)kg/s
Q <- 3/(365 * 24 * 60 * 60)
IF <- SummaPixCi * BR/Q
TalletusKohde <- paste(talletuspolku,"IF_Piste_",Piste,"_Distance.txt",sep="")
# TalletusKohde <- paste(talletuspolku,"IF_Piste_",Piste,"_Distance_Osuus_Kokeilu_SadeVali_1_Suurin_100.txt",sep="")
write.table(IF, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
}
#-----------------------------------------------------------------------------
# valitaan ruudut annetulla säteellä.
DataRuudut <- 10
MaxSade
# HUOM! tiedostoissa pääte .out on 50km ja .out.1 on 10km
if (DataRuudut == 10) paate <- ".out.1"
if (DataRuudut == 50) paate <- ".out"
SadeLisays <- 5
PieninSade <- 0
MaxSade <- 100
LOLisays <- 0.1981
LOLisaysArvio <- 0.199
LALisays <- 0.08995
LALisaysArvio <- 0.09
PisteKK <- c(25,60,22,60.5,22,63,26,62,30,63,24.5,65,29.5,65,24,67.5,28,69)
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Mallit\\IntakeFraction\\Distance\\"
PopulaatioData <- NULL
PopulaatioData <-
read.table(paste(polku,"Population_Europe\\",DataRuudut,"km\\ciesin_center_63lon_",DataRuudut,"km_2",sep=""))
apuPD <- NULL
apuPD <- PopulaatioData[PopulaatioData[,3]!=0 ,]
aputaulu <- NULL
for(Sade in seq(MaxSade, PieninSade, by= -SadeLisays )){
apu1 <- NULL
apu1[1] <- Sade
for(k in 1:dim(apuPD)[1]){
if( (apuPD[k,1] >= LO & apuPD[k,2] >= LA) | (apuPD[k,1] < LO & apuPD[k,2] < LA ) ){
apu1[k] <-
(Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2]),1) < Sade+25 &
Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2] + LALisays),1) < Sade+25)
}
if((apuPD[k,1] < LO & apuPD[k,2] >= LA) | (apuPD[k,1] >= LO & apuPD[k,2] < LA ) ){
apu1[k] <-
(Etaisyys(6371,LO, LA, (apuPD[k,1]), (apuPD[k,2]),1) < Sade+25 &
Etaisyys(6371,LO, LA, (apuPD[k,1] + LOLisays), (apuPD[k,2] + LALisays),1) < Sade+25)
}
}
aputaulu <- rbind(aputaulu, apu1)
TalletusKohde <- paste(talletuspolku,"PopulData_Sateittain.txt",sep="")
write.table(t(aputaulu), file = TalletusKohde, row.names = FALSE, col.names = FALSE)
}
#---------------------------------------------------------------------------------------------------------
kulmat <- c(24.9272,59.9944,24.9272+LOLisays,59.9944,
24.9272+LOLisays,59.9944 + LALisays,24.9272,59.9944+ LALisays)
kp <- c(25,60)
r <- 1
leik <- NULL
for(sade in 0:50){
leik[sade] <- YmpyraNelikulmioAla(c(25,60),sade, c(24.9272,59.9944,24.9272+LOLisays,59.9944,
24.9272+LOLisays,59.9944 + LALisays,24.9272,59.9944+ LALisays))[3]
}
leik
kulmat <- c(24.9272,59.9944,24.9272+LOLisays,59.9944,
24.9272+LOLisays,59.9944 + LALisays,24.9272,59.9944+ LALisays)
kp <- c(25,60)
r<- 7
leikkausData <- NULL
leikkausData <- t(t(kojetaulu2[,3]))
for(Sade in seq(1, 40, by= 1 )){
osuus <- NULL
for(i in 1:dim(leikkausData)[1]){
pisteLO <- NULL
pisteLA <- NULL
pisteLO <- kojetaulu2[i,1]
pisteLA <- kojetaulu2[i,2]
osuus[i] <- YmpyraNelikulmioAla(c(LO,LA),Sade, c(pisteLO,pisteLA,pisteLO+LOLisays,pisteLA,
pisteLO+LOLisays,pisteLA + LALisays,pisteLO,pisteLA+ LALisays))[3]
}
leikkausData <- cbind(leikkausData,kojetaulu2[,3]*t(t(osuus)))
}
alue <- NULL
for(i in 1:41){
alue[i] <- sum(leikkausData[,i])
}
#---------------------------------------------------------
#---------------------------------------------------------
#
# Vanhaa tavaraa!!!!
#-----------------------------------------------------------------------------
#Laskee ympyrän ja kolmion leikkauksen alan, kun ympyrän keskipiste on kolmion sisällä
AlaYmpyranKpKolmionSisalla <- function(kp1, kp2, r, kv1,kv2,kv3,kv4,kv5,kv6)
{
tilavek <- KolmionSisalla(kp1,kp2,kv1,kv2,kv3,kv4,kv5,kv6)
hNollaraja <- 10^(-6)
kp <- c(kp1,kp2)
A <- tilavek[2:3]
B <- tilavek[4:5]
C <- tilavek[6:7]
KolmionKulmat <- rbind(A,B,C,A)
#oikea kolmioa ja piste kolmion sisalla
if(tilavek[1]==0 & tilavek[8]==1)
{
AB <- Etaisyys(6371,tilavek[2],tilavek[3],tilavek[4],tilavek[5],1)
BC <- Etaisyys(6371,tilavek[4],tilavek[5],tilavek[6],tilavek[7],1)
AC <- Etaisyys(6371,tilavek[2],tilavek[3],tilavek[6],tilavek[7],1)
etA <- Etaisyys(6371, kp1, kp2, tilavek[2], tilavek[3],1)
etB <- Etaisyys(6371, kp1, kp2, tilavek[4], tilavek[5],1)
etC <- Etaisyys(6371, kp1, kp2, tilavek[6], tilavek[7],1)
etv <- c(etA, etB, etC, etA)
valit <- c(AB, BC, AC)
#h <- c(-100,-100,-100)
h <- c(1,1,1)
# Tapaukset, joissa ymp. kp on kolmion kulmapiste
if( (kp == A)[1] & (kp == A)[2]){ h[1] <- 0
h[3] <- 0
}
if( (kp == B)[1] & (kp == B)[2]){ h[1] <- 0
h[2] <- 0
}
if( (kp == C)[1] & (kp == C)[2]){ h[2] <- 0
h[3] <- 0
}
# Tapaukset, joissa ymp. kp on kolmion kulmapisteiden välisellä suoralla
if(A[1]!=B[1])
{
k <- (B[2] - A[2])/(B[1] - A[1])
y <- k* (kp[1] - A[1]) + A[2]
if( abs(kp[2] - y) < hNollaraja) h[1] <- 0
}
if( (A[1] == B[1] & kp[1]==A[1]))
{
h[1] <- 0
}
if(B[1]!=C[1])
{
k <- (C[2] - B[2])/(C[1] - B[1])
y <- k* (kp[1] - B[1]) + B[2]
if( abs(kp[2] - y) < hNollaraja) h[2] <- 0
}
if( (B[1] == C[1] & kp[1]==B[1]))
{
h[2] <- 0
}
if(C[1]!=A[1])
{
k <- (A[2] - C[2])/(A[1] - C[1])
y <- k* (kp[1] - C[1]) + C[2]
if( abs(kp[2] - y) < hNollaraja) h[3] <- 0
}
if( (C[1] == A[1] & kp[1]==C[1]))
{
h[3] <- 0
}
ala <- NULL
for(i in 1:3){
if(h[i] >= hNollaraja){
ala[i] <- AlaYmpyranKpKolmionKulma(kp[1], kp[2], r, KolmionKulmat[i,1],
KolmionKulmat[i,2],KolmionKulmat[(i+1),1],KolmionKulmat[(i+1),2])[1]
}
if(h[i] < hNollaraja){
ala[i] <- 0
}
}
tulos <- NULL
tulos[1] <- sum(ala)
tulos[2] <- KantaKorkeus(valit[1],
sqrt(valit[3]^2 - ((valit[1]^2 + valit[3]^2 - valit[2]^2)/(2*valit[1]))^2))
tulos[3] <- tulos[1]/tulos[2]
}
if(tilavek[1]==1 & length(tilavek)==1)
{
tulos <- c(0,0,0)
}
if(length(tilavek) > 1 & tilavek[8]== 0)
{
tulos <- c(-1,-1,-1)
}
tulos
}
#-------------------------------------------------------------------------------------------
#Pisteiden koordinaatit syötettävä kiertosuunnassa vastapäivään!!!!
YmpyraNelikulmioKpSisalla <- function(kpLO, kpLA, r, ALO, ALA, BLO, BLA, CLO, CLA, DLO, DLA)
{
tyovek <- c(kpLO, kpLA, r, ALO, ALA, BLO, BLA, CLO, CLA, DLO, DLA)
kp <- tyovek[1:2]
KolmionKulmat <- rbind(c(tyovek[4:5],tyovek[6:7]),
c(tyovek[6:7],tyovek[8:9]),
c(tyovek[8:9],tyovek[10:11]),
c(tyovek[10:11],tyovek[4:5]))
ala <- NULL
for(i in 1:4){
ala[i] <- AlaYmpyranKpKolmionKulma(kp[1], kp[2], r,
KolmionKulmat[i,1], KolmionKulmat[i,2],
KolmionKulmat[i,3],KolmionKulmat[i,4])[1]
}
AB <- Etaisyys(6371, tyovek[4],tyovek[5], tyovek[6], tyovek[7],1)
BC <- Etaisyys(6371, tyovek[6], tyovek[7], tyovek[8], tyovek[9],1)
CD <- Etaisyys(6371, tyovek[8], tyovek[9], tyovek[10], tyovek[11],1)
DA <- Etaisyys(6371, tyovek[10], tyovek[11], tyovek[4], tyovek[5],1)
AC <- Etaisyys(6371, tyovek[4], tyovek[5], tyovek[8], tyovek[9],1)
tulos <- NULL
tulos[1] <- sum(ala)
tulos[2] <- HeronAla(AB, BC, AC) + HeronAla(AC, CD, DA)
tulos[3] <- tulos[1]/tulos[2]
tulos
}
#----------------------------------------------------------------------------
#kv on vektori, joka sisaltaa kolmion kulmat
KolmioYmpyranAla <- function(kp1, kp2, r, kv1,kv2,kv3,kv4,kv5,kv6)
{
ala <- NULL
tilavek <- KolmionSisalla(kp1,kp2,kv1,kv2,kv3,kv4,kv5,kv6)
#oikea kolmioa ja piste kolmion sisalla
if(tilavek[1]==0 & tilavek[8]==1)
{
ala <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r, kv1,kv2,kv3,kv4,kv5,kv6)
}
#Oikea kolmio ja piste kolmion ulkopuolella
if(tilavek[1]==0 & tilavek[8]==0)
{
paikkavek <- PisteenAlue(kp1,kp2,tilavek[2:7])
A <- tilavek[2:3]
B <- tilavek[4:5]
C <- tilavek[6:7]
AB <- Etaisyys(6371,A[1],A[2],B[1],B[2],1)
BC <- Etaisyys(6371,B[1],B[2],C[1],C[2],1)
AC <- Etaisyys(6371,A[1],A[2],C[1],C[2],1)
if(paikkavek[1] == 1)
{
# kp = (1 - paikkavek[3])*C + paikkavek[3]*(paikkavek[2]*B + (1-paikkavek[2])*A)
uusiA <- (1 - (paikkavek[3] + 1))*C + (paikkavek[3] + 1)*A
uusiB <- (1 - (paikkavek[3] + 1))*C + (paikkavek[3] + 1)*B
uusiC <- C
#Pisteiden koordinaatit syötettävä kiertosuunnassa vastapäivään!!!!!!!!!!
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r,
uusiB[1], uusiB[2], uusiA[1], uusiA[2], uusiC[1], uusiC[2])[1]-
YmpyraNelikulmioKpSisalla(kp1, kp2, r,
uusiB[1], uusiB[2], uusiA[1], uusiA[2], A[1], A[2], B[1], B[2])[1]
}
if(paikkavek[1] == 2)
{
#iso kolmio - pikkukolmiot
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], C[1],C[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, B[1], B[2], C[1],C[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], B[1],B[2], kp1, kp2)[1]
}
if(paikkavek[1] == 3)
{
uusiA <- A
uusiB <- (1 - (paikkavek[3] + 1))*A + (paikkavek[3] + 1)*B
uusiC <- (1 - (paikkavek[3] + 1))*A + (paikkavek[3] + 1)*C
#Pisteiden koordinaatit syötettävä kiertosuunnassa vastapäivään!!!!!!!!!!
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r,
uusiB[1], uusiB[2], uusiA[1], uusiA[2], uusiC[1], uusiC[2])[1]-
YmpyraNelikulmioKpSisalla(kp1, kp2, r,
uusiB[1], uusiB[2], B[1], B[2], C[1], C[2], uusiC[1], uusiC[2])[1]
}
if(paikkavek[1] == 4)
{
#iso kolmio - pikkukolmiot
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], B[1],B[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, B[1], B[2], C[1],C[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], C[1],C[2], kp1, kp2)[1]
}
if(paikkavek[1] == 5)
{
uusiA <- (1 - (paikkavek[3] + 1))*B + (paikkavek[3] + 1)*A
uusiB <- B
uusiC <- (1 - (paikkavek[3] + 1))*B + (paikkavek[3] + 1)*C
#Pisteiden koordinaatit syötettävä kiertosuunnassa vastapäivään!!!!!!!!!!
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r,
uusiB[1], uusiB[2], uusiA[1], uusiA[2], uusiC[1], uusiC[2])[1]-
YmpyraNelikulmioKpSisalla(kp1, kp2, r,
uusiA[1], uusiA[2], uusiC[1], uusiC[2], C[1], C[2], A[1], A[2])[1]
}
if(paikkavek[1] == 6)
{
#iso kolmio - pikkukolmiot
ala[1] <- AlaYmpyranKpKolmionSisalla(kp1, kp2, r, C[1], C[2], B[1],B[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, B[1], B[2], A[1],A[2], kp1, kp2)[1]-
AlaYmpyranKpKolmionSisalla(kp1, kp2, r, A[1], A[2], C[1],C[2], kp1, kp2)[1]
}
ala[2] <- HeronAla(AB,BC,AC)
ala[3] <- ala[1]/ala[2]
}
ala
}
#----------------------------------------------------------------------------
# Syötä kulmat kiertojärjestyksessä vastapäivään!!!!!!
YmpyraNelikulmioAla <-function(kp,r,kulmat)
{
tulos <- NULL
apu1 <- KolmioYmpyranAla(kp[1], kp[2], r, kulmat[1],kulmat[2],kulmat[3],kulmat[4],kulmat[7],kulmat[8])
apu2 <- KolmioYmpyranAla(kp[1], kp[2], r, kulmat[5],kulmat[6],kulmat[3],kulmat[4],kulmat[7],kulmat[8])
tulos[1:2] <- apu1[1:2] + apu2[1:2]
tulos[3] <- tulos[1]/tulos[2]
tulos
}
|
- Code from U:\arkisto_kuopio\huippuyksikko\Tutkimus\R19_kopra\Mallit\Mika\Concentration_SOX (27.2.2006)
Piirto.r
Polku <- "C:\\omat\\dataa\\Concentration_SOX\\"
PaastoLahde <- "SO2"
postscript(file=paste(Polku,"IFConcentration",PaastoLahde,".ps",sep=""), width=10,height=10,horizontal=F)
IFmat10<-read.table(paste(Polku,"IFConcentration",PaastoLahde,"10km.txt",sep=""))
IFmat50<-read.table(paste(Polku,"IFConcentration",PaastoLahde,"50km.txt",sep=""))
#Ensin piirretään pisteen High sitten Low (=rgb(0,0,0) = musta)
for(i in 1:9){
piste <- paste("Piste ", i, sep="")
plot(c(1:12),IFmat10[,(2*i - 1)],type="b", xlim=c(1,12), ylim=c(0,1*10^(-5)),
xlab="kuukausi", ylab="IF 10km", col=rgb(1,0,0), main = paste(piste,"High = punainen, Low = musta"))
par(new=T)
plot(c(1:12),IFmat10[,(2*i)],type="b", xlim=c(1,12), ylim=c(0,1*10^(-5)),
xlab="kuukausi", ylab="IF 10km", col=rgb(0,0,0), main = paste(piste,"High = punainen, Low = musta"))
}
for(i in 1:9){
piste <- paste("Piste ", i, sep="")
plot(c(1:12),IFmat50[,(2*i - 1)],type="b", xlim=c(1,12), ylim=c(0,1*10^(-5)),
xlab="kuukausi", ylab="IF 50km", col=rgb(1,0,0), main = paste(piste,"High = punainen, Low = musta"))
par(new=T)
plot(c(1:12),IFmat50[,(2*i)],type="b", xlim=c(1,12), ylim=c(0,1*10^(-5)),
xlab="kuukausi", ylab="IF 50km", col=rgb(0,0,0), main = paste(piste,"High = punainen, Low = musta"))
}
# musta on 0,0,0
dev.off()
#_-----------------------------------------------------------
Polku <- "C:\\omat\\dataa\\Concentration_SOX\\"
PaastoLahde <- "SO4"
postscript(file=paste(Polku,"IFConcentration",PaastoLahde,".ps",sep=""), width=10,height=10,horizontal=F)
IFmat10<-read.table(paste(Polku,"IFConcentration",PaastoLahde,"10km.txt",sep=""))
IFmat50<-read.table(paste(Polku,"IFConcentration",PaastoLahde,"50km.txt",sep=""))
#Ensin piirretään pisteen High sitten Low (=rgb(0,0,0) = musta)
for(i in 1:9){
piste <- paste("Piste ", i, sep="")
plot(c(1:12),IFmat10[,(2*i - 1)],type="b", xlim=c(1,12), ylim=c(0,1.5*10^(-4)),
xlab="kuukausi", ylab="IF 10km", col=rgb(1,0,0), main = paste(piste,"High = punainen, Low = musta"))
par(new=T)
plot(c(1:12),IFmat10[,(2*i)],type="b", xlim=c(1,12), ylim=c(0,1.5*10^(-4)),
xlab="kuukausi", ylab="IF 10km", col=rgb(0,0,0), main = paste(piste,"High = punainen, Low = musta"))
}
for(i in 1:9){
piste <- paste("Piste ", i, sep="")
plot(c(1:12),IFmat50[,(2*i - 1)],type="b", xlim=c(1,12), ylim=c(0,1.5*10^(-4)),
xlab="kuukausi", ylab="IF 50km", col=rgb(1,0,0), main = paste(piste,"High = punainen, Low = musta"))
par(new=T)
plot(c(1:12),IFmat50[,(2*i)],type="b", xlim=c(1,12), ylim=c(0,1.5*10^(-4)),
xlab="kuukausi", ylab="IF 50km", col=rgb(0,0,0), main = paste(piste,"High = punainen, Low = musta"))
}
# musta on 0,0,0
dev.off()
|
koodi.r
# Datat -> N:\Huippuyksikko\Tutkimus\R19_Kopra\Data\Model_Converted
#--------------
# Concentration_SOX
# 50x50 ruudut
AlkuAika <- Sys.time()
SummaPixCi50kmSO2 <- matrix(nrow=12,ncol=18)
SummaPixCi50kmSO4 <- matrix(nrow=12,ncol=18)
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "c:\\omat\\dataa\\concentration_SOX\\"
PopulaatioData <- NULL
PopulaatioData<-read.table(paste(polku,"Population_Europe\\50km\\ciesin_center_63lon_50km_2",sep=""))
for(i in 1:12){
tiedosto1<-paste("joinfilelist_",i,"_SO2.joined.out",sep="")
tiedosto2<-paste("joinfilelist_",i,"_SO4.joined.out",sep="")
PitoisuusDataSO2 <- NULL
PitoisuusDataSO2 <-read.table(paste(polku,"Concentration_SOX\\",tiedosto1,sep=""))
PitoisuusDataSO4 <- NULL
PitoisuusDataSO4 <-read.table(paste(polku,"Concentration_SOX\\",tiedosto2,sep=""))
for(j in 1:18){
SummaPixCi50kmSO2[i,j] <- sum( PitoisuusDataSO2[,(j+2)] * PopulaatioData[,3] )
SummaPixCi50kmSO4[i,j] <- sum( PitoisuusDataSO4[,(j+2)] * PopulaatioData[,3] )
}
}
# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
BR <- 20/(24 * 60 * 60)
# Q <- 1 * 1/(365 * 24 * 60 * 60) mol/s
Q <- 1/(365 * 24 * 60 * 60)
IFConcentrationPoint50kmSO2 <- SummaPixCi50kmSO2 * BR/Q
IFConcentrationPoint50kmSO4 <- SummaPixCi50kmSO4 * BR/Q
TalletusKohde <- paste(talletuspolku,"IFConcentrationSO250km.txt",sep="")
write.table(IFConcentrationPoint50kmSO2, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
TalletusKohde <- paste(talletuspolku,"IFConcentrationSO450km.txt",sep="")
write.table(IFConcentrationPoint50kmSO4, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
LoppuAika <- Sys.time()
LoppuAika - AlkuAika
rm(SummaPixCi50kmSO2, SummaPixCi50kmSO4,polku, talletuspolku, PopulaatioData, AlkuAika, LoppuAika, tiedosto1, tiedosto2, PitoisuusDataSO2, PitoisuusDataSO4, Q,
BR, IFConcentrationPoint50kmSO2, IFConcentrationPoint50kmSO4, TalletusKohde)
#-----------------------------
# 10x10 ruudut
AlkuAika <- Sys.time()
SummaPixCi10kmSO2 <- matrix(nrow=12,ncol=18)
SummaPixCi10kmSO4 <- matrix(nrow=12,ncol=18)
polku <- "N:\\Huippuyksikko\\Tutkimus\\R19_Kopra\\Data\\Model_Converted\\"
talletuspolku <- "c:\\omat\\dataa\\concentration_SOX\\"
PopulaatioData <- NULL
PopulaatioData<-read.table(paste(polku,"Population_Europe\\10km\\ciesin_center_63lon_10km_2",sep=""))
for(i in 1:12){
tiedosto1<-paste("joinfilelist_",i,"_SO2.joined.out.1",sep="")
tiedosto2<-paste("joinfilelist_",i,"_SO4.joined.out.1",sep="")
PitoisuusDataSO2 <- NULL
PitoisuusDataSO2 <-read.table(paste(polku,"Concentration_SOX\\",tiedosto1,sep=""))
PitoisuusDataSO4 <- NULL
PitoisuusDataSO4 <-read.table(paste(polku,"Concentration_SOX\\",tiedosto2,sep=""))
for(j in 1:18){
SummaPixCi10kmSO2[i,j] <- sum( PitoisuusDataSO2[,(j+2)] * PopulaatioData[,3] )
SummaPixCi10kmSO4[i,j] <- sum( PitoisuusDataSO4[,(j+2)] * PopulaatioData[,3] )
}
}
# BR <- 20 m^3/päivä
# BR <- 20/(24 * 60 * 60) m^3/s
BR <- 20/(24 * 60 * 60)
# Q <- 1/(365 * 24 * 60 * 60) mol/s
Q <- 1/(365 * 24 * 60 * 60)
IFConcentrationPoint10kmSO2 <- SummaPixCi10kmSO2 * BR/Q
IFConcentrationPoint10kmSO4 <- SummaPixCi10kmSO4 * BR/Q
TalletusKohde <- paste(talletuspolku,"IFConcentrationSO210km.txt",sep="")
write.table(IFConcentrationPoint10kmSO2, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
TalletusKohde <- paste(talletuspolku,"IFConcentrationSO410km.txt",sep="")
write.table(IFConcentrationPoint10kmSO4, file = TalletusKohde, row.names = FALSE, col.names = FALSE)
LoppuAika <- Sys.time()
LoppuAika - AlkuAika
rm(SummaPixCi10kmSO2, SummaPixCi10kmSO4,polku, talletuspolku, PopulaatioData, AlkuAika, LoppuAika, tiedosto1, tiedosto2, PitoisuusDataSO2, PitoisuusDataSO4, Q,
BR, IFConcentrationPoint10kmSO2, IFConcentrationPoint10kmSO4, TalletusKohde)
|
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
#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)
|
Tainio_iF_emis_eur_pop_all_data-conbination.R
# 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)
|
Tainio_iF_emis_fin_pop_all_02.R
#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)
|
Tainio_iF_emis_fin_pop_all_data_conbination.R
# 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)
|
Tainio_population_01.R
#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)
|
Retrieved from "https://en.opasnet.org/index.php?title=Kopra&oldid=40405"