#spotkanie 08/11/2016
#saveRDS(calosc, file="data/polanka.rds")
library(openair)
## Loading required package: maps
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(lubridate)
##
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
##
## date
calosc <- readRDS( file="~/github/pm/data/polanka.rds")
ind <- which(year(calosc$date)>=2016)
calosc$PZ_WS[ind] <- calosc$PZ_WS[ind]*2
zima <- selectByDate(calosc, month =c(12,1,2))
timeAverage(zima, avg.time="season") %>% .[-1,] #%>% select(., date, PZ_Pol_pm10, PZ_Pol_pm25,PZ_Pol_t2m)
## # A tibble: 10 x 7
## date season PZ_Pol_pm10 PZ_Pol_pm25 PZ_T2M
## <time> <ord> <dbl> <dbl> <dbl>
## 1 2007-01-14 23:59:58 winter (DJF) 29.56613 NaN 3.47117075
## 2 2008-01-15 11:30:00 winter (DJF) 39.68802 33.81280 2.62083333
## 3 2009-01-14 23:59:58 winter (DJF) 56.05792 48.16566 -0.52827395
## 4 2010-01-14 23:59:58 winter (DJF) 60.98054 56.03106 -2.73304026
## 5 2011-01-14 23:59:58 winter (DJF) 55.49688 48.24546 -2.76992133
## 6 2012-01-15 11:30:00 winter (DJF) 44.98095 36.87825 0.07715201
## 7 2013-01-14 23:59:58 winter (DJF) 52.30306 44.34013 -1.14234151
## 8 2014-01-14 23:59:58 winter (DJF) 41.59804 29.17627 1.91813975
## 9 2015-01-15 00:59:50 winter (DJF) 40.06647 24.24244 1.72991216
## 10 2016-01-15 11:30:00 winter (DJF) 38.05686 30.93792 3.05819597
## # ... with 2 more variables: PZ_WS <dbl>, PZ_WD <dbl>
summaryPlot(calosc)
## date1 date2 PZ_Pol_pm10 PZ_Pol_pm25 PZ_T2M PZ_WS
## "POSIXct" "POSIXt" "integer" "integer" "numeric" "numeric"
## PZ_WD
## "numeric"
srednie_zima <- timeAverage(calosc, avg.time="day")
srednie_zima$yy <- lubridate::year(srednie_zima$date)
srednie_zima$mm <- lubridate::month(srednie_zima$date)
srednie_zima <- srednie_zima[which(srednie_zima$mm==11 | srednie_zima$mm==12 |
srednie_zima$mm==1 | srednie_zima$mm==2 |
srednie_zima$mm==3 ),] # grzewczy
srednie_zima$yy <- ifelse(srednie_zima$mm==12, srednie_zima$yy+1, srednie_zima$yy)
srednie_zima$yy <- ifelse(srednie_zima$mm==11, srednie_zima$yy+1, srednie_zima$yy)
ile50 <- function(x) sum(x>49, na.rm=T)
HDD <- function (x, x0 = 0, na.rm = TRUE) {
cold <- x < x0
hdd <- sum(x0 - x[cold], na.rm = na.rm)
return(hdd)
}
staty <- srednie_zima %>% dplyr::group_by(yy) %>%
dplyr::summarise( hdd0=round(HDD(PZ_T2M, 0)),
hdd10=round(HDD(PZ_T2M, 10)),
iledni_pow50=ile50(PZ_Pol_pm10), temp=round(mean(PZ_T2M),2), predkosc_wiatru=round(mean(PZ_WS),2)) %>% .[-c(1,12),]
# czy istnieja zwiazki statystyczne pomiedzy
cor(staty$temp, staty$iledni_pow50, method = "spearman")
## [1] -0.8389097
cor.test(staty$hdd10, staty$iledni_pow50)
##
## Pearson's product-moment correlation
##
## data: staty$hdd10 and staty$iledni_pow50
## t = 5.0779, df = 8, p-value = 0.0009556
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5422726 0.9698024
## sample estimates:
## cor
## 0.8736197
t.test(staty$hdd0, staty$iledni_pow50)
##
## Welch Two Sample t-test
##
## data: staty$hdd0 and staty$iledni_pow50
## t = 3.1715, df = 9.185, p-value = 0.01105
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 36.25888 214.74112
## sample estimates:
## mean of x mean of y
## 176.7 51.2
wilcox.test(staty$hdd10, staty$iledni_pow50)
## Warning in wilcox.test.default(staty$hdd10, staty$iledni_pow50): cannot
## compute exact p-value with ties
##
## Wilcoxon rank sum test with continuity correction
##
## data: staty$hdd10 and staty$iledni_pow50
## W = 100, p-value = 0.0001817
## alternative hypothesis: true location shift is not equal to 0
staty
## # A tibble: 10 x 6
## yy hdd0 hdd10 iledni_pow50 temp predkosc_wiatru
## <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 2007 47 796 38 4.76 4.52
## 2 2008 62 1068 34 2.98 4.13
## 3 2009 155 1273 56 1.60 3.80
## 4 2010 356 1458 65 0.42 3.74
## 5 2011 369 1520 62 -0.02 3.83
## 6 2012 224 1216 61 2.05 4.00
## 7 2013 292 1504 66 0.04 3.52
## 8 2014 99 968 47 3.64 3.70
## 9 2015 71 1013 49 3.30 3.85
## 10 2016 92 926 34 4.02 4.00
staty <- rbind(staty, colMeans(staty))
staty[11,1] <- 9999
model <- lm(staty$iledni_pow50~staty$hdd10+staty$temp+staty$predkosc_wiatru)
staty$pm10_theoret <- round(predict(model, staty),1)
staty
## # A tibble: 11 x 7
## yy hdd0 hdd10 iledni_pow50 temp predkosc_wiatru pm10_theoret
## * <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2007 47.0 796.0 38.0 4.760 4.520 34.0
## 2 2008 62.0 1068.0 34.0 2.980 4.130 46.0
## 3 2009 155.0 1273.0 56.0 1.600 3.800 56.2
## 4 2010 356.0 1458.0 65.0 0.420 3.740 62.2
## 5 2011 369.0 1520.0 62.0 -0.020 3.830 64.8
## 6 2012 224.0 1216.0 61.0 2.050 4.000 50.9
## 7 2013 292.0 1504.0 66.0 0.040 3.520 67.3
## 8 2014 99.0 968.0 47.0 3.640 3.700 45.0
## 9 2015 71.0 1013.0 49.0 3.300 3.850 46.8
## 10 2016 92.0 926.0 34.0 4.020 4.000 38.8
## 11 9999 176.7 1174.2 51.2 2.279 3.909 51.2
#summary(model)
plot(staty$yy[-11], staty$hdd0[-11], type='l', ylab='Heating Degree Days (0*C)', xlab="year", main="XI-III")
par(new=T)
plot(staty$yy[-11], staty$iledni_pow50[-11], type='l', xaxt='n', yaxt='n', ylab='', xlab='', col="blue")
axis(4)
mtext(side = 4, "liczba przekroczen progu 50uq/m3 ")
