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