library(maps); library(mapdata); library(maptools); library(chron); library(maptools)
##
## # ATTENTION: maps v3.0 has an updated 'world' map. #
## # Many country borders and names have changed since 1990. #
## # Type '?world' or 'news(package="maps")'. See README_v3. #
##
##
## Loading required package: sp
## Checking rgeos availability: TRUE
library(fields); library(lubridate); library(raster); library(rgdal); library(magic)
## Loading required package: spam
## Loading required package: grid
## Spam version 1.0-1 (2014-09-09) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attaching package: 'spam'
##
## Następujące obiekty zostały zakryte z 'package:base':
##
## backsolve, forwardsolve
##
##
## Attaching package: 'fields'
##
## Następujący obiekt został zakryty z 'package:maps':
##
## ozone
##
##
## Attaching package: 'lubridate'
##
## Następujące obiekty zostały zakryte z 'package:chron':
##
## days, hours, minutes, seconds, years
##
## rgdal: version: 1.0-7, (SVN revision 559)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
## Path to GDAL shared files: /usr/share/gdal/1.11
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: (autodetected)
## Linking to sp version: 1.2-0
## Loading required package: abind
##
## Attaching package: 'magic'
##
## Następujący obiekt został zakryty z 'package:raster':
##
## shift
library(modiscloud); library(openair); library(plyr); library(dplyr)
## Loading required package: date
## Loading required package: sfsmisc
## Loading required package: lazyeval
## Loading required package: dplyr
##
## Attaching package: 'dplyr'
##
## Następujący obiekt został zakryty z 'package:sfsmisc':
##
## last
##
## Następujące obiekty zostały zakryte z 'package:raster':
##
## intersect, select, union
##
## Następujące obiekty zostały zakryte z 'package:lubridate':
##
## intersect, setdiff, union
##
## Następujące obiekty zostały zakryte z 'package:stats':
##
## filter, lag
##
## Następujące obiekty zostały zakryte z 'package:base':
##
## intersect, setdiff, setequal, union
##
## -------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## -------------------------------------------------------------------------
##
## Attaching package: 'plyr'
##
## Następujące obiekty zostały zakryte z 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
##
## Następujący obiekt został zakryty z 'package:lubridate':
##
## here
##
## Następujący obiekt został zakryty z 'package:fields':
##
## ozone
##
## Następujący obiekt został zakryty z 'package:maps':
##
## ozone
#punkty_kasi <- readOGR(dsn="/home/bartosz/Pulpit/kasia/", layer="feno_1_1_centroidy_wgs84")
#pl=readOGR(dsn="/home/bartosz/Dokumenty/projekty/feno2012/gis/shp/", layer="atpol_duze")
#pl_woj=readOGR(dsn="/home/bartosz/Dokumenty/projekty/isok/shp", layer="POL_adm1")
#atpol=readOGR(dsn="/home/bartosz/Dokumenty/publikacje własne/2015/aerobiologia/gis/", layer="atpol_3_3_wgs84")
#rzeki=readOGR(dsn="/home/bartosz/Dokumenty/projekty/isok/shp", layer="rzeki")
#pl_granice=readOGR(dsn="/home/bartosz/Dokumenty/projekty/feno2012/gis/shp/", layer="pl_granice")
#
# #setwd("/home/bartosz/Pulpit/modis/przyciete/aqua/")
# setwd("/media/bartosz/toshiba/przyciete")
# #setwd("/media/bartosz/solarcity/modis/przyciete/")
#
# wynik <- NULL
# pliki=dir(pattern=".tif")
# daty <- unique(substr(pliki,1,8))
#
#
# for (i in 51:length(daty)){
#
# evi_tmp <- paste0(daty[i],".250m_16_days_EVI.tif")
# ndvi_tmp <- paste0(daty[i],".250m_16_days_NDVI.tif")
# doy_tmp <- paste0(daty[i],".250m_16_days_composite_day_of_the_year.tif")
# rel_tmp <- paste0(daty[i],".250m_16_days_pixel_reliability.tif")
#
#
# evi_raster <- raster(evi_tmp)
# evi_raster[which(evi_raster[]==-3000)]=NA
# ndvi_raster <- raster(ndvi_tmp)
# ndvi_raster[which(ndvi_raster[]==-3000)]=NA
# doy_raster <- raster(doy_tmp)
# doy_raster[which(doy_raster[]==-1)]=NA
# rel_raster <- raster(rel_tmp)
# rel_raster[which(rel_raster[]==255)]=NA
#
#
# rok <- as.numeric(substr(daty[i],1,4))
# miesiac <- as.numeric(substr(daty[i],5,6))
#
# tmp <- cbind(as.data.frame(punkty_kasi)[,2:3],extract(ndvi_raster,punkty_kasi), extract(evi_raster,punkty_kasi), extract(doy_raster,punkty_kasi), rok, extract(rel_raster,punkty_kasi))
# names(tmp) <- c("feno","meteo","ndvi","evi","doy","year", "reliability")
#
# if(miesiac==12) tmp$doy <- ifelse(tmp$doy<300,tmp$doy+365,tmp$doy)
#
#
#
# wynik=rbind(wynik,tmp)
# print(i)
#
# if(i%%50==0) saveRDS(wynik,file="/home/bartosz/Dokumenty/publikacje własne/2015/aerobiologia/NDVI_kasia.rds", compress=T)
# }
# saveRDS(wynik,file="/home/bartosz/Dokumenty/publikacje własne/2015/aerobiologia/NDVI_kasia.rds", compress=T)
wynik <- readRDS(file="/home/bartosz/Dokumenty/publikacje własne/2015/aerobiologia/NDVI_kasia.rds")
library(dplyr)
library(hydroTSM)
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## Następujące obiekty zostały zakryte z 'package:base':
##
## as.Date, as.Date.numeric
##
## Loading required package: xts
##
## Attaching package: 'xts'
##
## Następujące obiekty zostały zakryte z 'package:dplyr':
##
## first, last
##
## Następujący obiekt został zakryty z 'package:sfsmisc':
##
## last
##
##
## Attaching package: 'hydroTSM'
##
## Następujący obiekt został zakryty z 'package:raster':
##
## extract
wynik$data <- as.Date(paste0(wynik$year,"-01-01"))+wynik$doy-1
wynik$ndvi <- ifelse(wynik$reliability<=2,yes = wynik$ndvi,no=NA)
wynik$evi <- ifelse(wynik$reliability<=2,yes = wynik$evi,no=NA)
time_full=as.data.frame(seq(from=as.Date("2006-12-01"),to=as.Date("2015-01-31"),by="day"))
names(time_full) <- "data"
time <- seq(from=as.Date("2006-12-01"),to=as.Date("2015-01-31"),by="day")
calosc <- NULL
for (i in as.character(unique(wynik$feno))) {
tmp <- subset(wynik,feno==i)
tmp <- join(time_full,tmp,by="data")
stacja <- as.character(na.omit(unique(tmp$meteo))[1])
proste_ndvi <- as.data.frame(approx(tmp$data,tmp$ndvi, xout=time))
splajn_ndvi <- as.data.frame(spline(tmp$data,tmp$ndvi, xout=time,method = "natural"))
proste_evi <- as.data.frame(approx(tmp$data,tmp$evi, xout=time))
splajn_evi <- as.data.frame(spline(tmp$data,tmp$evi, xout=time,method = "natural"))
tmp <- cbind.data.frame(proste_ndvi,splajn_ndvi[,2], proste_evi[,2],splajn_evi[,2])
tmp <- subset(tmp,x>=as.Date("2007-01-01") & x<=as.Date("2014-12-31"))
names(tmp)=c("date","ndvi_simpl_1x1","ndvi_splajn_1x1","evi_simpl_1x1","evi_splajn_1x1")
# myplot <- timeVariation(tmp,pollutant = c("ndvi_splajn_1x1","ndvi_simpl_1x1"), normalise=F, labels=i)
# plot(myplot, subset="month", xlab="miesiac", legend=F)
#print(i)
#barplot(monthlyfunction(tmp[,c(1,2)], mean)/10000, main=i)
par(mfrow=c(1,2))
x <- zoo(tmp[,2],tmp[,1])
y <- zoo(tmp[,4],tmp[,1])
plot(aggregate(x, as.yearmon, mean)/10000, xlab='',ylab='',main=stacja, ylim=c(0.1,0.9)); lines(aggregate(y, as.yearmon, mean)/10000,col="red")
plot((aggregate(x, as.Date(as.Date("2012-01-01")+as.numeric(strftime(x, format="%j"))), mean))[-366] /10000, ylim=c(0.1000,0.8900), main=i, lwd=2, ylab="NDVI, EVI", xlim=c(as.Date("2012-01-14"),as.Date("2012-12-19")))
lines((aggregate(y, as.Date(as.Date("2012-01-01")+as.numeric(strftime(x, format="%j"))), mean))[-366] / 10000, col="red", lwd=2)
abline(v=c(seq.Date(from=as.Date("2012-01-01"), to=as.Date("2012-12-31"), by="month")), lty=3)
abline(h=c(1:5*0.2000), lty=3)
tmp$fenostacja=i
tmp$meteostacja=stacja
calosc <- rbind.data.frame(calosc,tmp)
}
# testy z agregacja i plotowaniem
#plot(aggregate(x, week, mean)); lines(aggregate(y, week, mean),col="red") # wykres tygodniowy
#axis(1, at=1:53, labels=as.Date("2000-01-01")+time(aggregate(x, week, mean))*7-4)
#plot(aggregate(x, as.Date(paste(substr(time(x),1,4),"01-01",sep="-")), mean)/10000, ylab='', main=i, ylim=c(0.2,0.9))
#lines(aggregate(y, as.Date(paste(substr(time(x),1,4),"01-01",sep="-")), mean)/10000, col='red')
#
#
#
# tmp2 <- join(time_full,tmp,by="data")
#
# bialy <- subset(wynik,meteo=="Białystok")
# bialy2 <- subset(wynik,meteo=="Białystok")
# bialy$ndvi <- ifelse(bialy$reliability<=2,yes = bialy$ndvi,no=NA)
# bialy$evi <- ifelse(bialy$reliability<=2,yes = bialy$evi,no=NA)
# plot(bialy$data,bialy$ndvi/10000, type='l',col='red')
# lines(bialy2$data,bialy$ndvi/10000, type='p',col='black',cex=0.4)
# #lines(bialy$data,bialy$evi/10000, type='l', col="red")
#
# library(ggplot2)
# qplot(bialy$data,bialy$evi/10000, type='l', col="red")
#
# feno <- readRDS("~/Pulpit/feno.rds")
#