Difference between revisions of "Kopra"

From Testiwiki
Jump to: navigation, search
(started to document Kopra project's R code)
(→‎Calculations: Mika's codes added)
Line 16: Line 16:
  
 
=== Calculations ===
 
=== 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
 +
 +
<rcode>
 +
# 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)
 +
}
 +
</rcode>
 +
 +
'''IFmaille_erikseen.r
 +
 +
<rcode>
 +
# 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)
 +
}
 +
 +
#---------------------
 +
</rcode>
 +
 +
:''Code from U:\arkisto_kuopio\huippuyksikko\Tutkimus\R19_kopra\Mallit\Mika\Concentration_Finland (date 1.3.2006)
 +
 +
<rcode>
 +
# 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()
 +
</rcode>
 +
 +
:'' 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)
 +
 +
<rcode>
 +
# 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)
 +
</rcode>
 +
 +
'''Distance_piirto.r (25.1.2006)
 +
 +
<rcode>
 +
 +
#-----------------------------------------------------------------------
 +
#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)
 +
}
 +
</rcode>
 +
 +
'''Distance.r (17.2.2006)
 +
Codes varaDistance.r (31.1.2006) and varmuuskopioDistance.r (13.2.3006) omitted.
 +
 +
<rcode>
 +
 +
#-----------------------------------------------------------------------
 +
#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]
 +
</rcode>
 +
 +
'''datansiirtoYMS.r (24.2.2006)
 +
 +
<rcode>
 +
 +
#-----------------------------------------------------------------------
 +
#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]
 +
</rcode>
 +
 +
'''ympyran_ja-monikulmion_leikkaus.r (16.6.2006)
 +
 +
<rcode>
 +
#-----------------------------------------------------------------------
 +
#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
 +
}
 +
</rcode>
 +
 +
:''Code from U:\arkisto_kuopio\huippuyksikko\Tutkimus\R19_kopra\Mallit\Mika\Concentration_SOX (27.2.2006)
 +
 +
'''Piirto.r
 +
 +
<rcode>
 +
 +
 +
 +
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()
 +
</rcode>
 +
 +
'''koodi.r
 +
 +
<rcode>
 +
# 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)
 +
</rcode>
  
 
==== Marko's iF calculations ====
 
==== Marko's iF calculations ====

Revision as of 09:02, 19 April 2017


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

Research question

Answer

Rationale

Data

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

Calculations

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

+ Show code

IFmaille_erikseen.r

+ Show code

Code from U:\arkisto_kuopio\huippuyksikko\Tutkimus\R19_kopra\Mallit\Mika\Concentration_Finland (date 1.3.2006)

+ Show code

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)

+ Show code

Distance_piirto.r (25.1.2006)

+ Show code

Distance.r (17.2.2006) Codes varaDistance.r (31.1.2006) and varmuuskopioDistance.r (13.2.3006) omitted.

+ Show code

datansiirtoYMS.r (24.2.2006)

+ Show code

ympyran_ja-monikulmion_leikkaus.r (16.6.2006)

+ Show code

Code from U:\arkisto_kuopio\huippuyksikko\Tutkimus\R19_kopra\Mallit\Mika\Concentration_SOX (27.2.2006)

Piirto.r

+ Show code

koodi.r

+ Show code

Marko's iF calculations

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

Tainio_iF_emis_eur_pop_all_01.R

+ Show code

Tainio_iF_emis_eur_pop_all_data-conbination.R

+ Show code

Tainio_iF_emis_fin_pop_all_02.R

+ Show code

Tainio_iF_emis_fin_pop_all_data_conbination.R

+ Show code

Tainio_population_01.R

+ Show code

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